Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2015-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 "multicolvar/MultiColvarShortcuts.h"
25 : #include "core/PlumedMain.h"
26 : #include "core/ActionSet.h"
27 : #include "core/ActionWithValue.h"
28 :
29 : namespace PLMD {
30 : namespace symfunc {
31 :
32 : //+PLUMEDOC MCOLVARF LOCAL_Q1
33 : /*
34 : Calculate the local degree of order around an atoms by taking the average dot product between the q_1 vector on the central atom and the q_3 vector on the atoms in the first coordination sphere.
35 :
36 : The [Q1](Q1.md) command allows one to calculate one complex vectors for each of the atoms in your system that describe the degree of order in the coordination sphere
37 : around a particular atom. The difficulty with these vectors comes when combining the order parameters from all of the individual atoms/molecules so as to get a
38 : measure of the global degree of order for the system. The simplest way of doing this - calculating the average Steinhardt parameter - can be problematic. If one is
39 : examining nucleation say only the order parameters for those atoms in the nucleus will change significantly when the nucleus forms. The order parameters for the
40 : atoms in the surrounding liquid will remain pretty much the same. As such if one models a small nucleus embedded in a very large amount of solution/melt any
41 : change in the average order parameter will be negligible. Substantial changes in the value of this average can be observed in simulations of nucleation but only
42 : because the number of atoms is relatively small.
43 :
44 : When the average [Q1](Q1.md) parameter is used to bias the dynamics a problems
45 : can occur. These averaged coordinates cannot distinguish between the correct,
46 : single-nucleus pathway and a concerted pathway in which all the atoms rearrange
47 : themselves into their solid-like configuration simultaneously. This second type
48 : of pathway would be impossible in reality because there is a large entropic
49 : barrier that prevents concerted processes like this from happening. However,
50 : in the finite sized systems that are commonly simulated this barrier is reduced
51 : substantially. As a result in simulations where average Steinhardt parameters
52 : are biased there are often quite dramatic system size effects
53 :
54 : If one wants to simulate nucleation using some form on
55 : biased dynamics what is really required is an order parameter that measures:
56 :
57 : - Whether or not the coordination spheres around atoms are ordered
58 : - Whether or not the atoms that are ordered are clustered together in a crystalline nucleus
59 :
60 : As discussed on [this page](https://www.plumed-tutorials.org/lessons/23/001/data/Steinhardt.html) a variety of variations on the Steinhardt parameters have been
61 : introduced to better describe nucleation. That page also shows how PLUMED provides you with flexibility that you can use to implement new combinations of the
62 : Steinhardt parameters. However, the inputs that you would need to write to implement common symmetry functions are rather complex so we also provide shortcuts
63 : like this one to help you compute CVs that have been widely used in the literature easily.
64 :
65 : This particular shortcut allows you to compute the LOCAL_Q1 parameter for a particular atom, which is a number that measures the extent to
66 : which the orientation of the atoms in the first coordination sphere of an atom match the orientation of the central atom. It does this by calculating the following
67 : quantity for each of the atoms in the system:
68 :
69 : $$
70 : s_i = \frac{ \sum_j \sigma( r_{ij} ) \sum_{m=-1}^1 q_{1m}^{*}(i)q_{1m}(j) }{ \sum_j \sigma( r_{ij} ) }
71 : $$
72 :
73 : where $q_{1m}(i)$ and $q_{1m}(j)$ are the 1st order Steinhardt vectors calculated for atom $i$ and atom $j$ respectively and the asterisk denotes complex
74 : conjugation. The function
75 : $\sigma( r_{ij} )$ is a switching function that acts on the distance between atoms $i$ and $j$. The parameters of this function should be set
76 : so that it the function is equal to one when atom $j$ is in the first coordination sphere of atom $i$ and is zero otherwise. The sum in the numerator
77 : of this expression is the dot product of the Steinhardt parameters for atoms $i$ and $j$ and thus measures the degree to which the orientations of these
78 : adjacent atoms are correlated.
79 :
80 : The following input shows how this works in practice. This input calculates the average value of the LOCAL_Q1 parameter for the 64 Lennard Jones atoms in the system under study and prints this
81 : quantity to a file called colvar.
82 :
83 : ```plumed
84 : q1: Q1 SPECIES=1-64 D_0=1.3 R_0=0.2
85 : lq1: LOCAL_Q1 SPECIES=q1 SWITCH={RATIONAL D_0=1.3 R_0=0.2} MEAN
86 : PRINT ARG=lq1.mean FILE=colvar
87 : ```
88 :
89 : The following input calculates the distribution of LOCAL_Q1 parameters at any given time and outputs this information to a file.
90 :
91 : ```plumed
92 : q1: Q1 SPECIES=1-64 D_0=1.3 R_0=0.2
93 : lq1: LOCAL_Q1 SPECIES=q1 SWITCH={RATIONAL D_0=1.3 R_0=0.2} HISTOGRAM={GAUSSIAN LOWER=0.0 UPPER=1.0 NBINS=20 SMEAR=0.1}
94 : PRINT ARG=lq1.* FILE=colvar
95 : ```
96 :
97 : The following calculates the LOCAL_Q1 parameters for atoms 1-5 only. For each of these atoms comparisons of the geometry of the coordination sphere
98 : are done with those of all the other atoms in the system. The final quantity is the average and is outputted to a file
99 :
100 : ```plumed
101 : q1a: Q1 SPECIESA=1-5 SPECIESB=1-64 D_0=1.3 R_0=0.2
102 : q1b: Q1 SPECIESA=6-64 SPECIESB=1-64 D_0=1.3 R_0=0.2
103 : w1: LOCAL_Q1 SPECIES=q1a,q1b SWITCH={RATIONAL D_0=1.3 R_0=0.2} MEAN
104 : PRINT ARG=w1.* FILE=colvar
105 : ```
106 :
107 : */
108 : //+ENDPLUMEDOC
109 :
110 : //+PLUMEDOC MCOLVARF LOCAL_Q3
111 : /*
112 : Calculate the local degree of order around an atoms by taking the average dot product between the q_3 vector on the central atom and the q_3 vector on the atoms in the first coordination sphere.
113 :
114 : The [Q3](Q3.md) command allows one to calculate one complex vectors for each of the atoms in your system that describe the degree of order in the coordination sphere
115 : around a particular atom. The difficulty with these vectors comes when combining the order parameters from all of the individual atoms/molecules so as to get a
116 : measure of the global degree of order for the system. The simplest way of doing this - calculating the average Steinhardt parameter - can be problematic. If one is
117 : examining nucleation say only the order parameters for those atoms in the nucleus will change significantly when the nucleus forms. The order parameters for the
118 : atoms in the surrounding liquid will remain pretty much the same. As such if one models a small nucleus embedded in a very large amount of solution/melt any
119 : change in the average order parameter will be negligible. Substantial changes in the value of this average can be observed in simulations of nucleation but only
120 : because the number of atoms is relatively small.
121 :
122 : When the average [Q3](Q3.md) parameter is used to bias the dynamics a problems
123 : can occur. These averaged coordinates cannot distinguish between the correct,
124 : single-nucleus pathway and a concerted pathway in which all the atoms rearrange
125 : themselves into their solid-like configuration simultaneously. This second type
126 : of pathway would be impossible in reality because there is a large entropic
127 : barrier that prevents concerted processes like this from happening. However,
128 : in the finite sized systems that are commonly simulated this barrier is reduced
129 : substantially. As a result in simulations where average Steinhardt parameters
130 : are biased there are often quite dramatic system size effects
131 :
132 : If one wants to simulate nucleation using some form on
133 : biased dynamics what is really required is an order parameter that measures:
134 :
135 : - Whether or not the coordination spheres around atoms are ordered
136 : - Whether or not the atoms that are ordered are clustered together in a crystalline nucleus
137 :
138 : As discussed on [this page](https://www.plumed-tutorials.org/lessons/23/001/data/Steinhardt.html) a variety of variations on the Steinhardt parameters have been
139 : introduced to better describe nucleation. That page also shows how PLUMED provides you with flexibility that you can use to implement new combinations of the
140 : Steinhardt parameters. However, the inputs that you would need to write to implement common symmetry functions are rather complex so we also provide shortcuts
141 : like this one to help you compute CVs that have been widely used in the literature easily.
142 :
143 : This particular shortcut allows you to compute the LOCAL_Q3 parameter for a particular atom, which is a number that measures the extent to
144 : which the orientation of the atoms in the first coordination sphere of an atom match the orientation of the central atom. It does this by calculating the following
145 : quantity for each of the atoms in the system:
146 :
147 : $$
148 : s_i = \frac{ \sum_j \sigma( r_{ij} ) \sum_{m=-3}^3 q_{3m}^{*}(i)q_{3m}(j) }{ \sum_j \sigma( r_{ij} ) }
149 : $$
150 :
151 : where $q_{3m}(i)$ and $q_{3m}(j)$ are the 3rd order Steinhardt vectors calculated for atom $i$ and atom $j$ respectively and the asterisk denotes complex
152 : conjugation. The function
153 : $\sigma( r_{ij} )$ is a switching function that acts on the distance between atoms $i$ and $j$. The parameters of this function should be set
154 : so that it the function is equal to one when atom $j$ is in the first coordination sphere of atom $i$ and is zero otherwise. The sum in the numerator
155 : of this expression is the dot product of the Steinhardt parameters for atoms $i$ and $j$ and thus measures the degree to which the orientations of these
156 : adjacent atoms are correlated.
157 :
158 : The following input shows how this works in practice. This input calculates the average value of the LOCAL_Q3 parameter for the 64 Lennard Jones atoms in the system under study and prints this
159 : quantity to a file called colvar.
160 :
161 : ```plumed
162 : q3: Q3 SPECIES=1-64 D_0=1.3 R_0=0.2
163 : lq3: LOCAL_Q3 SPECIES=q3 SWITCH={RATIONAL D_0=1.3 R_0=0.2} MEAN
164 : PRINT ARG=lq3.mean FILE=colvar
165 : ```
166 :
167 : The following input calculates the distribution of LOCAL_Q3 parameters at any given time and outputs this information to a file.
168 :
169 : ```plumed
170 : q3: Q3 SPECIES=1-64 D_0=1.3 R_0=0.2
171 : lq3: LOCAL_Q3 SPECIES=q3 SWITCH={RATIONAL D_0=1.3 R_0=0.2} HISTOGRAM={GAUSSIAN LOWER=0.0 UPPER=1.0 NBINS=20 SMEAR=0.1}
172 : PRINT ARG=lq3.* FILE=colvar
173 : ```
174 :
175 : The following calculates the LOCAL_Q3 parameters for atoms 1-5 only. For each of these atoms comparisons of the geometry of the coordination sphere
176 : are done with those of all the other atoms in the system. The final quantity is the average and is outputted to a file
177 :
178 : ```plumed
179 : q3a: Q3 SPECIESA=1-5 SPECIESB=1-64 D_0=1.3 R_0=0.2
180 : q3b: Q3 SPECIESA=6-64 SPECIESB=1-64 D_0=1.3 R_0=0.2
181 : w3: LOCAL_Q3 SPECIES=q3a,q3b SWITCH={RATIONAL D_0=1.3 R_0=0.2} MEAN
182 : PRINT ARG=w3.* FILE=colvar
183 : ```
184 :
185 : */
186 : //+ENDPLUMEDOC
187 :
188 : //+PLUMEDOC MCOLVARF LOCAL_Q4
189 : /*
190 : Calculate the local degree of order around an atoms by taking the average dot product between the q_4 vector on the central atom and the q_4 vector on the atoms in the first coordination sphere.
191 :
192 : The [Q4](Q4.md) command allows one to calculate one complex vectors for each of the atoms in your system that describe the degree of order in the coordination sphere
193 : around a particular atom. The difficulty with these vectors comes when combining the order parameters from all of the individual atoms/molecules so as to get a
194 : measure of the global degree of order for the system. The simplest way of doing this - calculating the average Steinhardt parameter - can be problematic. If one is
195 : examining nucleation say only the order parameters for those atoms in the nucleus will change significantly when the nucleus forms. The order parameters for the
196 : atoms in the surrounding liquid will remain pretty much the same. As such if one models a small nucleus embedded in a very large amount of solution/melt any
197 : change in the average order parameter will be negligible. Substantial changes in the value of this average can be observed in simulations of nucleation but only
198 : because the number of atoms is relatively small.
199 :
200 : When the average [Q4](Q4.md) parameter is used to bias the dynamics a problems
201 : can occur. These averaged coordinates cannot distinguish between the correct,
202 : single-nucleus pathway and a concerted pathway in which all the atoms rearrange
203 : themselves into their solid-like configuration simultaneously. This second type
204 : of pathway would be impossible in reality because there is a large entropic
205 : barrier that prevents concerted processes like this from happening. However,
206 : in the finite sized systems that are commonly simulated this barrier is reduced
207 : substantially. As a result in simulations where average Steinhardt parameters
208 : are biased there are often quite dramatic system size effects
209 :
210 : If one wants to simulate nucleation using some form on
211 : biased dynamics what is really required is an order parameter that measures:
212 :
213 : - Whether or not the coordination spheres around atoms are ordered
214 : - Whether or not the atoms that are ordered are clustered together in a crystalline nucleus
215 :
216 : As discussed on [this page](https://www.plumed-tutorials.org/lessons/23/001/data/Steinhardt.html) a variety of variations on the Steinhardt parameters have been
217 : introduced to better describe nucleation. That page also shows how PLUMED provides you with flexibility that you can use to implement new combinations of the
218 : Steinhardt parameters. However, the inputs that you would need to write to implement common symmetry functions are rather complex so we also provide shortcuts
219 : like this one to help you compute CVs that have been widely used in the literature easily.
220 :
221 : This particular shortcut allows you to compute the LOCAL_Q4 parameter for a particular atom, which is a number that measures the extent to
222 : which the orientation of the atoms in the first coordination sphere of an atom match the orientation of the central atom. It does this by calculating the following
223 : quantity for each of the atoms in the system:
224 :
225 : $$
226 : s_i = \frac{ \sum_j \sigma( r_{ij} ) \sum_{m=-4}^4 q_{4m}^{*}(i)q_{4m}(j) }{ \sum_j \sigma( r_{ij} ) }
227 : $$
228 :
229 : where $q_{4m}(i)$ and $q_{4m}(j)$ are the 4th order Steinhardt vectors calculated for atom $i$ and atom $j$ respectively and the asterisk denotes complex
230 : conjugation. The function
231 : $\sigma( r_{ij} )$ is a switching function that acts on the distance between atoms $i$ and $j$. The parameters of this function should be set
232 : so that it the function is equal to one when atom $j$ is in the first coordination sphere of atom $i$ and is zero otherwise. The sum in the numerator
233 : of this expression is the dot product of the Steinhardt parameters for atoms $i$ and $j$ and thus measures the degree to which the orientations of these
234 : adjacent atoms are correlated.
235 :
236 : The following input shows how this works in practice. This input calculates the average value of the LOCAL_Q4 parameter for the 64 Lennard Jones atoms in the system under study and prints this
237 : quantity to a file called colvar.
238 :
239 : ```plumed
240 : q4: Q4 SPECIES=1-64 D_0=1.3 R_0=0.2
241 : lq4: LOCAL_Q4 SPECIES=q4 SWITCH={RATIONAL D_0=1.3 R_0=0.2} MEAN
242 : PRINT ARG=lq4.mean FILE=colvar
243 : ```
244 :
245 : The following input calculates the distribution of LOCAL_Q4 parameters at any given time and outputs this information to a file.
246 :
247 : ```plumed
248 : q4: Q4 SPECIES=1-64 D_0=1.3 R_0=0.2
249 : lq4: LOCAL_Q4 SPECIES=q4 SWITCH={RATIONAL D_0=1.3 R_0=0.2} HISTOGRAM={GAUSSIAN LOWER=0.0 UPPER=1.0 NBINS=20 SMEAR=0.1}
250 : PRINT ARG=lq4.* FILE=colvar
251 : ```
252 :
253 : The following calculates the LOCAL_Q4 parameters for atoms 1-5 only. For each of these atoms comparisons of the geometry of the coordination sphere
254 : are done with those of all the other atoms in the system. The final quantity is the average and is outputted to a file
255 :
256 : ```plumed
257 : q4a: Q4 SPECIESA=1-5 SPECIESB=1-64 D_0=1.3 R_0=0.2
258 : q4b: Q4 SPECIESA=6-64 SPECIESB=1-64 D_0=1.3 R_0=0.2
259 : w4: LOCAL_Q4 SPECIES=q4a,q4b SWITCH={RATIONAL D_0=1.3 R_0=0.2} MEAN
260 : PRINT ARG=w4.* FILE=colvar
261 : ```
262 :
263 : */
264 : //+ENDPLUMEDOC
265 :
266 : //+PLUMEDOC MCOLVARF LOCAL_Q6
267 : /*
268 : Calculate the local degree of order around an atoms by taking the average dot product between the q_6 vector on the central atom and the q_6 vector on the atoms in the first coordination sphere.
269 :
270 : The [Q6](Q6.md) command allows one to calculate one complex vectors for each of the atoms in your system that describe the degree of order in the coordination sphere
271 : around a particular atom. The difficulty with these vectors comes when combining the order parameters from all of the individual atoms/molecules so as to get a
272 : measure of the global degree of order for the system. The simplest way of doing this - calculating the average Steinhardt parameter - can be problematic. If one is
273 : examining nucleation say only the order parameters for those atoms in the nucleus will change significantly when the nucleus forms. The order parameters for the
274 : atoms in the surrounding liquid will remain pretty much the same. As such if one models a small nucleus embedded in a very large amount of solution/melt any
275 : change in the average order parameter will be negligible. Substantial changes in the value of this average can be observed in simulations of nucleation but only
276 : because the number of atoms is relatively small.
277 :
278 : When the average [Q6](Q6.md) parameter is used to bias the dynamics a problems
279 : can occur. These averaged coordinates cannot distinguish between the correct,
280 : single-nucleus pathway and a concerted pathway in which all the atoms rearrange
281 : themselves into their solid-like configuration simultaneously. This second type
282 : of pathway would be impossible in reality because there is a large entropic
283 : barrier that prevents concerted processes like this from happening. However,
284 : in the finite sized systems that are commonly simulated this barrier is reduced
285 : substantially. As a result in simulations where average Steinhardt parameters
286 : are biased there are often quite dramatic system size effects
287 :
288 : If one wants to simulate nucleation using some form on
289 : biased dynamics what is really required is an order parameter that measures:
290 :
291 : - Whether or not the coordination spheres around atoms are ordered
292 : - Whether or not the atoms that are ordered are clustered together in a crystalline nucleus
293 :
294 : As discussed on [this page](https://www.plumed-tutorials.org/lessons/23/001/data/Steinhardt.html) a variety of variations on the Steinhardt parameters have been
295 : introduced to better describe nucleation. That page also shows how PLUMED provides you with flexibility that you can use to implement new combinations of the
296 : Steinhardt parameters. However, the inputs that you would need to write to implement common symmetry functions are rather complex so we also provide shortcuts
297 : like this one to help you compute CVs that have been widely used in the literature easily.
298 :
299 : This particular shortcut allows you to compute the LOCAL_Q6 parameter for a particular atom, which is a number that measures the extent to
300 : which the orientation of the atoms in the first coordination sphere of an atom match the orientation of the central atom. It does this by calculating the following
301 : quantity for each of the atoms in the system:
302 :
303 : $$
304 : s_i = \frac{ \sum_j \sigma( r_{ij} ) \sum_{m=-6}^6 q_{6m}^{*}(i)q_{6m}(j) }{ \sum_j \sigma( r_{ij} ) }
305 : $$
306 :
307 : where $q_{6m}(i)$ and $q_{6m}(j)$ are the 6th order Steinhardt vectors calculated for atom $i$ and atom $j$ respectively and the asterisk denotes complex
308 : conjugation. The function
309 : $\sigma( r_{ij} )$ is a switching function that acts on the distance between atoms $i$ and $j$. The parameters of this function should be set
310 : so that it the function is equal to one when atom $j$ is in the first coordination sphere of atom $i$ and is zero otherwise. The sum in the numerator
311 : of this expression is the dot product of the Steinhardt parameters for atoms $i$ and $j$ and thus measures the degree to which the orientations of these
312 : adjacent atoms are correlated.
313 :
314 : The following input shows how this works in practice. This input calculates the average value of the LOCAL_Q6 parameter for the 64 Lennard Jones atoms in the system under study and prints this
315 : quantity to a file called colvar.
316 :
317 : ```plumed
318 : q6: Q6 SPECIES=1-64 D_0=1.3 R_0=0.2
319 : lq6: LOCAL_Q6 SPECIES=q6 SWITCH={RATIONAL D_0=1.3 R_0=0.2} MEAN
320 : PRINT ARG=lq6.mean FILE=colvar
321 : ```
322 :
323 : The following input calculates the distribution of LOCAL_Q6 parameters at any given time and outputs this information to a file.
324 :
325 : ```plumed
326 : q6: Q6 SPECIES=1-64 D_0=1.3 R_0=0.2
327 : lq6: LOCAL_Q6 SPECIES=q6 SWITCH={RATIONAL D_0=1.3 R_0=0.2} HISTOGRAM={GAUSSIAN LOWER=0.0 UPPER=1.0 NBINS=20 SMEAR=0.1}
328 : PRINT ARG=lq6.* FILE=colvar
329 : ```
330 :
331 : The following calculates the LOCAL_Q6 parameters for atoms 1-5 only. For each of these atoms comparisons of the geometry of the coordination sphere
332 : are done with those of all the other atoms in the system. The final quantity is the average and is outputted to a file
333 :
334 : ```plumed
335 : q6a: Q6 SPECIESA=1-5 SPECIESB=1-64 D_0=1.3 R_0=0.2
336 : q6b: Q6 SPECIESA=6-64 SPECIESB=1-64 D_0=1.3 R_0=0.2
337 : w6: LOCAL_Q6 SPECIES=q6a,q6b SWITCH={RATIONAL D_0=1.3 R_0=0.2} MEAN
338 : PRINT ARG=w6.* FILE=colvar
339 : ```
340 :
341 : */
342 : //+ENDPLUMEDOC
343 :
344 : class LocalSteinhardt : public ActionShortcut {
345 : private:
346 : std::string getSymbol( const int& m ) const ;
347 : std::string getArgsForStack( const int& l, const std::string& lab ) const;
348 : public:
349 : static void registerKeywords( Keywords& keys );
350 : explicit LocalSteinhardt(const ActionOptions&);
351 : };
352 :
353 : PLUMED_REGISTER_ACTION(LocalSteinhardt,"LOCAL_Q1")
354 : PLUMED_REGISTER_ACTION(LocalSteinhardt,"LOCAL_Q3")
355 : PLUMED_REGISTER_ACTION(LocalSteinhardt,"LOCAL_Q4")
356 : PLUMED_REGISTER_ACTION(LocalSteinhardt,"LOCAL_Q6")
357 :
358 20 : void LocalSteinhardt::registerKeywords( Keywords& keys ) {
359 20 : ActionShortcut::registerKeywords( keys );
360 20 : keys.add("optional","SPECIES","the label of the action that computes the Steinhardt parameters for which you would like to calculate local steinhardt parameters");
361 20 : keys.add("optional","SPECIESA","the label of the action that computes the Steinhardt parameters for which you would like to calculate local steinhardt parameters");
362 20 : keys.add("optional","SPECIESB","the label of the action that computes the Steinhardt parameters that you would like to use when calculating the loal steinhardt parameters");
363 20 : keys.add("optional","SWITCH","This keyword is used if you want to employ an alternative to the continuous swiching function defined above. "
364 : "The following provides information on the \\ref switchingfunction that are available. "
365 : "When this keyword is present you no longer need the NN, MM, D_0 and R_0 keywords.");
366 20 : keys.addFlag("LOWMEM",false,"this flag does nothing and is present only to ensure back-compatibility");
367 40 : keys.setValueDescription("vector","the values of the local steinhardt parameters for the input atoms");
368 20 : multicolvar::MultiColvarShortcuts::shortcutKeywords( keys );
369 20 : keys.needsAction("CONTACT_MATRIX");
370 20 : keys.needsAction("MATRIX_PRODUCT");
371 20 : keys.needsAction("GROUP");
372 20 : keys.needsAction("ONES");
373 20 : keys.needsAction("OUTER_PRODUCT");
374 20 : keys.needsAction("VSTACK");
375 20 : keys.needsAction("CONCATENATE");
376 20 : keys.needsAction("CUSTOM");
377 20 : keys.needsAction("TRANSPOSE");
378 20 : keys.needsAction("MATRIX_VECTOR_PRODUCT");
379 20 : }
380 :
381 98 : std::string LocalSteinhardt::getSymbol( const int& m ) const {
382 98 : if( m<0 ) {
383 : std::string num;
384 40 : Tools::convert( -1*m, num );
385 40 : return "n" + num;
386 58 : } else if( m>0 ) {
387 : std::string num;
388 49 : Tools::convert( m, num );
389 49 : return "p" + num;
390 : }
391 9 : return "0";
392 : }
393 :
394 9 : std::string LocalSteinhardt::getArgsForStack( const int& l, const std::string& sp_lab ) const {
395 : std::string numstr;
396 9 : Tools::convert( l, numstr );
397 18 : std::string data_mat = " ARG=" + sp_lab + "_sp.rm-n" + numstr + "," + sp_lab + "_sp.im-n" + numstr;
398 107 : for(int i=-l+1; i<=l; ++i) {
399 98 : numstr = getSymbol( i );
400 196 : data_mat += "," + sp_lab + "_sp.rm-" + numstr + "," + sp_lab + "_sp.im-" + numstr;
401 : }
402 9 : return data_mat;
403 : }
404 :
405 7 : LocalSteinhardt::LocalSteinhardt(const ActionOptions& ao):
406 : Action(ao),
407 7 : ActionShortcut(ao) {
408 : bool lowmem;
409 7 : parseFlag("LOWMEM",lowmem);
410 7 : if( lowmem ) {
411 0 : warning("LOWMEM flag is deprecated and is no longer required for this action");
412 : }
413 : // Get the Q value
414 : int l;
415 14 : Tools::convert( getName().substr(7), l);
416 : // Create a vector filled with ones
417 : std::string twolplusone;
418 7 : Tools::convert( 2*(2*l+1), twolplusone );
419 14 : readInputLine( getShortcutLabel() + "_uvec: ONES SIZE=" + twolplusone );
420 : // Read in species keyword
421 : std::string sp_str;
422 14 : parse("SPECIES",sp_str);
423 : std::string spa_str;
424 14 : parse("SPECIESA",spa_str);
425 7 : if( sp_str.length()>0 ) {
426 : // Create a group with these atoms
427 12 : readInputLine( getShortcutLabel() + "_grp: GROUP ATOMS=" + sp_str );
428 6 : std::vector<std::string> sp_lab = Tools::getWords(sp_str, "\t\n ,");
429 : // This creates the stash to hold all the vectors
430 6 : if( sp_lab.size()==1 ) {
431 : // The lengths of all the vectors in a vector
432 12 : readInputLine( getShortcutLabel() + "_nmat: OUTER_PRODUCT ARG=" + sp_lab[0] + "_norm," + getShortcutLabel() + "_uvec");
433 : // The unormalised vectors
434 12 : readInputLine( getShortcutLabel() + "_uvecs: VSTACK" + getArgsForStack( l, sp_lab[0] ) );
435 : } else {
436 0 : std::string len_vec = getShortcutLabel() + "_mags: CONCATENATE ARG=" + sp_lab[0] + "_norm";
437 0 : for(unsigned i=1; i<sp_lab.size(); ++i) {
438 0 : len_vec += "," + sp_lab[i] + "_norm";
439 : }
440 : // This is the vector that contains all the magnitudes
441 0 : readInputLine( len_vec );
442 0 : std::string concat_str = getShortcutLabel() + "_uvecs: CONCATENATE";
443 0 : for(unsigned i=0; i<sp_lab.size(); ++i) {
444 : std::string snum;
445 0 : Tools::convert( i+1, snum );
446 0 : concat_str += " MATRIX" + snum + "1=" + getShortcutLabel() + "_uvecs" + snum;
447 0 : readInputLine( getShortcutLabel() + "_uvecs" + snum + ": VSTACK" + getArgsForStack( l, sp_lab[i] ) );
448 : }
449 : // And the normalising matrix by taking the column vector of magnitudes and multiplying by the row vector of ones
450 0 : readInputLine( getShortcutLabel() + "_nmat: OUTER_PRODUCT ARG=" + getShortcutLabel() + "_mags," + getShortcutLabel() + "_uvec");
451 : // The unormalised vectors
452 0 : readInputLine( concat_str );
453 : }
454 : // Now normalise all the vectors by doing Hadammard "product" with normalising matrix
455 12 : readInputLine( getShortcutLabel() + "_vecs: CUSTOM ARG=" + getShortcutLabel() + "_uvecs," + getShortcutLabel() + "_nmat FUNC=x/y PERIODIC=NO");
456 : // And transpose the matrix
457 12 : readInputLine( getShortcutLabel() + "_vecsT: TRANSPOSE ARG=" + getShortcutLabel() + "_vecs" );
458 : std::string sw_str;
459 6 : parse("SWITCH",sw_str);
460 12 : readInputLine( getShortcutLabel() + "_cmap: CONTACT_MATRIX GROUP=" + sp_str + " SWITCH={" + sw_str + "}");
461 : // And the matrix of dot products
462 12 : readInputLine( getShortcutLabel() + "_dpmat: MATRIX_PRODUCT ARG=" + getShortcutLabel() + "_vecs," + getShortcutLabel() + "_vecsT" );
463 7 : } else if( spa_str.length()>0 ) {
464 : // Create a group with these atoms
465 2 : readInputLine( getShortcutLabel() + "_grp: GROUP ATOMS=" + spa_str );
466 : std::string spb_str;
467 2 : parse("SPECIESB",spb_str);
468 1 : if( spb_str.length()==0 ) {
469 0 : plumed_merror("need both SPECIESA and SPECIESB in input");
470 : }
471 1 : std::vector<std::string> sp_laba = Tools::getWords(spa_str, "\t\n ,");
472 1 : std::vector<std::string> sp_labb = Tools::getWords(spb_str, "\t\n ,");
473 1 : if( sp_laba.size()==1 ) {
474 : // The matrix that is used for normalising
475 2 : readInputLine( getShortcutLabel() + "_nmatA: OUTER_PRODUCT ARG=" + sp_laba[0] + "_norm," + getShortcutLabel() + "_uvec");
476 : // The unormalised vectors
477 2 : readInputLine( getShortcutLabel() + "_uvecsA: VSTACK" + getArgsForStack( l, sp_laba[0] ) );
478 : } else {
479 0 : std::string len_vec = getShortcutLabel() + "_magsA: CONCATENATE ARG=" + sp_laba[0] + "_norm";
480 0 : for(unsigned i=1; i<sp_laba.size(); ++i) {
481 0 : len_vec += "," + sp_laba[i] + "_norm";
482 : }
483 : // This is the vector that contains all the magnitudes
484 0 : readInputLine( len_vec );
485 0 : std::string concat_str = getShortcutLabel() + "_uvecsA: CONCATENATE";
486 0 : for(unsigned i=0; i<sp_laba.size(); ++i) {
487 : std::string snum;
488 0 : Tools::convert( i+1, snum );
489 0 : concat_str += " MATRIX" + snum + "1=" + getShortcutLabel() + "_uvecsA" + snum;
490 0 : readInputLine( getShortcutLabel() + "_uvecsA" + snum + ": VSTACK" + getArgsForStack( l, sp_laba[i] ) );
491 : }
492 : // And the normalising matrix by taking the column vector of magnitudes and multiplying by the row vector of ones
493 0 : readInputLine( getShortcutLabel() + "_nmatA: OUTER_PRODUCT ARG=" + getShortcutLabel() + "_magsA," + getShortcutLabel() + "_uvec");
494 : // The unormalised vector
495 0 : readInputLine( concat_str );
496 : }
497 : // Now normalise all the vectors by doing Hadammard "product" with normalising matrix
498 2 : readInputLine( getShortcutLabel() + "_vecsA: CUSTOM ARG=" + getShortcutLabel() + "_uvecsA," + getShortcutLabel() + "_nmatA FUNC=x/y PERIODIC=NO");
499 : // Now do second matrix
500 1 : if( sp_labb.size()==1 ) {
501 0 : readInputLine( getShortcutLabel() + "_nmatB: OUTER_PRODUCT ARG=" + sp_laba[0] + "_norm," + getShortcutLabel() + "_uvec");
502 0 : readInputLine( getShortcutLabel() + "_uvecsBT: VSTACK" + getArgsForStack( l, sp_labb[0] ) );
503 0 : readInputLine( getShortcutLabel() + "_uvecsB: TRANSPOSE ARG=" + getShortcutLabel() + "_uvecsBT");
504 : } else {
505 2 : std::string len_vec = getShortcutLabel() + "_magsB: CONCATENATE ARG=" + sp_labb[0] + "_norm";
506 2 : for(unsigned i=1; i<sp_labb.size(); ++i) {
507 2 : len_vec += "," + sp_labb[i] + "_norm";
508 : }
509 : // This is the vector that contains all the magnitudes
510 1 : readInputLine( len_vec );
511 1 : std::string concat_str = getShortcutLabel() + "_uvecsB: CONCATENATE";
512 3 : for(unsigned i=0; i<sp_labb.size(); ++i) {
513 : std::string snum;
514 2 : Tools::convert( i+1, snum );
515 4 : concat_str += " MATRIX1" + snum + "=" + getShortcutLabel() + "_uvecsB" + snum;
516 4 : readInputLine( getShortcutLabel() + "_uvecsBT" + snum + ": VSTACK" + getArgsForStack( l, sp_labb[i] ) );
517 4 : readInputLine( getShortcutLabel() + "_uvecsB" + snum + ": TRANSPOSE ARG=" + getShortcutLabel() + "_uvecsBT" + snum );
518 : }
519 : // And the normalising matrix
520 2 : readInputLine( getShortcutLabel() + "_nmatB: OUTER_PRODUCT ARG=" + getShortcutLabel() + "_uvec," + getShortcutLabel() + "_magsB");
521 : // The unormalised vectors
522 1 : readInputLine( concat_str );
523 : }
524 : // Now normalise all the vectors by doing Hadammard "product" with normalising matrix
525 2 : readInputLine( getShortcutLabel() + "_vecsB: CUSTOM ARG=" + getShortcutLabel() + "_uvecsB," + getShortcutLabel() + "_nmatB FUNC=x/y PERIODIC=NO");
526 : std::string sw_str;
527 1 : parse("SWITCH",sw_str);
528 2 : readInputLine( getShortcutLabel() + "_cmap: CONTACT_MATRIX GROUPA=" + spa_str + " GROUPB=" + spb_str + " SWITCH={" + sw_str + "}");
529 2 : readInputLine( getShortcutLabel() + "_dpmat: MATRIX_PRODUCT ARG=" + getShortcutLabel() + "_vecsA," + getShortcutLabel() + "_vecsB" );
530 1 : }
531 :
532 : // Now create the product matrix
533 14 : readInputLine( getShortcutLabel() + "_prod: CUSTOM ARG=" + getShortcutLabel() + "_cmap," + getShortcutLabel() + "_dpmat FUNC=x*y PERIODIC=NO");
534 : // Now the sum of coordination numbers times the switching functions
535 7 : ActionWithValue* av = plumed.getActionSet().selectWithLabel<ActionWithValue*>( getShortcutLabel() + "_cmap");
536 7 : plumed_assert( av && av->getNumberOfComponents()>0 && (av->copyOutput(0))->getRank()==2 );
537 : std::string size;
538 7 : Tools::convert( (av->copyOutput(0))->getShape()[1], size );
539 14 : readInputLine( getShortcutLabel() + "_ones: ONES SIZE=" + size );
540 14 : readInputLine( getShortcutLabel() + ": MATRIX_VECTOR_PRODUCT ARG=" + getShortcutLabel() +"_prod," + getShortcutLabel() +"_ones");
541 : // And just the sum of the coordination numbers
542 14 : readInputLine( getShortcutLabel() + "_denom: MATRIX_VECTOR_PRODUCT ARG=" + getShortcutLabel() + "_cmap," + getShortcutLabel() +"_ones");
543 : // And matheval to get the final quantity
544 14 : readInputLine( getShortcutLabel() + "_av: CUSTOM ARG=" + getShortcutLabel() + "," + getShortcutLabel() + "_denom FUNC=x/y PERIODIC=NO");
545 : // And this expands everything
546 14 : multicolvar::MultiColvarShortcuts::expandFunctions( getShortcutLabel(), getShortcutLabel() + "_av", "", this );
547 7 : }
548 :
549 : }
550 : }
|
