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 "CoordinationNumbers.h"
23 : #include "multicolvar/MultiColvarShortcuts.h"
24 : #include "core/ActionWithValue.h"
25 : #include "core/ActionRegister.h"
26 : #include "core/PlumedMain.h"
27 : #include "core/ActionSet.h"
28 :
29 : #include <string>
30 : #include <cmath>
31 :
32 : namespace PLMD {
33 : namespace symfunc {
34 :
35 : //+PLUMEDOC MCOLVAR COORDINATIONNUMBER
36 : /*
37 : Calculate the coordination numbers of atoms so that you can then calculate functions of the distribution of
38 : coordination numbers such as the minimum, the number less than a certain quantity and so on.
39 :
40 : The coordination number of a atom $i$ is the number of atoms that are within a certain cutoff distance of it.
41 : This quantity can be calculated as follows:
42 :
43 : $$
44 : s_i = \sum_j \sigma(r_{ij})
45 : $$
46 :
47 : where $r_{ij}$ is the distance between atoms $i$ and $j$ and $\sigma$ is a switching function. The typical switching
48 : function that is used in metadynamics is this one:
49 :
50 : $$
51 : s(r) = \frac{ 1 - \left(\frac{r-d_0}{r_0}\right)^n } { 1 - \left(\frac{r-d_0}{r_0}\right)^m }
52 : $$
53 :
54 : The following example shows how you can use this shortcut action to calculate and print the coordination numbers of
55 : one hundred atoms with themselves:
56 :
57 : ```plumed
58 : c: COORDINATIONNUMBER SPECIES=1-100 R_0=1.0
59 : DUMPATOMS ATOMS=c ARG=c FILE=coords.xyz
60 : ```
61 :
62 : This input will produce an output file called coords that contains the coordination numbers of the 100 input atoms. The cutoff
63 : that is used to calculate the coordination number in this case is 1.0. In the input above we use a rational [switching function](LESS_THAN.md)
64 : with the parameters above. We would recommend using SWITCH syntax
65 : rather than the syntax above when giving the parameters for the switching function as you can then use any of the switching functions described
66 : in the documentation for [LESS_THAN](LESS_THAN.md). More importantly, however, using this syntax allows you to set the D_MAX parameter for the
67 : switching function as demonstrated below:
68 :
69 : ```plumed
70 : c: COORDINATIONNUMBER SPECIES=1-100 SWITCH={RATIONAL R_0=1.0}
71 : DUMPATOMS ATOMS=c ARG=c FILE=coords.xyz
72 : ```
73 :
74 : Setting the `D_MAX` can substantially improve PLUMED performance as it turns on the linked list algorithm that is discussed in the optimisation details part
75 : of the documentation for [CONTACT_MATRIX](CONTACT_MATRIX.md).
76 :
77 : The vectors that are output by the COORDINATIONNUMBER shortcut can be used in the input for many other functions that are within
78 : PLUMED. In addition, in order to ensure compatibility with older versions of PLUMED you can add additional keywords on the input
79 : line for COORDINATIONNUMBER in order to calculate various functions of the input. For example, the following input tells plumed
80 : ato calculate the coordination numbers of atoms 1-100 with themselves. The minimum coordination number is then calculated.
81 :
82 : ```plumed
83 : c: COORDINATIONNUMBER SPECIES=1-100 R_0=1.0 MIN={BETA=0.1}
84 : ```
85 :
86 : By constrast, this input tells plumed to calculate how many atoms from 1-100 are within 3.0 of each of the atoms
87 : from 101-110. In the first 101 is the central atom, in the second 102 is the central atom and so on. The
88 : number of coordination numbers that are more than 6 is then computed.
89 :
90 : ```plumed
91 : c: COORDINATIONNUMBER SPECIESA=101-110 SPECIESB=1-100 R_0=3.0
92 : mt: MORE_THAN ARG=c SWITCH={RATIONAL R_0=6.0 NN=6 MM=12 D_0=0}
93 : s: SUM ARG=mt PERIODIC=NO
94 : ```
95 :
96 : ## The MASK keyword
97 :
98 : Supppose that you want to calculate the average coordination number for the atoms that are within a sphere in the center of your simulation box. You can do so by exploiting an input similar to the one shown
99 : below:
100 :
101 : ```plumed
102 : # Fixed virtual atom which serves as the probe volume's center (pos. in nm)
103 : center: FIXEDATOM AT=2.5,2.5,2.5
104 : # Vector in which element i is one if atom i is in sphere of interest and zero otherwise
105 : sphere: INSPHERE ATOMS=1-400 CENTER=center RADIUS={GAUSSIAN D_0=0.5 R_0=0.01 D_MAX=0.52}
106 : # Calculate the coordination numbers
107 : cc: COORDINATIONNUMBER SPECIES=1-400 SWITCH={RATIONAL D_0=3.0 R_0=1.5 D_MAX=6.0}
108 : # Multiply fccubic parameters numbers by sphere vector
109 : prod: CUSTOM ARG=cc,sphere FUNC=x*y PERIODIC=NO
110 : # Sum of coordination numbers for atoms that are in the sphere of interest
111 : numer: SUM ARG=prod PERIODIC=NO
112 : # Number of atoms that are in sphere of interest
113 : denom: SUM ARG=sphere PERIODIC=NO
114 : # Average coordination number for atoms in sphere of interest
115 : av: CUSTOM ARG=numer,denom FUNC=x/y PERIODIC=NO
116 : # And print out final CV to a file
117 : PRINT ARG=av FILE=colvar STRIDE=1
118 : ```
119 :
120 : This calculation is slow because you have to calculate the coordination numbers of all the atoms even though only a small subset of these quanitties are required to compute the average coordination number in the
121 : sphere. To avoid all these unnecessary calculations you use the `MASK` keyword as shown below:
122 :
123 : ```plumed
124 : # Fixed virtual atom which serves as the probe volume's center (pos. in nm)
125 : center: FIXEDATOM AT=2.5,2.5,2.5
126 : # Vector in which element i is one if atom i is in sphere of interest and zero otherwise
127 : sphere: INSPHERE ATOMS=1-400 CENTER=center RADIUS={GAUSSIAN D_0=0.5 R_0=0.01 D_MAX=0.52}
128 : # Calculate the coordination numbers
129 : cc: COORDINATIONNUMBER SPECIES=1-400 MASK=sphere SWITCH={RATIONAL D_0=3.0 R_0=1.5 D_MAX=6.0}
130 : # Multiply fccubic parameters numbers by sphere vector
131 : prod: CUSTOM ARG=cc,sphere FUNC=x*y PERIODIC=NO
132 : # Sum of coordination numbers for atoms that are in the sphere of interest
133 : numer: SUM ARG=prod PERIODIC=NO
134 : # Number of atoms that are in sphere of interest
135 : denom: SUM ARG=sphere PERIODIC=NO
136 : # Average coordination number for atoms in sphere of interest
137 : av: CUSTOM ARG=numer,denom FUNC=x/y PERIODIC=NO
138 : # And print out final CV to a file
139 : PRINT ARG=av FILE=colvar STRIDE=1
140 : ```
141 :
142 : Adding the instruction `MASK=sphere` to the CONTACT_MATRIX line in this input tells PLUMED to only calculate the $i$th row in the adjacency matrix if the $i$th element of the vector `sphere` is non-zero.
143 : In other words, by adding this command we have ensured that we are not calculating coordination numbers for atoms that are not in the sphere that is of interest. In this way we can thus reduce the computational
144 : expense of the calculation enormously.
145 :
146 : Notice, that there are other places where we can use this same trick. Some further examples are given in the documentation for [CONTACT_MATRIX](CONTACT_MATRIX.md).
147 :
148 : */
149 : //+ENDPLUMEDOC
150 :
151 : //+PLUMEDOC MCOLVAR COORDINATION_MOMENTS
152 : /*
153 : Calculate moments of the distribution of distances in the first coordination sphere
154 :
155 : This is the CV that was developed by White and Voth and is described in the paper in the bibliograhy below. This action provides a way of indirectly biasing radial distribution functions and computes the following function
156 :
157 : $$
158 : s_i = \sum_j \sigma(r_{ij})r_{ij}^k
159 : $$
160 :
161 : where $k$ is the value that is input using the R_POWER keyword, $r_{ij}$ is the distance between atoms $i$ and $j$ and $\sigma$ is a switching function.
162 :
163 : The following example shows how this action can be used.
164 :
165 : ```plumed
166 : cn1: COORDINATION_MOMENTS SPECIES=1-10 R_0=1.0 R_POWER=1
167 : cn1_mean: MEAN ARG=cn1 PERIODIC=NO
168 : PRINT ARG=cn1_mean FILE=colvar
169 : ```
170 :
171 : As you can see, the action works similarlly to [COORDINATIONNUMBER](COORDINATIONNUMBER.md).
172 :
173 : In the input above we use a rational [switching function](LESS_THAN.md)
174 : with the parameters above. We would recommend using SWITCH syntax
175 : rather than the syntax above when giving the parameters for the switching function as you can then use any of the switching functions described
176 : in the documentation for [LESS_THAN](LESS_THAN.md). More importantly, however, using this syntax allows you to set the D_MAX parameter for the
177 : switching function as demonstrated below:
178 :
179 :
180 : ```plumed
181 : cn0: COORDINATIONNUMBER SPECIES=1-10 SWITCH={RATIONAL R_0=1.0 D_MAX=8}
182 : cn0_mean: MEAN ARG=cn0 PERIODIC=NO
183 : cn1: COORDINATION_MOMENTS SPECIES=1-10 SWITCH={RATIONAL R_0=1.0 D_MAX=8} R_POWER=1
184 : cn1_mean: MEAN ARG=cn1 PERIODIC=NO
185 : cn2: COORDINATION_MOMENTS SPECIES=1-10 SWITCH={RATIONAL R_0=1.0 D_MAX=8} R_POWER=2
186 : cn2_mean: MEAN ARG=cn2 PERIODIC=NO
187 : PRINT ARG=cn0_mean,cn1_mean,cn2_mean STRIDE=1 FILE=cn_out
188 : ```
189 :
190 : Setting the `D_MAX` can substantially improve PLUMED performance as it turns on the linked list algorithm that is discussed in the optimisation details part
191 : of the documentation for [CONTACT_MATRIX](CONTACT_MATRIX.md).
192 :
193 : ## Working with two types of atom
194 :
195 : If you would like a way of indirectly biasing the radial distribution function that describes how the atoms in GROUPB are arranged around the atoms in GROUPA you use an input like the one
196 : shown below:
197 :
198 : ```plumed
199 : d: COORDINATION_MOMENTS ...
200 : SPECIESA=1-64 SPECIESB=65-200 R_POWER=1
201 : SWITCH={RATIONAL D_0=3.0 R_0=1.5 D_MAX=6.0}
202 : ...
203 : s: MEAN ARG=d PERIODIC=NO
204 : PRINT ARG=s FILE=colv
205 : ```
206 :
207 : ## The MASK keyword
208 :
209 : You can use the MASK keyword with this action in the same way that it is used with [COORDINATIONNUMBER](COORDINATIONNUMBER.md). This keyword thus expects a vector in
210 : input, which tells COORDINATION_MOMENTS that it is safe not to calculate the COORDINATION_MOMENTS parameter for some of the atoms. As illustrated below, this is useful if you are using functionality
211 : from the [volumes module](module_volumes.md) to calculate the average value of the COORDINATION_MOMENTS parameter for only those atoms that lie in a certain part of the simulation box.
212 :
213 : ```plumed
214 : # Fixed virtual atom which serves as the probe volume's center (pos. in nm)
215 : center: FIXEDATOM AT=2.5,2.5,2.5
216 : # Vector in which element i is one if atom i is in sphere of interest and zero otherwise
217 : sphere: INSPHERE ATOMS=1-400 CENTER=center RADIUS={GAUSSIAN D_0=0.5 R_0=0.01 D_MAX=0.52}
218 : # Calculate the coordination moments of the atoms
219 : cc: COORDINATION_MOMENTS ...
220 : SPECIES=1-400 MASK=sphere R_POWER=1
221 : SWITCH={RATIONAL D_0=3.0 R_0=1.5 D_MAX=6.0}
222 : ...
223 : # Multiply coordination moments by sphere vector
224 : prod: CUSTOM ARG=cc,sphere FUNC=x*y PERIODIC=NO
225 : # Sum of coordination numbers for atoms that are in the sphere of interest
226 : numer: SUM ARG=prod PERIODIC=NO
227 : # Number of atoms that are in sphere of interest
228 : denom: SUM ARG=sphere PERIODIC=NO
229 : # Average coordination number for atoms in sphere of interest
230 : av: CUSTOM ARG=numer,denom FUNC=x/y PERIODIC=NO
231 : # And print out final CV to a file
232 : PRINT ARG=av FILE=colvar STRIDE=1
233 : ```
234 :
235 : This input calculate the average value of the COORDINATION_MOMENTS parameter for only those atoms that are within a spherical region that is centered on the point
236 : $(2.5,2.5,2.5)$.
237 :
238 : ## Deprecated syntax
239 :
240 : More information on the deprecated keywords that are given below is available in the documentation for the [DISTANCES](DISTANCES.md) command.
241 :
242 :
243 : */
244 : //+ENDPLUMEDOC
245 :
246 :
247 : PLUMED_REGISTER_ACTION(CoordinationNumbers,"COORDINATIONNUMBER")
248 : PLUMED_REGISTER_ACTION(CoordinationNumbers,"COORDINATION_MOMENTS")
249 :
250 214 : void CoordinationNumbers::shortcutKeywords( Keywords& keys ) {
251 214 : ActionShortcut::registerKeywords( keys );
252 214 : keys.add("atoms-3","SPECIES","the list of atoms for which the symmetry function is being calculated and the atoms that can be in the environments");
253 214 : keys.add("atoms-4","SPECIESA","the list of atoms for which the symmetry function is being calculated. This keyword must be used in conjunction with SPECIESB, which specifies the atoms that are in the environment.");
254 214 : keys.add("atoms-4","SPECIESB","the list of atoms that can be in the environments of each of the atoms for which the symmetry function is being calculated. This keyword must be used in conjunction with SPECIESA, which specifies the atoms for which the symmetry function is being calculated.");
255 214 : keys.add("compulsory","NN","6","The n parameter of the switching function ");
256 214 : keys.add("compulsory","MM","0","The m parameter of the switching function; 0 implies 2*NN");
257 214 : keys.add("compulsory","D_0","0.0","The d_0 parameter of the switching function");
258 214 : keys.add("compulsory","R_0","The r_0 parameter of the switching function");
259 214 : keys.add("optional","SWITCH","the switching function that it used in the construction of the contact matrix");
260 214 : keys.add("optional","MASK","the label for a vector that is used to determine which rows of the matrix are computed");
261 428 : keys.linkActionInDocs("SWITCH","LESS_THAN");
262 214 : multicolvar::MultiColvarShortcuts::shortcutKeywords( keys );
263 214 : keys.needsAction("CONTACT_MATRIX");
264 214 : keys.needsAction("GROUP");
265 214 : keys.addDOI("10.1021/ct500320c");
266 214 : }
267 :
268 114 : void CoordinationNumbers::expandMatrix( const bool& components, const std::string& lab, const std::string& sp_str,
269 : const std::string& spa_str, const std::string& spb_str, ActionShortcut* action ) {
270 114 : if( sp_str.length()==0 && spa_str.length()==0 ) {
271 0 : return;
272 : }
273 :
274 114 : std::string matinp = lab + "_mat: CONTACT_MATRIX";
275 114 : if( sp_str.length()>0 ) {
276 75 : matinp += " GROUP=" + sp_str;
277 150 : action->readInputLine( lab + "_grp: GROUP ATOMS=" + sp_str );
278 39 : } else if( spa_str.length()>0 ) {
279 78 : matinp += " GROUPA=" + spa_str + " GROUPB=" + spb_str;
280 78 : action->readInputLine( lab + "_grp: GROUP ATOMS=" + spa_str );
281 : }
282 :
283 : std::string sw_str;
284 228 : action->parse("SWITCH",sw_str);
285 114 : if( sw_str.length()>0 ) {
286 86 : matinp += " SWITCH={" + sw_str + "}";
287 : } else {
288 : std::string r0;
289 56 : action->parse("R_0",r0);
290 : std::string d0;
291 56 : action->parse("D_0",d0);
292 28 : if( r0.length()==0 ) {
293 0 : action->error("missing switching function parameters use SWITCH/R_0");
294 : }
295 : std::string nn;
296 56 : action->parse("NN",nn);
297 : std::string mm;
298 28 : action->parse("MM",mm);
299 112 : matinp += " R_0=" + r0 + " D_0=" + d0 + " NN=" + nn + " MM=" + mm;
300 : }
301 114 : if( components ) {
302 : matinp += " COMPONENTS";
303 : }
304 : std::string maskstr;
305 228 : action->parse("MASK",maskstr);
306 114 : if( maskstr.length()>0 ) {
307 4 : matinp += " MASK=" + maskstr;
308 : }
309 114 : action->readInputLine( matinp );
310 118 : }
311 :
312 74 : void CoordinationNumbers::registerKeywords( Keywords& keys ) {
313 74 : shortcutKeywords( keys );
314 74 : if( keys.getDisplayName()=="COORDINATION_MOMENTS" ) {
315 9 : keys.add("compulsory","R_POWER","the power to which you want to raise the distance");
316 : }
317 148 : keys.addDeprecatedFlag("LOWMEM","");
318 74 : keys.add("optional","MOMENTS","the list of moments that you would like to calculate");
319 148 : keys.addOutputComponent("moment","MOMENTS","scalar","the moments of the distribution");
320 148 : keys.reset_style("MOMENTS","deprecated");
321 74 : keys.needsAction("MATRIX_VECTOR_PRODUCT");
322 74 : keys.needsAction("ONES");
323 74 : keys.needsAction("MOMENTS");
324 148 : keys.setValueDescription("vector","the coordination numbers of the specified atoms");
325 74 : }
326 :
327 48 : CoordinationNumbers::CoordinationNumbers(const ActionOptions& ao):
328 : Action(ao),
329 48 : ActionShortcut(ao) {
330 : bool lowmem;
331 48 : parseFlag("LOWMEM",lowmem);
332 48 : if( lowmem ) {
333 0 : warning("LOWMEM flag is deprecated and is no longer required for this action");
334 : }
335 : // Setup the contract matrix if that is what is needed
336 : std::string matlab, sp_str, specA, specB;
337 48 : parse("SPECIES",sp_str);
338 48 : parse("SPECIESA",specA);
339 96 : parse("SPECIESB",specB);
340 48 : if( sp_str.length()>0 || specA.length()>0 ) {
341 48 : matlab = getShortcutLabel() + "_mat";
342 48 : bool comp=false;
343 48 : if( getName()=="COORDINATION_MOMENTS" ) {
344 3 : comp=true;
345 6 : matlab = getShortcutLabel() + "_mat";
346 : }
347 48 : expandMatrix( comp, getShortcutLabel(), sp_str, specA, specB, this );
348 : } else {
349 0 : error("missing atoms input use SPECIES or SPECIESA/SPECIESB");
350 : }
351 48 : ActionWithValue* mb=plumed.getActionSet().selectWithLabel<ActionWithValue*>( matlab );
352 48 : if( !mb ) {
353 0 : error("could not find action with name " + matlab );
354 : }
355 48 : Value* arg=mb->copyOutput(0);
356 48 : if( arg->getRank()!=2 || arg->hasDerivatives() ) {
357 0 : error("the input to this action should be a matrix or scalar");
358 : }
359 : // Create vector of ones to multiply input matrix by
360 : std::string nones;
361 48 : Tools::convert( arg->getShape()[1], nones );
362 48 : readInputLine( getShortcutLabel() + "_ones: ONES SIZE=" + nones );
363 48 : if( getName()=="COORDINATION_MOMENTS" ) {
364 : // Calculate the lengths of the vectors
365 : std::string r_power;
366 3 : parse("R_POWER",r_power);
367 15 : readInputLine( getShortcutLabel() + "_pow: CUSTOM ARG=" + matlab + ".x," + matlab + ".y," + matlab + ".z," + matlab + ".w VAR=x,y,z,w "
368 3 : + "PERIODIC=NO FUNC=w*(sqrt(x*x+y*y+z*z)^" + r_power +")");
369 6 : matlab = getShortcutLabel() + "_pow";
370 : }
371 : // Calcualte coordination numbers as matrix vector times vector of ones
372 144 : readInputLine( getShortcutLabel() + ": MATRIX_VECTOR_PRODUCT ARG=" + matlab + "," + getShortcutLabel() + "_ones");
373 48 : std::vector<std::string> moments;
374 48 : parseVector("MOMENTS",moments);
375 48 : Tools::interpretRanges( moments );
376 48 : if( moments.size()>0 ) {
377 2 : readInputLine( getShortcutLabel() + "_caverage: MEAN ARG=" + getShortcutLabel() + " PERIODIC=NO");
378 3 : for(unsigned i=0; i<moments.size(); ++i) {
379 10 : readInputLine( getShortcutLabel() + "_diffpow-" + moments[i] + ": CUSTOM ARG=" + getShortcutLabel() + "," + getShortcutLabel() + "_caverage PERIODIC=NO FUNC=(x-y)^" + moments[i] );
380 10 : readInputLine( getShortcutLabel() + "_moment-" + moments[i] + ": MEAN ARG=" + getShortcutLabel() + "_diffpow-" + moments[i] + " PERIODIC=NO");
381 : }
382 : }
383 : // Read in all the shortcut stuff
384 : std::map<std::string,std::string> keymap;
385 48 : multicolvar::MultiColvarShortcuts::readShortcutKeywords( keymap, this );
386 96 : multicolvar::MultiColvarShortcuts::expandFunctions( getShortcutLabel(), getShortcutLabel(), "", keymap, this );
387 96 : }
388 :
389 :
390 : }
391 : }
|
