Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2012-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 "core/ActionShortcut.h"
23 : #include "core/ActionRegister.h"
24 : #include "MultiColvarShortcuts.h"
25 : #include <string>
26 : #include <cmath>
27 :
28 : //+PLUMEDOC MCOLVAR DISTANCES
29 : /*
30 : Calculate the distances between multiple piars of atoms
31 :
32 : __This shortcut action allows you to calculate function of the distribution of distances and reproduces the syntax in older PLUMED versions.
33 : We think that the syntax in newer versions is more flexible and easier to use and that you should try to avoid using this old syntax. If you look at the example inputs below you can
34 : see how the new syntax can be used to replace this old one. We would strongly encourage you to use the newer syntax.__
35 :
36 : The following input tells plumed to calculate the distances between atoms 3 and 5 and
37 : between atoms 1 and 2 and to print the shorter of these two distances.
38 :
39 : ```plumed
40 : d1: DISTANCES ATOMS1=3,5 ATOMS2=1,2 LOWEST
41 : PRINT ARG=d1.lowest
42 : ```
43 :
44 : The following input tells plumed to calculate the distances between atoms 3 and 5 and between atoms 1 and 2
45 : and then to calculate the number of these distances that are less than 0.1 nm. The number of distances
46 : less than 0.1nm is then printed to a file.
47 :
48 : ```plumed
49 : d1: DISTANCES ATOMS1=3,5 ATOMS2=1,2 LESS_THAN={RATIONAL R_0=0.1}
50 : PRINT ARG=d1.lessthan
51 : ```
52 :
53 : The following input tells plumed to calculate all 3 distances between atoms 1, 2 and 3 (i.e. the distances between atoms
54 : 1 and 2, atoms 1 and 3 and atoms 2 and 3). The average of these distances is then calculated in two ways.
55 :
56 : ```plumed
57 : d1: DISTANCES GROUP=1-3 MEAN SUM
58 : d1mean: CUSTOM ARG=d1.sum FUNC=x/3 PERIODIC=NO
59 : PRINT ARG=d1.mean,d1mean
60 : ```
61 :
62 : The following input tells plumed to calculate all the distances between the atoms in GROUPA and the atoms in GROUPB.
63 : In other words the distances between atoms 1 and 2 and the distance between atoms 1 and 3. The number of distances
64 : more than 0.1 is then printed to a file.
65 :
66 : ```plumed
67 : d1: DISTANCES GROUPA=1 GROUPB=2,3 MORE_THAN={RATIONAL R_0=0.1}
68 : PRINT ARG=d1.morethan
69 : ```
70 :
71 : ## Calculating minimum distances
72 :
73 : There are various way for calculating the minimum distance between two groups of atoms. The simplest one
74 : is the one we saw above; namely:
75 :
76 : ```plumed
77 : d1: DISTANCES GROUPA=1-10 GROUPB=11-20 LOWEST
78 : PRINT ARG=d1.lowest
79 : ```
80 :
81 : A disadvantage of this approach is that the derivatives of the resulting value are not continuous. To be clear,
82 : the value of the shortest distance is continuous when you use this approach so you will likely be fine if you use
83 : the quantity output by the command above as an input in a bias. However, if you would like to use a quantity with
84 : continuous derivatives for the minimum distance you can use the following input:
85 :
86 : ```plumed
87 : d1: DISTANCES GROUPA=1-10 GROUPB=11-20 MIN={BETA=500.}
88 : PRINT ARG=d1.min FILE=colvar STRIDE=10
89 : ```
90 :
91 : In order to ensure that the minimum value has continuous derivatives we use the following function:
92 :
93 : $$
94 : d_{min} = \frac{\beta}{ \log \sum_i \exp\left( \frac{\beta}{d_i} \right) }
95 : $$
96 :
97 : where $\beta$ is a user specified parameter.
98 :
99 : This input is used rather than a separate MINDIST colvar so that the same routine and the same input style can be
100 : used to calculate minimum coordination numbers (see [COORDINATIONNUMBER](COORDINATIONNUMBER.md)), minimum
101 : angles (see [ANGLES](ANGLES.md)) and many other variables. If you expand the input above you will notice that
102 : the DISTANCES command is a shortcut. By expanding this shortcut you can see how PLUMED calculates the minimum distance.
103 :
104 : An alternative implementation of the minimum distance is provided through the ALT_MIN keyword that is demonstrated below:
105 :
106 : ```plumed
107 : d1: DISTANCES GROUPA=1-10 GROUPB=11-20 ALT_MIN={BETA=500.}
108 : PRINT ARG=d1.altmin FILE=colvar STRIDE=10
109 : ```
110 :
111 : The minimum value output by this command also has continuous derivatives and is evaluated using the following function:
112 :
113 : $$
114 : d_{min} = -\frac{1}{\beta} \log \left( \sum_i \exp(-\beta d_i) \right)
115 : $$
116 :
117 : ## Calculating maximum distances
118 :
119 : There are two ways for calculating the maximum distance between two groups of atom. The simplest one is:
120 :
121 : ```plumed
122 : d1: DISTANCES GROUPA=1-10 GROUPB=11-20 HIGHEST
123 : PRINT ARG=d1.highest
124 : ```
125 :
126 : As with LOWEST above a disadvantage of this approach is that the derivatives of the resulting value are not continuous.
127 : If you would like to use a quantity with continuous derivatives for the maximum distance you can use the following input:
128 :
129 : ```plumed
130 : d1: DISTANCES GROUPA=1-10 GROUPB=11-20 MAX={BETA=500.}
131 : PRINT ARG=d1.max
132 : ```
133 :
134 : The maximum value output by this command is evaluated using the following function:
135 :
136 : $$
137 : d_{max} = \beta \log \left( \sum_i \exp\left( \frac{x}{\beta} \right) \right)
138 : $$
139 :
140 : ## Other parts of the old multicolvar syntax
141 :
142 : The ways of calculating the minimum and maximum distance described above are part of plumed 2's multicolvar functionality. In older
143 : versions of PLUMED these special actions
144 : allowed you to calculate multiple functions of a distribution of simple collective variables. As an example you
145 : can calculate the number of distances less than 1.0, the minimum distance, the number of distances more than
146 : 2.0, the number of distances between 1.0 and 2.0 and histogram that esimates the distribution of distances by using the following command:
147 :
148 : ```plumed
149 : d1: DISTANCES ...
150 : GROUPA=1-10 GROUPB=11-20
151 : LESS_THAN={RATIONAL R_0=1.0}
152 : MORE_THAN={RATIONAL R_0=2.0}
153 : BETWEEN={GAUSSIAN LOWER=1.0 UPPER=2.0}
154 : HISTOGRAM={GAUSSIAN LOWER=0.0 UPPER=4.0 NBINS=10}
155 : MIN={BETA=500.}
156 : ...
157 : PRINT ARG=d1.* FILE=colvar STRIDE=10
158 : ```
159 :
160 : A calculation performed this way is fast because the expensive part of the calculation - the calculation of all the distances - is only done once per step.
161 : If you expand the shortcut in the input above you can see how this result is achieved using a combination of many actions. We recommend you use these combinations of
162 : actions rather than the shortcut when using newer versions of PLUMED as you will have a more control over what you are doing. This newer input syntax is easier to
163 : understand especially when you start evaluating more complicated CVs.
164 :
165 :
166 : */
167 : //+ENDPLUMEDOC
168 :
169 : //+PLUMEDOC MCOLVAR XDISTANCES
170 : /*
171 : Calculate the x components of the vectors connecting one or many pairs of atoms.
172 :
173 : __This shortcut action allows you to calculate functions of the distribution of the x components of the vectors connecting pairs of atoms and
174 : reproduces the syntax in older PLUMED versions.
175 : If you look at the example inputs below you can
176 : see how the new syntax operates. We would strongly encourage you to use the newer syntax as it offers greater flexibility.__
177 :
178 : The following input tells plumed to calculate the x-component of the vector connecting atom 3 to atom 5 and
179 : the x-component of the vector connecting atom 1 to atom 2. The minimum of these two quantities is then
180 : printed
181 :
182 : ```plumed
183 : d1: XDISTANCES ATOMS1=3,5 ATOMS2=1,2 MIN={BETA=0.1}
184 : PRINT ARG=d1.min
185 : ```
186 :
187 : The following input tells plumed to calculate the x-component of the vector connecting atom 3 to atom 5 and
188 : the x-component of the vector connecting atom 1 to atom 2. The number of values that are
189 : less than 0.1nm is then printed to a file.
190 :
191 : ```plumed
192 : d1: XDISTANCES ATOMS1=3,5 ATOMS2=1,2 LESS_THAN={RATIONAL R_0=0.1}
193 : PRINT ARG=d1.lessthan
194 : ```
195 :
196 : The following input tells plumed to calculate the x-components of all the distinct vectors that can be created
197 : between atoms 1, 2 and 3 (i.e. the vectors between atoms 1 and 2, atoms 1 and 3 and atoms 2 and 3).
198 : The average of these quantities is then calculated.
199 :
200 : ```plumed
201 : d1: XDISTANCES GROUP=1-3 MEAN
202 : PRINT ARG=d1.mean
203 : ```
204 :
205 : The following input tells plumed to calculate all the vectors connecting the the atoms in GROUPA to the atoms in GROUPB.
206 : In other words the vector between atoms 1 and 2 and the vector between atoms 1 and 3. The number of values
207 : more than 0.1 is then printed to a file.
208 :
209 : ```plumed
210 : d1: XDISTANCES GROUPA=1 GROUPB=2,3 MORE_THAN={RATIONAL R_0=0.1}
211 : PRINT ARG=d1.morethan
212 : ```
213 :
214 : */
215 : //+ENDPLUMEDOC
216 :
217 :
218 : //+PLUMEDOC MCOLVAR YDISTANCES
219 : /*
220 : Calculate the y components of the vectors connecting one or many pairs of atoms.
221 :
222 : __This shortcut action allows you to calculate functions of the distribution of the y components of the vectors connecting pairs of atoms and
223 : reproduces the syntax in older PLUMED versions.
224 : If you look at the example inputs below you can
225 : see how the new syntax operates. We would strongly encourage you to use the newer syntax as it offers greater flexibility.__
226 :
227 : The following input tells plumed to calculate the y-component of the vector connecting atom 3 to atom 5 and
228 : the y-component of the vector connecting atom 1 to atom 2. The minimum of these two quantities is then
229 : printed
230 :
231 : ```plumed
232 : d1: YDISTANCES ATOMS1=3,5 ATOMS2=1,2 MIN={BETA=0.1}
233 : PRINT ARG=d1.min
234 : ```
235 :
236 : The following input tells plumed to calculate the y-component of the vector connecting atom 3 to atom 5 and
237 : the y-component of the vector connecting atom 1 to atom 2. The number of values that are
238 : less than 0.1nm is then printed to a file.
239 :
240 : ```plumed
241 : d1: YDISTANCES ATOMS1=3,5 ATOMS2=1,2 LESS_THAN={RATIONAL R_0=0.1}
242 : PRINT ARG=d1.lessthan
243 : ```
244 :
245 : The following input tells plumed to calculate the y-components of all the distinct vectors that can be created
246 : between atoms 1, 2 and 3 (i.e. the vectors between atoms 1 and 2, atoms 1 and 3 and atoms 2 and 3).
247 : The average of these quantities is then calculated.
248 :
249 : ```plumed
250 : d1: YDISTANCES GROUP=1-3 MEAN
251 : PRINT ARG=d1.mean
252 : ```
253 :
254 : The following input tells plumed to calculate all the vectors connecting the the atoms in GROUPA to the atoms in GROUPB.
255 : In other words the vector between atoms 1 and 2 and the vector between atoms 1 and 3. The number of values
256 : more than 0.1 is then printed to a file.
257 :
258 : ```plumed
259 : d1: YDISTANCES GROUPA=1 GROUPB=2,3 MORE_THAN={RATIONAL R_0=0.1}
260 : PRINT ARG=d1.morethan
261 : ```
262 :
263 : */
264 : //+ENDPLUMEDOC
265 :
266 : //+PLUMEDOC MCOLVAR ZDISTANCES
267 : /*
268 : Calculate the z components of the vectors connecting one or many pairs of atoms.
269 :
270 : __This shortcut action allows you to calculate functions of the distribution of the z components of the vectors connecting pairs of atoms and
271 : reproduces the syntax in older PLUMED versions.
272 : If you look at the example inputs below you can
273 : see how the new syntax operates. We would strongly encourage you to use the newer syntax as it offers greater flexibility.__
274 :
275 : The following input tells plumed to calculate the z-component of the vector connecting atom 3 to atom 5 and
276 : the z-component of the vector connecting atom 1 to atom 2. The minimum of these two quantities is then
277 : printed
278 :
279 : ```plumed
280 : d1: ZDISTANCES ATOMS1=3,5 ATOMS2=1,2 MIN={BETA=0.1}
281 : PRINT ARG=d1.min
282 : ```
283 :
284 : The following input tells plumed to calculate the z-component of the vector connecting atom 3 to atom 5 and
285 : the z-component of the vector connecting atom 1 to atom 2. The number of values that are
286 : less than 0.1nm is then printed to a file.
287 :
288 : ```plumed
289 : d1: ZDISTANCES ATOMS1=3,5 ATOMS2=1,2 LESS_THAN={RATIONAL R_0=0.1}
290 : PRINT ARG=d1.lessthan
291 : ```
292 :
293 : The following input tells plumed to calculate the z-components of all the distinct vectors that can be created
294 : between atoms 1, 2 and 3 (i.e. the vectors between atoms 1 and 2, atoms 1 and 3 and atoms 2 and 3).
295 : The average of these quantities is then calculated.
296 :
297 : ```plumed
298 : d1: ZDISTANCES GROUP=1-3 MEAN
299 : PRINT ARG=d1.mean
300 : ```
301 :
302 : The following input tells plumed to calculate all the vectors connecting the the atoms in GROUPA to the atoms in GROUPB.
303 : In other words the vector between atoms 1 and 2 and the vector between atoms 1 and 3. The number of values
304 : more than 0.1 is then printed to a file.
305 :
306 : ```plumed
307 : d1: ZDISTANCES GROUPA=1 GROUPB=2,3 MORE_THAN={RATIONAL R_0=0.1}
308 : PRINT ARG=d1.morethan
309 : ```
310 :
311 : */
312 : //+ENDPLUMEDOC
313 :
314 :
315 : namespace PLMD {
316 : namespace multicolvar {
317 :
318 : class Distances : public ActionShortcut {
319 : public:
320 : static void registerKeywords(Keywords& keys);
321 : explicit Distances(const ActionOptions&);
322 : };
323 :
324 : PLUMED_REGISTER_ACTION(Distances,"DISTANCES")
325 : PLUMED_REGISTER_ACTION(Distances,"XDISTANCES")
326 : PLUMED_REGISTER_ACTION(Distances,"YDISTANCES")
327 : PLUMED_REGISTER_ACTION(Distances,"ZDISTANCES")
328 :
329 115 : void Distances::registerKeywords(Keywords& keys) {
330 115 : ActionShortcut::registerKeywords( keys );
331 230 : keys.setDeprecated("DISTANCE");
332 115 : keys.add("atoms-1","GROUP","Calculate the distance between each distinct pair of atoms in the group");
333 115 : keys.add("atoms-2","GROUPA","Calculate the distances between all the atoms in GROUPA and all "
334 : "the atoms in GROUPB. This must be used in conjunction with GROUPB.");
335 115 : keys.add("atoms-2","GROUPB","Calculate the distances between all the atoms in GROUPA and all the atoms "
336 : "in GROUPB. This must be used in conjunction with GROUPA.");
337 115 : keys.add("numbered","ATOMS","the pairs of atoms that you would like to calculate the angles for");
338 115 : keys.addFlag("NOPBC",false,"ignore the periodic boundary conditions when calculating distances");
339 115 : keys.addFlag("COMPONENTS",false,"calculate the x, y and z components of the distance separately and store them as label.x, label.y and label.z");
340 115 : keys.addFlag("SCALED_COMPONENTS",false,"calculate the a, b and c scaled components of the distance separately and store them as label.a, label.b and label.c");
341 230 : keys.addDeprecatedFlag("LOWMEM","");
342 230 : keys.reset_style("ATOMS","atoms");
343 115 : MultiColvarShortcuts::shortcutKeywords( keys );
344 115 : keys.add("atoms","ORIGIN","calculate the distance of all the atoms specified using the ATOMS keyword from this point");
345 115 : keys.add("numbered","LOCATION","the location at which the CV is assumed to be in space");
346 230 : keys.reset_style("LOCATION","atoms");
347 230 : keys.setValueDescription("vector","the DISTANCES between the each pair of atoms that were specified");
348 230 : keys.addOutputComponent("x","COMPONENTS","vector","the x-components of the distance vectors");
349 230 : keys.addOutputComponent("y","COMPONENTS","vector","the y-components of the distance vectors");
350 230 : keys.addOutputComponent("z","COMPONENTS","vector","the z-components of the distance vectors");
351 115 : keys.needsAction("GROUP");
352 115 : keys.needsAction("DISTANCE");
353 115 : keys.needsAction("CENTER");
354 115 : }
355 :
356 58 : Distances::Distances(const ActionOptions& ao):
357 : Action(ao),
358 58 : ActionShortcut(ao) {
359 : // Create distances
360 : bool lowmem;
361 58 : parseFlag("LOWMEM",lowmem);
362 58 : if( lowmem ) {
363 0 : warning("LOWMEM flag is deprecated and is no longer required for this action");
364 : }
365 58 : std::string dline = getShortcutLabel() + ": DISTANCE";
366 : bool nopbc;
367 58 : parseFlag("NOPBC",nopbc);
368 58 : if( nopbc ) {
369 : dline += " NOPBC";
370 : }
371 58 : if( getName()=="DISTANCES" ) {
372 : bool comp;
373 55 : parseFlag("COMPONENTS",comp);
374 55 : if( comp ) {
375 : dline += " COMPONENTS";
376 : }
377 : bool scomp;
378 55 : parseFlag("SCALED_COMPONENTS",scomp);
379 55 : if( scomp ) {
380 : dline += " SCALED_COMPONENTS";
381 : }
382 : } else {
383 : dline += " COMPONENTS";
384 : }
385 : // Parse origin
386 : std::string ostr;
387 116 : parse("ORIGIN",ostr);
388 : std::string num;
389 58 : if( ostr.length()>0 ) {
390 : // Parse atoms
391 : std::vector<std::string> afstr;
392 23 : MultiColvarShortcuts::parseAtomList("ATOMS",afstr,this);
393 17118 : for(unsigned i=0; i<afstr.size(); ++i) {
394 17095 : Tools::convert( i+1, num );
395 34190 : dline += " ATOMS" + num + "=" + ostr + "," + afstr[i];
396 : }
397 23 : } else {
398 : std::vector<std::string> grp;
399 70 : MultiColvarShortcuts::parseAtomList("GROUP",grp,this);
400 : std::vector<std::string> grpa;
401 70 : MultiColvarShortcuts::parseAtomList("GROUPA",grpa,this);
402 35 : if( grp.size()>0 ) {
403 3 : if( grpa.size()>0 ) {
404 0 : error("should not be using GROUPA in tandem with GROUP");
405 : }
406 : unsigned n=0;
407 203 : for(unsigned i=1; i<grp.size(); ++i) {
408 10103 : for(unsigned j=0; j<i; ++j) {
409 9903 : Tools::convert( n+1, num );
410 : n++;
411 19806 : dline += " ATOMS" + num + "=" + grp[i] + "," + grp[j];
412 : }
413 : }
414 32 : } else if( grpa.size()>0 ) {
415 : std::vector<std::string> grpb;
416 2 : MultiColvarShortcuts::parseAtomList("GROUPB",grpb,this);
417 1 : if( grpb.size()==0 ) {
418 0 : error("found GROUPA but no corresponding GROUPB");
419 : }
420 4 : for(unsigned i=0; i<grpa.size(); ++i) {
421 294 : for(unsigned j=0; j<grpb.size(); ++j) {
422 291 : Tools::convert( i*grpb.size() + j + 1, num );
423 582 : dline += " ATOMS" + num + "=" + grpa[i] + "," + grpb[j];
424 : }
425 : }
426 1 : } else {
427 31 : std::string grpstr = getShortcutLabel() + "_grp: GROUP ATOMS=";
428 31 : for(unsigned i=1;; ++i) {
429 : std::string atstring;
430 1844 : parseNumbered("ATOMS",i,atstring);
431 922 : if( atstring.length()==0 ) {
432 : break;
433 : }
434 : std::string locstr;
435 1782 : parseNumbered("LOCATION",i,locstr);
436 891 : if( locstr.length()==0 ) {
437 281 : Tools::convert( i, num );
438 562 : readInputLine( getShortcutLabel() + "_vatom" + num + ": CENTER ATOMS=" + atstring );
439 281 : if( i==1 ) {
440 54 : grpstr += getShortcutLabel() + "_vatom" + num;
441 : } else {
442 508 : grpstr += "," + getShortcutLabel() + "_vatom" + num;
443 : }
444 : } else {
445 610 : if( i==1 ) {
446 : grpstr += locstr;
447 : } else {
448 1212 : grpstr += "," + locstr;
449 : }
450 : }
451 891 : Tools::convert( i, num );
452 1782 : dline += " ATOMS" + num + "=" + atstring;
453 891 : }
454 31 : readInputLine( grpstr );
455 : }
456 35 : }
457 58 : readInputLine( dline );
458 : // Add shortcuts to label
459 58 : if( getName()=="DISTANCES" ) {
460 110 : MultiColvarShortcuts::expandFunctions( getShortcutLabel(), getShortcutLabel(), "", this );
461 3 : } else if( getName()=="XDISTANCES" ) {
462 2 : MultiColvarShortcuts::expandFunctions( getShortcutLabel(), getShortcutLabel() + ".x", "", this );
463 2 : } else if( getName()=="YDISTANCES" ) {
464 2 : MultiColvarShortcuts::expandFunctions( getShortcutLabel(), getShortcutLabel() + ".y", "", this );
465 1 : } else if( getName()=="ZDISTANCES" ) {
466 2 : MultiColvarShortcuts::expandFunctions( getShortcutLabel(), getShortcutLabel() + ".z", "", this );
467 : }
468 58 : }
469 :
470 : }
471 : }
|
