Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2017-2019 The plumed team
3 : (see the PEOPLE file at the root of the distribution for a list of names)
4 :
5 : See http://www.plumed.org for more information.
6 :
7 : This file is part of plumed, version 2.
8 :
9 : plumed is free software: you can redistribute it and/or modify
10 : it under the terms of the GNU Lesser General Public License as published by
11 : the Free Software Foundation, either version 3 of the License, or
12 : (at your option) any later version.
13 :
14 : plumed is distributed in the hope that it will be useful,
15 : but WITHOUT ANY WARRANTY; without even the implied warranty of
16 : MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
17 : GNU Lesser General Public License for more details.
18 :
19 : You should have received a copy of the GNU Lesser General Public License
20 : along with plumed. If not, see <http://www.gnu.org/licenses/>.
21 : +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
22 :
23 :
24 : #include "Colvar.h"
25 : #include "ActionRegister.h"
26 : #include "core/PlumedMain.h"
27 :
28 : #include <string>
29 : #include <cmath>
30 : #include <cassert>
31 : #include <iostream>
32 : #include <vector>
33 :
34 : using namespace std;
35 :
36 : namespace PLMD {
37 : namespace colvar {
38 :
39 : //+PLUMEDOC COLVAR DIMER
40 : /*
41 : This CV computes the Dimer interaction energy for a collection of Dimers.
42 :
43 : Each Dimer represents an atom, as described in the Dimer paper,
44 : JCTC 13, 425 (2017). A system of N atoms is thus represented with N Dimers, each
45 : Dimer being composed of two beads and eventually a virtual site representing its center of mass.
46 :
47 : A typical configuration for a dimerized system has the following ordering of atoms:
48 :
49 : 1 TAG1 X Y Z N atoms representing the first bead of each Dimer
50 :
51 : 2 TAG2 X Y Z
52 :
53 : ...
54 :
55 : N TAGN X Y Z N atoms representing the second bead of each Dimer
56 :
57 : N+1 TAG1 X Y Z
58 :
59 : N+2 TAG2 X Y Z
60 :
61 : ...
62 :
63 : 2N TAGN X Y Z Optional: N atoms representing the center of mass of each Dimer
64 :
65 : 2N+1 TAG1 X Y Z
66 :
67 : 2N+2 TAG2 X Y Z
68 :
69 : ...
70 :
71 : 3N TAGN X Y Z The configuration might go on with un-dimerized atoms (like a solvent)
72 :
73 : 3N+1
74 :
75 : 3N+2
76 :
77 : ...
78 :
79 :
80 : The Dimer interaction energy is defined between atoms x and N+x, for x=1,...,N and is
81 : characterized by two parameters Q and DSIGMA. These are passed as mandatory arguments along with
82 : the temperature of the system.
83 :
84 : \par Examples
85 :
86 : This line tells Plumed to compute the Dimer interaction energy for every dimer in the system.
87 :
88 : \plumedfile
89 : dim: DIMER TEMP=300 Q=0.5 ALLATOMS DSIGMA=0.002
90 : \endplumedfile
91 :
92 : If the simulation doesn't use virtual sites for the dimers centers of mass,
93 : Plumed has to know in order to determine correctly the total number of dimers from
94 : the total number of atoms:
95 : \plumedfile
96 : dim: DIMER TEMP=300 Q=0.5 ALLATOMS DSIGMA=0.002 NOVSITES
97 : \endplumedfile
98 :
99 : The NOVSITES flag is not required if one provides the atom serials of each Dimer. These are
100 : defined through two atomlists provided __instead__ of the ALLATOMS keyword.
101 : For example, the Dimer interaction energy of dimers specified by beads (1;23),(5;27),(7;29) is:
102 : \plumedfile
103 : dim: DIMER TEMP=300 Q=0.5 ATOMS1=1,5,7 ATOMS2=23,27,29 DSIGMA=0.002
104 : \endplumedfile
105 :
106 : Note that the ATOMS1,ATOMS2 keywords can support atom groups and
107 : interval notation as defined in \ref GROUP.
108 :
109 :
110 : In a Replica Exchange simulation the keyword DSIGMA can be used in two ways:
111 : if a plumed.n.dat file is provided for each replica, then DSIGMA is passed as a single value,
112 : like in the previous examples, and each replica will read its own DSIGMA value. If
113 : a unique plumed.dat is given, DSIGMA has to be a list containing a value for each replica.
114 : For 4 replicas:
115 : \plumedfile
116 : dim: DIMER TEMP=300 Q=0.5 ATOMS1=1,5,7 ATOMS2=23,27,29 DSIGMA=0.002,0.002,0.004,0.01
117 : \endplumedfile
118 :
119 :
120 : \par Usage of the CV
121 :
122 : The dimer interaction is not coded in the driver program and has to be inserted
123 : in the hamiltonian of the system as a linear RESTRAINT (see \ref RESTRAINT):
124 : \plumedfile
125 : dim: DIMER TEMP=300 Q=0.5 ALLATOMS DSIGMA=0.002
126 : RESTRAINT ARG=dim AT=0 KAPPA=0 SLOPE=1 LABEL=dimforces
127 : \endplumedfile
128 :
129 : In a replica exchange, Metadynamics (see \ref METAD) can be used on the Dimer CV to reduce
130 : the number of replicas. Just keep in mind that METAD SIGMA values should be tuned
131 : in the standard way for each replica according to the value of DSIGMA.
132 : */
133 : //+ENDPLUMEDOC
134 :
135 6 : class Dimer : public Colvar {
136 : public:
137 : static void registerKeywords( Keywords& keys);
138 : explicit Dimer(const ActionOptions&);
139 : virtual void calculate();
140 : protected:
141 : bool trimer,useall;
142 : int myrank, nranks;
143 : double qexp,temperature,beta,dsigma;
144 : vector<double> dsigmas;
145 : private:
146 : void consistencyCheck();
147 : vector<AtomNumber> usedatoms1;
148 : vector<AtomNumber> usedatoms2;
149 :
150 : };
151 :
152 6454 : PLUMED_REGISTER_ACTION(Dimer, "DIMER")
153 :
154 :
155 :
156 3 : void Dimer::registerKeywords( Keywords& keys) {
157 3 : Colvar::registerKeywords(keys);
158 :
159 12 : keys.add("compulsory","DSIGMA","The interaction strength of the dimer bond.");
160 12 : keys.add("compulsory", "Q", "The exponent of the dimer potential.");
161 12 : keys.add("compulsory", "TEMP", "The temperature (in Kelvin) of the simulation.");
162 12 : keys.add("atoms", "ATOMS1", "The list of atoms representing the first bead of each Dimer being considered by this CV. Used if ALLATOMS flag is missing");
163 12 : keys.add("atoms", "ATOMS2", "The list of atoms representing the second bead of each Dimer being considered by this CV. Used if ALLATOMS flag is missing");
164 9 : keys.addFlag("ALLATOMS", false, "Use EVERY atom of the system. Overrides ATOMS keyword.");
165 9 : keys.addFlag("NOVSITES", false, "If present the configuration is without virtual sites at the centroids.");
166 :
167 3 : }
168 :
169 :
170 :
171 2 : Dimer::Dimer(const ActionOptions& ao):
172 6 : PLUMED_COLVAR_INIT(ao)
173 : {
174 :
175 6 : log<<" Bibliography "<<plumed.cite("M Nava, F. Palazzesi, C. Perego and M. Parrinello, J. Chem. Theory Comput. 13, 425(2017)")<<"\n";
176 4 : parseVector("DSIGMA",dsigmas);
177 4 : parse("Q",qexp);
178 4 : parse("TEMP",temperature);
179 :
180 :
181 : vector<AtomNumber> atoms;
182 4 : parseFlag("ALLATOMS",useall);
183 2 : trimer=true;
184 : bool notrim;
185 4 : parseFlag("NOVSITES",notrim);
186 2 : trimer=!notrim;
187 :
188 2 : nranks=multi_sim_comm.Get_size();
189 2 : myrank=multi_sim_comm.Get_rank();
190 2 : if(dsigmas.size()==1)
191 2 : dsigma=dsigmas[0];
192 : else
193 0 : dsigma=dsigmas[myrank];
194 :
195 :
196 :
197 :
198 2 : if(useall)
199 : {
200 : // go with every atom in the system but not the virtuals...
201 : int natoms;
202 1 : if(trimer)
203 1 : natoms= 2*getTotAtoms()/3;
204 : else
205 0 : natoms=getTotAtoms()/2;
206 :
207 89 : for(unsigned int i=0; i<((unsigned int)natoms); i++)
208 : {
209 : AtomNumber ati;
210 : ati.setIndex(i);
211 44 : atoms.push_back(ati);
212 : }
213 : }
214 : else // serials for the first beads of each dimer are given
215 : {
216 2 : parseAtomList("ATOMS1",usedatoms1);
217 2 : parseAtomList("ATOMS2",usedatoms2);
218 :
219 : int isz1 = usedatoms1.size();
220 :
221 9 : for(unsigned int i=0; i<isz1; i++)
222 : {
223 : AtomNumber ati;
224 4 : ati.setIndex(usedatoms1[i].index());
225 4 : atoms.push_back(ati);
226 : }
227 :
228 : int isz2 = usedatoms2.size();
229 9 : for(unsigned int i=0; i<isz2; i++)
230 : {
231 : AtomNumber atip2;
232 4 : atip2.setIndex(usedatoms2[i].index());
233 4 : atoms.push_back(atip2);
234 : }
235 :
236 : }
237 2 : consistencyCheck();
238 2 : checkRead();
239 2 : beta = 1./(kBoltzmann*temperature);
240 :
241 2 : addValueWithDerivatives(); // allocate
242 2 : requestAtoms(atoms);
243 2 : setNotPeriodic();
244 2 : }
245 :
246 4 : void Dimer::calculate()
247 : {
248 4 : double cv_val=0;
249 4 : Tensor virial;
250 : vector<Vector> derivatives;
251 4 : vector<Vector> my_pos=getPositions();
252 4 : int atms = my_pos.size();
253 : vector<Vector> der_b2;
254 72 : for(int i=0; i<atms/2; i++)
255 : {
256 34 : Vector dist;
257 102 : dist = pbcDistance(my_pos[i],my_pos[i+atms/2]);
258 : double distquad=0;
259 238 : for(int j=0; j<3; j++)
260 102 : distquad += dist[j]*dist[j];
261 :
262 34 : double dsigquad = dsigma*dsigma;
263 34 : double fac1 = 1.0 + distquad/(2*qexp*dsigquad);
264 34 : double fac1qm1 = pow(fac1,qexp-1);
265 :
266 :
267 34 : cv_val += (fac1*fac1qm1-1.0)/beta;
268 34 : Vector der_val;
269 34 : Vector mder_val;
270 238 : for(int j=0; j<3; j++)
271 : {
272 102 : der_val[j] = -fac1qm1*dist[j]/(dsigquad*beta);
273 102 : mder_val[j]=-der_val[j];
274 : }
275 34 : derivatives.push_back(der_val);
276 34 : der_b2.push_back(mder_val);
277 :
278 : // virial part: each dimer contributes -x_{ij}*ds/dx_{ij} (s is the CV)
279 34 : double dfunc = fac1qm1/(beta*dsigquad);
280 34 : Vector dd(dfunc*dist);
281 34 : Tensor vv(dd,dist);
282 34 : virial -= vv;
283 :
284 : }
285 :
286 4 : derivatives.insert(derivatives.end(), der_b2.begin(), der_b2.end());
287 :
288 212 : for(unsigned int i=0; i<derivatives.size(); i++)
289 136 : setAtomsDerivatives(i,derivatives[i]);
290 :
291 4 : setValue(cv_val);
292 4 : setBoxDerivatives(virial);
293 :
294 4 : }
295 :
296 :
297 :
298 : /*****************
299 : There are some conditions that a valid input should satisfy.
300 : These are checked here and PLUMED error handlers are (eventually) called.
301 : ******************/
302 2 : void Dimer::consistencyCheck()
303 : {
304 3 : if(!useall && usedatoms1.size()!=usedatoms2.size())
305 0 : error("The provided atom lists are of different sizes.");
306 :
307 2 : if(qexp<0.5 || qexp>1)
308 0 : warning("Dimer CV is meant to be used with q-exponents between 0.5 and 1. We are not responsible for any black hole. :-)");
309 :
310 2 : if(dsigma<0)
311 0 : error("Please use positive sigma values for the Dimer strength constant");
312 :
313 2 : if(temperature<0)
314 0 : error("Please, use a positive value for the temperature...");
315 :
316 : // if dsigmas has only one element means that either
317 : // you are using different plumed.x.dat files or a plumed.dat with a single replica
318 2 : if(dsigmas.size()!=nranks && dsigmas.size()!=1)
319 0 : error("Mismatch between provided sigmas and number of replicas");
320 :
321 2 : }
322 :
323 :
324 : }
325 4839 : }
326 :
|
