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