Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2012-2017 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 "core/ActionWithValue.h"
25 : #include "multicolvar/MultiColvarShortcuts.h"
26 : #include "core/PlumedMain.h"
27 : #include "core/ActionSet.h"
28 : #include "CoordinationNumbers.h"
29 : #include <string>
30 : #include <cmath>
31 :
32 : namespace PLMD {
33 : namespace symfunc {
34 :
35 : //+PLUMEDOC MCOLVARF LOCAL_AVERAGE
36 : /*
37 : Calculate averages over spherical regions centered on atoms
38 :
39 : As is explained in <a href="http://www.youtube.com/watch?v=iDvZmbWE5ps"> this video </a> certain PLUMED actions
40 : calculate one scalar quantity or one vector for each of the atoms in the system. For example
41 : [COORDINATIONNUMBER](COORDINATIONNUMBER.md) measures the coordination number of each of the atoms in the system and [Q4](Q4.md) measures
42 : the 4th order Steinhardt parameter for each of the atoms in the system. These quantities provide tell us something about
43 : the disposition of the atoms in the first coordination sphere of each of the atoms of interest. In the paper in the bibliography Lechner and Dellago
44 : have suggested that one can probe local order in a system by taking the average value of such symmetry functions over
45 : the atoms within a spherical cutoff of each of these atoms in the systems. When this is done with Steinhardt parameters
46 : they claim this gives a coordinate that is better able to distinguish solid and liquid configurations of Lennard-Jones atoms.
47 :
48 : You can calculate such locally averaged quantities within plumed by using the LOCAL_AVERAGE command. This command calculates
49 : the following atom-centered quantities:
50 :
51 : $$
52 : s_i = \frac{ c_i + \sum_j \sigma(r_{ij})c_j }{ 1 + \sum_j \sigma(r_{ij}) }
53 : $$
54 :
55 : where the $c_i$ and $c_j$ values can be any vector of [symmetry functions](https://www.plumed-tutorials.org/lessons/23/001/data/SymmetryFunction.html)
56 : that can be calculated using plumed multicolvars. The function $\sigma( r_{ij} )$ is a switching function that acts on the distance between
57 : atoms $i$ and $j$. Lechner and Dellago suggest that the parameters of this function should be set so that it the function is equal to one
58 : when atom $j$ is in the first coordination sphere of atom $i$ and is zero otherwise.
59 :
60 : To see how this works in practice consider the following example input.
61 :
62 : ```plumed
63 : d1: COORDINATIONNUMBER SPECIES=1-64 SWITCH={RATIONAL D_0=1.3 R_0=0.2 D_MAX=3.0}
64 : la: LOCAL_AVERAGE SPECIES=d1 D_0=1.3 R_0=0.2 NN=6 MM=12
65 : mtf: MORE_THAN ARG=la SWITCH={RATIONAL R_0=4}
66 : la_morethan: SUM ARG=mtf PERIODIC=NO
67 : PRINT ARG=la_morethan FILE=colvar
68 : ```
69 :
70 : This input calculates the coordination numbers for all the atoms in the system. These coordination numbers are then averaged over
71 : spherical regions. The number of averaged coordination numbers that are greater than 4 is then output to a file. Furthermore, if you
72 : expand the input above you can see how the LOCAL_AVERAGE command is a shortcut action that expands to a longer input that you should be able to
73 : interpret.
74 :
75 : In the input above we use a rational [switching function](LESS_THAN.md) with the parameters above. We would recommend using SWITCH syntax
76 : rather than the syntax above when giving the parameters for the switching function as you can then use any of the switching functions described
77 : 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
78 : switching function as demonstrated below:
79 :
80 : ```plumed
81 : d1: COORDINATIONNUMBER SPECIES=1-64 SWITCH={RATIONAL D_0=1.3 R_0=0.2 D_MAX=3.0}
82 : la: LOCAL_AVERAGE SPECIES=d1 SWITCH={RATIONAL D_0=1.3 R_0=0.2 D_MAX=3.0}
83 : mtf: MORE_THAN ARG=la SWITCH={RATIONAL R_0=4}
84 : la_morethan: SUM ARG=mtf PERIODIC=NO
85 : PRINT ARG=la_morethan FILE=colvar
86 : ```
87 :
88 : 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
89 : of the documentation for [CONTACT_MATRIX](CONTACT_MATRIX.md).
90 :
91 : ## Locally averaged Steinhardt Parameters
92 :
93 : What Lechner and Dellago did in their paper was a little more complicated than the example in the previous section. To reproduce what they did you would use
94 : an input something like this:
95 :
96 : ```plumed
97 : q4: Q4 SPECIES=1-64 SWITCH={RATIONAL D_0=1.3 R_0=0.2}
98 : la: LOCAL_AVERAGE SPECIES=q4 SWITCH={RATIONAL D_0=1.3 R_0=0.2}
99 : la_mean: MEAN ARG=la PERIODIC=NO
100 : PRINT ARG=la_mean FILE=colvar
101 : ```
102 :
103 : This example input calculates the [Q4](Q4.md) vectors for each of the atoms in the system. These vectors are then averaged
104 : component by component over a spherical region. The average value for this quantity is then outputeed to a file. If you want
105 : to understand more about the shortcut that is used here you can read [this page](https://www.plumed-tutorials.org/lessons/23/001/data/Steinhardt.html).
106 :
107 : ## Using SPECIESA/SPECIESB
108 :
109 : If you only want to calculate the local averages for a subset of the atoms you might use an input like the one shown below:
110 :
111 : ```plumed
112 : # Calculate the coordination numbers for the two types of atom
113 : d1a: COORDINATIONNUMBER SPECIES=1-5 SWITCH={RATIONAL D_0=1.3 R_0=0.2 D_MAX=3.0}
114 : d1b: COORDINATIONNUMBER SPECIES=1-64 SWITCH={RATIONAL D_0=1.3 R_0=0.2 D_MAX=3.0}
115 : # Now compute the local average
116 : la: LOCAL_AVERAGE SPECIESA=d1a SPECIESB=d1b SWITCH={RATIONAL D_0=1.3 R_0=0.2 D_MAX=3.0}
117 : DUMPATOMS ARG=la ATOMS=1-5 FILE=atoms.xyz
118 : ```
119 :
120 : In the input above the local averages are computed for atoms 1-5 only. However, to calculate these local averages we still need to calculate the coordination numbers of all
121 : the atoms in the system.
122 :
123 : ## The MASK keyword
124 :
125 : 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
126 : input, which tells PLUMED the atoms for which you do not need to calculate the function. As illustrated below, this is useful if you are using functionality
127 : from the [volumes module](module_volumes.md) to calculate the average value of the local averages for only those atoms that
128 : lie in a certain part of the simulation box.
129 :
130 : ```plumed
131 : # Calculate the coordination number for all the atoms
132 : cc: COORDINATIONNUMBER SPECIES=1-400 SWITCH={RATIONAL D_0=1.3 R_0=0.2 D_MAX=3.0}
133 : # Fixed virtual atom which serves as the probe volume's center (pos. in nm)
134 : center: FIXEDATOM AT=2.5,2.5,2.5
135 : # Vector in which element i is one if atom i is in sphere of interest and zero otherwise
136 : sphere: INSPHERE ATOMS=1-400 CENTER=center RADIUS={GAUSSIAN D_0=0.5 R_0=0.01 D_MAX=0.52}
137 : # Calculate the local average
138 : lcc: LOCAL_AVERAGE SPECIES=cc SWITCH={RATIONAL D_0=1.3 R_0=0.2 D_MAX=3.0} MASK=sphere
139 : # Multiply local averages parameters numbers by sphere vector
140 : prod: CUSTOM ARG=lcc,sphere FUNC=x*y PERIODIC=NO
141 : # Sum of coordination numbers for atoms that are in the sphere of interest
142 : numer: SUM ARG=prod PERIODIC=NO
143 : # Number of atoms that are in sphere of interest
144 : denom: SUM ARG=sphere PERIODIC=NO
145 : # Average coordination number for atoms in sphere of interest
146 : av: CUSTOM ARG=numer,denom FUNC=x/y PERIODIC=NO
147 : # And print out final CV to a file
148 : PRINT ARG=av FILE=colvar STRIDE=1
149 : ```
150 :
151 : This input calculates the average value of the local average for those atoms that are within a spherical region that is centered on the point $(2.5,2.5,2.5)$.
152 : By including the MASK keyword in the LOCAL_AVERAGE line we reduce the number of local averages that are computed. These quantities are only calculated for those
153 : atoms that are within the spherical region of interest. However, we are still asking PLUMED to calculate the coordination numbers for many atoms that
154 : will not contribute to the final averaged quantity in the COORDINATIONNUMBER commmand with label `cc`. If you use an input like the one shown below you can ensure
155 : that only those coordination numbers that compute to the final average are computed.
156 :
157 : ```plumed
158 : # Fixed virtual atom which serves as the probe volume's center (pos. in nm)
159 : center: FIXEDATOM AT=2.5,2.5,2.5
160 : # Vector in which element i is one if atom i is in sphere of interest and zero otherwise
161 : sphere: INSPHERE ATOMS=1-400 CENTER=center RADIUS={GAUSSIAN D_0=0.5 R_0=0.01 D_MAX=0.52}
162 : ones: ONES SIZE=400
163 : # Calculate a elements of adjacency matrix for those elements that are inside the sphere
164 : vol_cmap: CONTACT_MATRIX GROUP=1-400 SWITCH={RATIONAL D_0=1.3 R_0=0.2 D_MAX=3.0} MASK=sphere
165 : # Create a matrix that for the sphere stuff
166 : sphere_mat: OUTER_PRODUCT ARG=sphere,ones MASK=vol_cmap
167 : # Multiply this contact matrix by the map
168 : bonds_mat: CUSTOM ARG=vol_cmap,sphere_mat FUNC=x*y PERIODIC=NO
169 : # Transpose the above matrix
170 : bonds_matT: TRANSPOSE ARG=bonds_mat
171 : # And multiply by the volume
172 : bonds: MATRIX_VECTOR_PRODUCT ARG=bonds_matT,sphere
173 : # And make sure we calculate the coordination numbers for the atoms in the volume too
174 : mask: CUSTOM ARG=bonds,sphere FUNC=x+y PERIODIC=NO
175 : # Calculate the coordination number for all the atoms
176 : cc: COORDINATIONNUMBER SPECIES=1-400 SWITCH={RATIONAL D_0=1.3 R_0=0.2 D_MAX=3.0} MASK=mask
177 : # Calculate the local average
178 : lcc: LOCAL_AVERAGE SPECIES=cc SWITCH={RATIONAL D_0=1.3 R_0=0.2 D_MAX=3.0} MASK=sphere
179 : # Multiply local averages parameters numbers by sphere vector
180 : prod: CUSTOM ARG=lcc,sphere FUNC=x*y PERIODIC=NO
181 : # Sum of coordination numbers for atoms that are in the sphere of interest
182 : numer: SUM ARG=prod PERIODIC=NO
183 : # Number of atoms that are in sphere of interest
184 : denom: SUM ARG=sphere PERIODIC=NO
185 : # Average coordination number for atoms in sphere of interest
186 : av: CUSTOM ARG=numer,denom FUNC=x/y PERIODIC=NO
187 : # And print out final CV to a file
188 : PRINT ARG=av FILE=colvar STRIDE=1
189 : ```
190 :
191 : In the above input the only elements of the [CONTACT_MATRIX](CONTACT_MATRIX.md) with label `vol_cmap` that are non-zero are those that are bonded to one of the atoms
192 : in the sphere of interest. By transposing this matrix and multiplying by a vector of ones we thus get a vector with a length equal to total number of coordination
193 : numbers in the system. The only elements in this vector that are non-zero are those that are bonded to atoms that are within the sphere of interest. Consequently,
194 : if we add the vector `sphere` to this we can a vector in which element $i$ is only non-zero if we have to calculate that coordination number to evaluate the final
195 : quantity of interest.
196 :
197 : */
198 : //+ENDPLUMEDOC
199 :
200 : class LocalAverage : public ActionShortcut {
201 : private:
202 : std::string getMomentumSymbol( const int& m ) const ;
203 : public:
204 : static void registerKeywords( Keywords& keys );
205 : explicit LocalAverage(const ActionOptions&);
206 : };
207 :
208 : PLUMED_REGISTER_ACTION(LocalAverage,"LOCAL_AVERAGE")
209 :
210 15 : void LocalAverage::registerKeywords( Keywords& keys ) {
211 15 : CoordinationNumbers::shortcutKeywords( keys );
212 15 : keys.needsAction("ONES");
213 15 : keys.needsAction("MATRIX_VECTOR_PRODUCT");
214 15 : keys.needsAction("VSTACK");
215 15 : keys.needsAction("CUSTOM");
216 15 : keys.needsAction("OUTER_PRODUCT");
217 30 : keys.setValueDescription("vector","the values of the local averages");
218 30 : keys.addDOI("10.1063/1.2977970");
219 15 : keys.addFlag("USEGPU",false,"run part of this calculation on the GPU");
220 15 : }
221 :
222 11 : LocalAverage::LocalAverage(const ActionOptions&ao):
223 : Action(ao),
224 11 : ActionShortcut(ao) {
225 : std::string sp_str, specA, specB;
226 11 : parse("SPECIES",sp_str);
227 11 : parse("SPECIESA",specA);
228 11 : parse("SPECIESB",specB);
229 :
230 : bool usegpu;
231 11 : parseFlag("USEGPU",usegpu);
232 22 : const std::string doUSEGPU = usegpu?" USEGPU":"";
233 :
234 11 : CoordinationNumbers::expandMatrix( false,
235 : getShortcutLabel(),
236 : sp_str,
237 : specA,
238 : specB,
239 : this );
240 : std::map<std::string,std::string> keymap;
241 11 : multicolvar::MultiColvarShortcuts::readShortcutKeywords( keymap, this );
242 11 : if( sp_str.length()>0 ) {
243 : specA=specB=sp_str;
244 : }
245 11 : const std::string matLab=getShortcutLabel() + "_mat";
246 11 : const std::string onesLab=getShortcutLabel() + "_ones";
247 11 : const std::string coordLab=getShortcutLabel() + "_coord";
248 : // Calculate the coordination numbers
249 11 : ActionWithValue* av = plumed.getActionSet().selectWithLabel<ActionWithValue*>(matLab);
250 11 : plumed_assert( av && av->getNumberOfComponents()>0 && (av->copyOutput(0))->getRank()==2 );
251 : std::string size;
252 11 : Tools::convert( (av->copyOutput(0))->getShape()[1], size );
253 11 : readInputLine(onesLab + ": ONES SIZE=" + size );
254 22 : readInputLine(coordLab + ": MATRIX_VECTOR_PRODUCT ARG=" + matLab + "," + onesLab + doUSEGPU);
255 :
256 11 : int l=-1;
257 11 : std::vector<ActionShortcut*> shortcuts=plumed.getActionSet().select<ActionShortcut*>();
258 104 : for(unsigned i=0; i<shortcuts.size(); ++i) {
259 102 : if( specA==shortcuts[i]->getShortcutLabel() ) {
260 22 : std::string sname = shortcuts[i]->getName();
261 22 : if( sname=="Q1" || sname=="Q3" || sname=="Q4" || sname=="Q6" ) {
262 18 : Tools::convert( sname.substr(1), l );
263 : break;
264 : }
265 : }
266 : }
267 :
268 11 : const std::string prodLab=getShortcutLabel() + "_prod";
269 11 : if( l>0 ) {
270 9 : std::string vargs="";
271 9 : std::string svargs="";
272 9 : std::string comma="ARG=";
273 9 : std::string sargs = "ARG=" + matLab;
274 194 : auto stringDiv=[&](const std::string& num, char realImm) {
275 388 : return specB + "_"+realImm+"mn-" + num + ": CUSTOM ARG="
276 582 : + specB + "_sp."+realImm+"m-" + num + ","
277 582 : + specB + "_denom FUNC=x/y PERIODIC=NO";
278 9 : };
279 106 : for(int i=-l; i<=l; ++i) {
280 97 : std::string num = getMomentumSymbol(i);
281 194 : if( !plumed.getActionSet().selectWithLabel<ActionWithValue*>(specB + "_rmn-" + num) ) {
282 194 : readInputLine( stringDiv(num,'r') );
283 : }
284 194 : if( !plumed.getActionSet().selectWithLabel<ActionWithValue*>(specB + "_imn-" + num) ) {
285 194 : readInputLine( stringDiv(num,'i') );
286 : }
287 291 : vargs += comma + specB + "_rmn-" + num + "," + specB + "_imn-" + num;
288 194 : svargs += comma + prodLab + "." + specB + "_rmn-" + num
289 388 : + "," + prodLab + "." + specB + "_imn-" + num;
290 : comma=",";
291 485 : sargs += "," + specB + "_rmn-" + num + "," + specB + "_imn-" + num;
292 : }
293 :
294 9 : const std::string vstackLab=getShortcutLabel() + "_vstack";
295 9 : const std::string vpstackLab=getShortcutLabel() + "_vpstack";
296 9 : readInputLine(vstackLab + ": VSTACK " + vargs );
297 9 : readInputLine(prodLab + ": MATRIX_VECTOR_PRODUCT " + sargs + doUSEGPU);
298 18 : readInputLine(vpstackLab + ": VSTACK " + svargs );
299 : std::string twolplusone;
300 9 : Tools::convert( 2*(2*l+1), twolplusone );
301 :
302 9 : const std::string lonesLab=getShortcutLabel() + "_lones";
303 9 : readInputLine(lonesLab + ": ONES SIZE=" + twolplusone );
304 :
305 9 : const std::string unormLab=getShortcutLabel() + "_unorm";
306 18 : readInputLine(unormLab + ": OUTER_PRODUCT ARG=" +coordLab + "," + lonesLab );
307 :
308 9 : const std::string avLab=getShortcutLabel() + "_av";
309 27 : readInputLine(avLab + ": CUSTOM ARG=" + vpstackLab + "," + vstackLab + ","
310 9 : + unormLab + " FUNC=(x+y)/(1+z) PERIODIC=NO");
311 :
312 9 : const std::string av2Lab=getShortcutLabel() + "_av2";
313 18 : readInputLine(av2Lab + ": CUSTOM ARG=" +avLab + " FUNC=x*x PERIODIC=NO");
314 :
315 9 : const std::string twoLab=getShortcutLabel() + "_2";
316 27 : readInputLine( twoLab+": MATRIX_VECTOR_PRODUCT ARG=" +av2Lab + "," + lonesLab
317 9 : + doUSEGPU);
318 18 : readInputLine( getShortcutLabel() + ": CUSTOM "
319 18 : "ARG=" + twoLab + " FUNC=sqrt(x) PERIODIC=NO");
320 : } else {
321 4 : readInputLine( prodLab + ": MATRIX_VECTOR_PRODUCT "
322 4 : "ARG=" +matLab + ","+ specB
323 2 : + doUSEGPU);
324 4 : readInputLine( getShortcutLabel() + ": CUSTOM "
325 6 : "ARG=" + prodLab + "," + specA + "," +coordLab
326 2 : + " FUNC=(x+y)/(1+z) PERIODIC=NO");
327 : }
328 22 : multicolvar::MultiColvarShortcuts::expandFunctions( getShortcutLabel(),
329 : getShortcutLabel(),
330 : "",
331 : keymap, this );
332 205 : }
333 :
334 97 : std::string LocalAverage::getMomentumSymbol( const int& m ) const {
335 97 : if( m<0 ) {
336 : std::string num;
337 44 : Tools::convert( -1*m, num );
338 44 : return "n" + num;
339 53 : } else if( m>0 ) {
340 : std::string num;
341 44 : Tools::convert( m, num );
342 44 : return "p" + num;
343 : }
344 9 : return "0";
345 : }
346 :
347 : }
348 : }
|
