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 unecessary 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 178 : void CoordinationNumbers::shortcutKeywords( Keywords& keys ) {
251 178 : ActionShortcut::registerKeywords( keys );
252 178 : 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 178 : 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 178 : 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 178 : keys.add("compulsory","NN","6","The n parameter of the switching function ");
256 178 : keys.add("compulsory","MM","0","The m parameter of the switching function; 0 implies 2*NN");
257 178 : keys.add("compulsory","D_0","0.0","The d_0 parameter of the switching function");
258 178 : keys.add("compulsory","R_0","The r_0 parameter of the switching function");
259 178 : keys.add("optional","SWITCH","the switching function that it used in the construction of the contact matrix");
260 178 : keys.add("optional","MASK","the label for a vector that is used to determine which rows of the matrix are computed");
261 356 : keys.linkActionInDocs("SWITCH","LESS_THAN");
262 178 : multicolvar::MultiColvarShortcuts::shortcutKeywords( keys );
263 178 : keys.needsAction("CONTACT_MATRIX");
264 178 : keys.needsAction("GROUP");
265 178 : keys.addDOI("10.1021/ct500320c");
266 178 : }
267 :
268 111 : 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 111 : if( sp_str.length()==0 && spa_str.length()==0 ) {
271 0 : return;
272 : }
273 :
274 111 : std::string matinp = lab + "_mat: CONTACT_MATRIX";
275 111 : if( sp_str.length()>0 ) {
276 72 : matinp += " GROUP=" + sp_str;
277 144 : 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 222 : action->parse("SWITCH",sw_str);
285 111 : if( sw_str.length()>0 ) {
286 168 : matinp += " SWITCH={" + sw_str + "}";
287 : } else {
288 : std::string r0;
289 54 : action->parse("R_0",r0);
290 : std::string d0;
291 54 : action->parse("D_0",d0);
292 27 : if( r0.length()==0 ) {
293 0 : action->error("missing switching function parameters use SWITCH/R_0");
294 : }
295 : std::string nn;
296 54 : action->parse("NN",nn);
297 : std::string mm;
298 27 : action->parse("MM",mm);
299 54 : matinp += " R_0=" + r0 + " D_0=" + d0 + " NN=" + nn + " MM=" + mm;
300 : }
301 111 : if( components ) {
302 : matinp += " COMPONENTS";
303 : }
304 : std::string maskstr;
305 222 : action->parse("MASK",maskstr);
306 111 : if( maskstr.length()>0 ) {
307 8 : matinp += " MASK=" + maskstr;
308 : }
309 111 : action->readInputLine( matinp );
310 : }
311 :
312 68 : void CoordinationNumbers::registerKeywords( Keywords& keys ) {
313 68 : shortcutKeywords( keys );
314 136 : if( keys.getDisplayName()=="COORDINATION_MOMENTS" ) {
315 7 : keys.add("compulsory","R_POWER","the power to which you want to raise the distance");
316 : }
317 136 : keys.addDeprecatedFlag("LOWMEM","");
318 68 : keys.add("optional","MOMENTS","the list of moments that you would like to calculate");
319 136 : keys.addOutputComponent("moment","MOMENTS","scalar","the moments of the distribution");
320 136 : keys.reset_style("MOMENTS","deprecated");
321 68 : keys.needsAction("MATRIX_VECTOR_PRODUCT");
322 68 : keys.needsAction("ONES");
323 68 : keys.needsAction("MOMENTS");
324 136 : keys.setValueDescription("vector","the coordination numbers of the specified atoms");
325 68 : }
326 :
327 47 : CoordinationNumbers::CoordinationNumbers(const ActionOptions& ao):
328 : Action(ao),
329 47 : ActionShortcut(ao) {
330 : bool lowmem;
331 47 : parseFlag("LOWMEM",lowmem);
332 47 : 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 47 : parse("SPECIES",sp_str);
338 47 : parse("SPECIESA",specA);
339 94 : parse("SPECIESB",specB);
340 47 : if( sp_str.length()>0 || specA.length()>0 ) {
341 47 : matlab = getShortcutLabel() + "_mat";
342 47 : bool comp=false;
343 47 : if( getName()=="COORDINATION_MOMENTS" ) {
344 3 : comp=true;
345 6 : matlab = getShortcutLabel() + "_mat";
346 : }
347 47 : expandMatrix( comp, getShortcutLabel(), sp_str, specA, specB, this );
348 : } else {
349 0 : error("missing atoms input use SPECIES or SPECIESA/SPECIESB");
350 : }
351 47 : ActionWithValue* mb=plumed.getActionSet().selectWithLabel<ActionWithValue*>( matlab );
352 47 : if( !mb ) {
353 0 : error("could not find action with name " + matlab );
354 : }
355 47 : Value* arg=mb->copyOutput(0);
356 47 : 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 47 : Tools::convert( arg->getShape()[1], nones );
362 94 : readInputLine( getShortcutLabel() + "_ones: ONES SIZE=" + nones );
363 47 : if( getName()=="COORDINATION_MOMENTS" ) {
364 : // Calculate the lengths of the vectors
365 : std::string r_power;
366 3 : parse("R_POWER",r_power);
367 9 : readInputLine( getShortcutLabel() + "_pow: CUSTOM ARG=" + matlab + ".x," + matlab + ".y," + matlab + ".z," + matlab + ".w VAR=x,y,z,w "
368 6 : + "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 94 : readInputLine( getShortcutLabel() + ": MATRIX_VECTOR_PRODUCT ARG=" + matlab + "," + getShortcutLabel() + "_ones");
373 : std::vector<std::string> moments;
374 47 : parseVector("MOMENTS",moments);
375 47 : Tools::interpretRanges( moments );
376 47 : 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 4 : readInputLine( getShortcutLabel() + "_diffpow-" + moments[i] + ": CUSTOM ARG=" + getShortcutLabel() + "," + getShortcutLabel() + "_caverage PERIODIC=NO FUNC=(x-y)^" + moments[i] );
380 4 : 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 47 : multicolvar::MultiColvarShortcuts::readShortcutKeywords( keymap, this );
386 94 : multicolvar::MultiColvarShortcuts::expandFunctions( getShortcutLabel(), getShortcutLabel(), "", keymap, this );
387 94 : }
388 :
389 :
390 : }
391 : }
|
