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 vector 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}
86 : lq1_mean: MEAN ARG=lq1 PERIODIC=NO
87 : PRINT ARG=lq1_mean FILE=colvar
88 : ```
89 :
90 : The following input calculates the distribution of LOCAL_Q1 parameters at any given time and outputs this information to a file.
91 :
92 : ```plumed
93 : q1: Q1 SPECIES=1-64 D_0=1.3 R_0=0.2
94 : 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}
95 : PRINT ARG=lq1.* FILE=colvar
96 : ```
97 :
98 : 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
99 : are done with those of all the other atoms in the system. The final quantity is the average and is outputted to a file
100 :
101 : ```plumed
102 : q1a: Q1 SPECIESA=1-5 SPECIESB=1-64 D_0=1.3 R_0=0.2
103 : q1b: Q1 SPECIESA=6-64 SPECIESB=1-64 D_0=1.3 R_0=0.2
104 : w1: LOCAL_Q1 SPECIESA=q1a SPECIESB=q1b SWITCH={RATIONAL D_0=1.3 R_0=0.2} MEAN
105 : PRINT ARG=w1.* FILE=colvar
106 : ```
107 :
108 : ## The MASK keyword
109 :
110 : 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
111 : 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
112 : from the [volumes module](module_volumes.md) to calculate the average value of the $s_i$ parameter that is defined by the equation above for only those atoms that
113 : lie in a certain part of the simulation box.
114 :
115 : ```plumed
116 : # Calculate the Q1 parameters for all the atoms
117 : q1: Q1 SPECIES=1-400 SWITCH={RATIONAL D_0=1.3 R_0=0.2 D_MAX=3.0}
118 : # Fixed virtual atom which serves as the probe volume's center (pos. in nm)
119 : center: FIXEDATOM AT=2.5,2.5,2.5
120 : # Vector in which element i is one if atom i is in sphere of interest and zero otherwise
121 : sphere: INSPHERE ATOMS=1-400 CENTER=center RADIUS={GAUSSIAN D_0=0.5 R_0=0.01 D_MAX=0.52}
122 : # Calculate the local_q1 parameters
123 : lq1: LOCAL_Q1 SPECIES=q1 SWITCH={RATIONAL D_0=1.3 R_0=0.2 D_MAX=3.0} MASK=sphere
124 : # Multiply fccubic parameters numbers by sphere vector
125 : prod: CUSTOM ARG=lq1,sphere FUNC=x*y PERIODIC=NO
126 : # Sum of coordination numbers for atoms that are in the sphere of interest
127 : numer: SUM ARG=prod PERIODIC=NO
128 : # Number of atoms that are in sphere of interest
129 : denom: SUM ARG=sphere PERIODIC=NO
130 : # Average coordination number for atoms in sphere of interest
131 : av: CUSTOM ARG=numer,denom FUNC=x/y PERIODIC=NO
132 : # And print out final CV to a file
133 : PRINT ARG=av FILE=colvar STRIDE=1
134 : ```
135 :
136 : This input calculates the average $s_i$ parameter for those atoms that are within a spherical region that is centered on the point $(2.5,2.5,2.5)$. By including the MASK
137 : keyword in the LOCAL_Q1 line we reduce the number of $s_i$ values we have to compute using the expression above. However, we are still asking PLUMED to calculate the $q_{lm}$ vectors for many atoms that
138 : will not contribute to the final averaged quantity. The documentation for [LOCAL_AVERAGE](LOCAL_AVERAGE.md) discusses how you can create a mask vector that can act upon
139 : the [Q1](Q1.md) action here that will ensure that you are not calculating $q_{lm}$ parameters that are not required.
140 :
141 : */
142 : //+ENDPLUMEDOC
143 :
144 : //+PLUMEDOC MCOLVARF LOCAL_Q3
145 : /*
146 : 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.
147 :
148 : 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
149 : 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
150 : 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
151 : 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
152 : 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
153 : 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
154 : because the number of atoms is relatively small.
155 :
156 : When the average [Q3](Q3.md) parameter is used to bias the dynamics a problems
157 : can occur. These averaged coordinates cannot distinguish between the correct,
158 : single-nucleus pathway and a concerted pathway in which all the atoms rearrange
159 : themselves into their solid-like configuration simultaneously. This second type
160 : of pathway would be impossible in reality because there is a large entropic
161 : barrier that prevents concerted processes like this from happening. However,
162 : in the finite sized systems that are commonly simulated this barrier is reduced
163 : substantially. As a result in simulations where average Steinhardt parameters
164 : are biased there are often quite dramatic system size effects
165 :
166 : If one wants to simulate nucleation using some form on
167 : biased dynamics what is really required is an order parameter that measures:
168 :
169 : - Whether or not the coordination spheres around atoms are ordered
170 : - Whether or not the atoms that are ordered are clustered together in a crystalline nucleus
171 :
172 : 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
173 : 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
174 : Steinhardt parameters. However, the inputs that you would need to write to implement common symmetry functions are rather complex so we also provide shortcuts
175 : like this one to help you compute CVs that have been widely used in the literature easily.
176 :
177 : This particular shortcut allows you to compute the LOCAL_Q3 parameter for a particular atom, which is a number that measures the extent to
178 : 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
179 : quantity for each of the atoms in the system:
180 :
181 : $$
182 : s_i = \frac{ \sum_j \sigma( r_{ij} ) \sum_{m=-3}^3 q_{3m}^{*}(i)q_{3m}(j) }{ \sum_j \sigma( r_{ij} ) }
183 : $$
184 :
185 : 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
186 : conjugation. The function
187 : $\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
188 : 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
189 : 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
190 : adjacent atoms are correlated.
191 :
192 : 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
193 : quantity to a file called colvar.
194 :
195 : ```plumed
196 : q3: Q3 SPECIES=1-64 D_0=1.3 R_0=0.2
197 : lq3: LOCAL_Q3 SPECIES=q3 SWITCH={RATIONAL D_0=1.3 R_0=0.2}
198 : lq3_mean: MEAN ARG=lq3 PERIODIC=NO
199 : PRINT ARG=lq3.mean FILE=colvar
200 : ```
201 :
202 : The following input calculates the distribution of LOCAL_Q3 parameters at any given time and outputs this information to a file.
203 :
204 : ```plumed
205 : q3: Q3 SPECIES=1-64 D_0=1.3 R_0=0.2
206 : 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}
207 : PRINT ARG=lq3.* FILE=colvar
208 : ```
209 :
210 : 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
211 : are done with those of all the other atoms in the system. The final quantity is the average and is outputted to a file
212 :
213 : ```plumed
214 : q3a: Q3 SPECIESA=1-5 SPECIESB=1-64 D_0=1.3 R_0=0.2
215 : q3b: Q3 SPECIESA=6-64 SPECIESB=1-64 D_0=1.3 R_0=0.2
216 : w3: LOCAL_Q3 SPECIESA=q3a SPECIESB=q3b SWITCH={RATIONAL D_0=1.3 R_0=0.2} MEAN
217 : PRINT ARG=w3.* FILE=colvar
218 : ```
219 :
220 : ## The MASK keyword
221 :
222 : 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
223 : 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
224 : from the [volumes module](module_volumes.md) to calculate the average value of the $s_i$ parameter that is defined by the equation above for only those atoms that
225 : lie in a certain part of the simulation box.
226 :
227 : ```plumed
228 : # Calculate the Q3 parameters for all the atoms
229 : q3: Q3 SPECIES=1-400 SWITCH={RATIONAL D_0=1.3 R_0=0.2 D_MAX=3.0}
230 : # Fixed virtual atom which serves as the probe volume's center (pos. in nm)
231 : center: FIXEDATOM AT=2.5,2.5,2.5
232 : # Vector in which element i is one if atom i is in sphere of interest and zero otherwise
233 : sphere: INSPHERE ATOMS=1-400 CENTER=center RADIUS={GAUSSIAN D_0=0.5 R_0=0.01 D_MAX=0.52}
234 : # Calculate the local_q3 parameters
235 : lq3: LOCAL_Q3 SPECIES=q3 SWITCH={RATIONAL D_0=1.3 R_0=0.2 D_MAX=3.0} MASK=sphere
236 : # Multiply fccubic parameters numbers by sphere vector
237 : prod: CUSTOM ARG=lq3,sphere FUNC=x*y PERIODIC=NO
238 : # Sum of coordination numbers for atoms that are in the sphere of interest
239 : numer: SUM ARG=prod PERIODIC=NO
240 : # Number of atoms that are in sphere of interest
241 : denom: SUM ARG=sphere PERIODIC=NO
242 : # Average coordination number for atoms in sphere of interest
243 : av: CUSTOM ARG=numer,denom FUNC=x/y PERIODIC=NO
244 : # And print out final CV to a file
245 : PRINT ARG=av FILE=colvar STRIDE=1
246 : ```
247 :
248 : This input calculates the average $s_i$ parameter for those atoms that are within a spherical region that is centered on the point $(2.5,2.5,2.5)$. By including the MASK
249 : keyword in the LOCAL_Q3 line we reduce the number of $s_i$ values we have to compute using the expression above. However, we are still asking PLUMED to calculate the $q_{lm}$ vectors for many atoms that
250 : will not contribute to the final averaged quantity. The documentation for [LOCAL_AVERAGE](LOCAL_AVERAGE.md) discusses how you can create a mask vector that can act upon
251 : the [Q3](Q3.md) action here that will ensure that you are not calculating $q_{lm}$ vectors that are not required.
252 :
253 : */
254 : //+ENDPLUMEDOC
255 :
256 : //+PLUMEDOC MCOLVARF LOCAL_Q4
257 : /*
258 : 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.
259 :
260 : 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
261 : 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
262 : 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
263 : 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
264 : 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
265 : 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
266 : because the number of atoms is relatively small.
267 :
268 : When the average [Q4](Q4.md) parameter is used to bias the dynamics a problems
269 : can occur. These averaged coordinates cannot distinguish between the correct,
270 : single-nucleus pathway and a concerted pathway in which all the atoms rearrange
271 : themselves into their solid-like configuration simultaneously. This second type
272 : of pathway would be impossible in reality because there is a large entropic
273 : barrier that prevents concerted processes like this from happening. However,
274 : in the finite sized systems that are commonly simulated this barrier is reduced
275 : substantially. As a result in simulations where average Steinhardt parameters
276 : are biased there are often quite dramatic system size effects
277 :
278 : If one wants to simulate nucleation using some form on
279 : biased dynamics what is really required is an order parameter that measures:
280 :
281 : - Whether or not the coordination spheres around atoms are ordered
282 : - Whether or not the atoms that are ordered are clustered together in a crystalline nucleus
283 :
284 : 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
285 : 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
286 : Steinhardt parameters. However, the inputs that you would need to write to implement common symmetry functions are rather complex so we also provide shortcuts
287 : like this one to help you compute CVs that have been widely used in the literature easily.
288 :
289 : This particular shortcut allows you to compute the LOCAL_Q4 parameter for a particular atom, which is a number that measures the extent to
290 : 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
291 : quantity for each of the atoms in the system:
292 :
293 : $$
294 : s_i = \frac{ \sum_j \sigma( r_{ij} ) \sum_{m=-4}^4 q_{4m}^{*}(i)q_{4m}(j) }{ \sum_j \sigma( r_{ij} ) }
295 : $$
296 :
297 : 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
298 : conjugation. The function
299 : $\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
300 : 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
301 : 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
302 : adjacent atoms are correlated.
303 :
304 : 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
305 : quantity to a file called colvar.
306 :
307 : ```plumed
308 : q4: Q4 SPECIES=1-64 D_0=1.3 R_0=0.2
309 : lq4: LOCAL_Q4 SPECIES=q4 SWITCH={RATIONAL D_0=1.3 R_0=0.2}
310 : lq4_mean: MEAN ARG=lq4 PERIODIC=NO
311 : PRINT ARG=lq4_mean FILE=colvar
312 : ```
313 :
314 : The following input calculates the distribution of LOCAL_Q4 parameters at any given time and outputs this information to a file.
315 :
316 : ```plumed
317 : q4: Q4 SPECIES=1-64 D_0=1.3 R_0=0.2
318 : 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}
319 : PRINT ARG=lq4.* FILE=colvar
320 : ```
321 :
322 : 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
323 : are done with those of all the other atoms in the system. The final quantity is the average and is outputted to a file
324 :
325 : ```plumed
326 : q4a: Q4 SPECIESA=1-5 SPECIESB=1-64 D_0=1.3 R_0=0.2
327 : q4b: Q4 SPECIESA=6-64 SPECIESB=1-64 D_0=1.3 R_0=0.2
328 : w4: LOCAL_Q4 SPECIESA=q4a SPECIESB=q4b SWITCH={RATIONAL D_0=1.3 R_0=0.2} MEAN
329 : PRINT ARG=w4.* FILE=colvar
330 : ```
331 :
332 : ## The MASK keyword
333 :
334 : 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
335 : 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
336 : from the [volumes module](module_volumes.md) to calculate the average value of the $s_i$ parameter that is defined by the equation above for only those atoms that
337 : lie in a certain part of the simulation box.
338 :
339 : ```plumed
340 : # Calculate the Q4 parameters for all the atoms
341 : q4: Q4 SPECIES=1-400 SWITCH={RATIONAL D_0=1.3 R_0=0.2 D_MAX=3.0}
342 : # Fixed virtual atom which serves as the probe volume's center (pos. in nm)
343 : center: FIXEDATOM AT=2.5,2.5,2.5
344 : # Vector in which element i is one if atom i is in sphere of interest and zero otherwise
345 : sphere: INSPHERE ATOMS=1-400 CENTER=center RADIUS={GAUSSIAN D_0=0.5 R_0=0.01 D_MAX=0.52}
346 : # Calculate the local_q4 parameters
347 : lq4: LOCAL_Q4 SPECIES=q4 SWITCH={RATIONAL D_0=1.3 R_0=0.2 D_MAX=3.0} MASK=sphere
348 : # Multiply fccubic parameters numbers by sphere vector
349 : prod: CUSTOM ARG=lq4,sphere FUNC=x*y PERIODIC=NO
350 : # Sum of coordination numbers for atoms that are in the sphere of interest
351 : numer: SUM ARG=prod PERIODIC=NO
352 : # Number of atoms that are in sphere of interest
353 : denom: SUM ARG=sphere PERIODIC=NO
354 : # Average coordination number for atoms in sphere of interest
355 : av: CUSTOM ARG=numer,denom FUNC=x/y PERIODIC=NO
356 : # And print out final CV to a file
357 : PRINT ARG=av FILE=colvar STRIDE=1
358 : ```
359 :
360 : This input calculates the average $s_i$ parameter for those atoms that are within a spherical region that is centered on the point $(2.5,2.5,2.5)$. By including the MASK
361 : keyword in the LOCAL_Q4 line we reduce the number of $s_i$ values we have to compute using the expression above. However, we are still asking PLUMED to calculate the $q_{lm}$ vectors for many atoms that
362 : will not contribute to the final averaged quantity. The documentation for [LOCAL_AVERAGE](LOCAL_AVERAGE.md) discusses how you can create a mask vector that can act upon
363 : the [Q4](Q4.md) action here that will ensure that you are not calculating $q_{lm}$ vectors that are not required.
364 :
365 : */
366 : //+ENDPLUMEDOC
367 :
368 : //+PLUMEDOC MCOLVARF LOCAL_Q6
369 : /*
370 : 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.
371 :
372 : 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
373 : 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
374 : 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
375 : 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
376 : 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
377 : 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
378 : because the number of atoms is relatively small.
379 :
380 : When the average [Q6](Q6.md) parameter is used to bias the dynamics a problems
381 : can occur. These averaged coordinates cannot distinguish between the correct,
382 : single-nucleus pathway and a concerted pathway in which all the atoms rearrange
383 : themselves into their solid-like configuration simultaneously. This second type
384 : of pathway would be impossible in reality because there is a large entropic
385 : barrier that prevents concerted processes like this from happening. However,
386 : in the finite sized systems that are commonly simulated this barrier is reduced
387 : substantially. As a result in simulations where average Steinhardt parameters
388 : are biased there are often quite dramatic system size effects
389 :
390 : If one wants to simulate nucleation using some form on
391 : biased dynamics what is really required is an order parameter that measures:
392 :
393 : - Whether or not the coordination spheres around atoms are ordered
394 : - Whether or not the atoms that are ordered are clustered together in a crystalline nucleus
395 :
396 : 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
397 : 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
398 : Steinhardt parameters. However, the inputs that you would need to write to implement common symmetry functions are rather complex so we also provide shortcuts
399 : like this one to help you compute CVs that have been widely used in the literature easily.
400 :
401 : This particular shortcut allows you to compute the LOCAL_Q6 parameter for a particular atom, which is a number that measures the extent to
402 : 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
403 : quantity for each of the atoms in the system:
404 :
405 : $$
406 : s_i = \frac{ \sum_j \sigma( r_{ij} ) \sum_{m=-6}^6 q_{6m}^{*}(i)q_{6m}(j) }{ \sum_j \sigma( r_{ij} ) }
407 : $$
408 :
409 : 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
410 : conjugation. The function
411 : $\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
412 : 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
413 : 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
414 : adjacent atoms are correlated.
415 :
416 : 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
417 : quantity to a file called colvar.
418 :
419 : ```plumed
420 : q6: Q6 SPECIES=1-64 D_0=1.3 R_0=0.2
421 : lq6: LOCAL_Q6 SPECIES=q6 SWITCH={RATIONAL D_0=1.3 R_0=0.2}
422 : lq6_mean: MEAN ARG=lq6 PERIODIC=NO
423 : PRINT ARG=lq6_mean FILE=colvar
424 : ```
425 :
426 : The following input calculates the distribution of LOCAL_Q6 parameters at any given time and outputs this information to a file.
427 :
428 : ```plumed
429 : q6: Q6 SPECIES=1-64 D_0=1.3 R_0=0.2
430 : 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}
431 : PRINT ARG=lq6.* FILE=colvar
432 : ```
433 :
434 : 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
435 : are done with those of all the other atoms in the system. The final quantity is the average and is outputted to a file
436 :
437 : ```plumed
438 : q6a: Q6 SPECIESA=1-5 SPECIESB=1-64 D_0=1.3 R_0=0.2
439 : q6b: Q6 SPECIESA=6-64 SPECIESB=1-64 D_0=1.3 R_0=0.2
440 : w6: LOCAL_Q6 SPECIESA=q6a SPECIESB=q6b SWITCH={RATIONAL D_0=1.3 R_0=0.2} MEAN
441 : PRINT ARG=w6.* FILE=colvar
442 : ```
443 :
444 : ## The MASK keyword
445 :
446 : 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
447 : 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
448 : from the [volumes module](module_volumes.md) to calculate the average value of the $s_i$ parameter that is defined by the equation above for only those atoms that
449 : lie in a certain part of the simulation box.
450 :
451 : ```plumed
452 : # Calculate the Q6 parameters for all the atoms
453 : q6: Q6 SPECIES=1-400 SWITCH={RATIONAL D_0=1.3 R_0=0.2 D_MAX=3.0}
454 : # Fixed virtual atom which serves as the probe volume's center (pos. in nm)
455 : center: FIXEDATOM AT=2.5,2.5,2.5
456 : # Vector in which element i is one if atom i is in sphere of interest and zero otherwise
457 : sphere: INSPHERE ATOMS=1-400 CENTER=center RADIUS={GAUSSIAN D_0=0.5 R_0=0.01 D_MAX=0.52}
458 : # Calculate the local_q6 parameters
459 : lq6: LOCAL_Q6 SPECIES=q6 SWITCH={RATIONAL D_0=1.3 R_0=0.2 D_MAX=3.0} MASK=sphere
460 : # Multiply fccubic parameters numbers by sphere vector
461 : prod: CUSTOM ARG=lq6,sphere FUNC=x*y PERIODIC=NO
462 : # Sum of coordination numbers for atoms that are in the sphere of interest
463 : numer: SUM ARG=prod PERIODIC=NO
464 : # Number of atoms that are in sphere of interest
465 : denom: SUM ARG=sphere PERIODIC=NO
466 : # Average coordination number for atoms in sphere of interest
467 : av: CUSTOM ARG=numer,denom FUNC=x/y PERIODIC=NO
468 : # And print out final CV to a file
469 : PRINT ARG=av FILE=colvar STRIDE=1
470 : ```
471 :
472 : This input calculates the average $s_i$ parameter for those atoms that are within a spherical region that is centered on the point $(2.5,2.5,2.5)$. By including the MASK
473 : keyword in the LOCAL_Q6 line we reduce the number of $s_i$ values we have to compute using the expression above. However, we are still asking PLUMED to calculate the $q_{lm}$ vectors for many atoms that
474 : will not contribute to the final averaged quantity. The documentation for [LOCAL_AVERAGE](LOCAL_AVERAGE.md) discusses how you can create a mask vector that can act upon
475 : the [Q6](Q6.md) action here that will ensure that you are not calculating $q_{lm}$ vectors that are not required.
476 :
477 : */
478 : //+ENDPLUMEDOC
479 :
480 : class LocalSteinhardt : public ActionShortcut {
481 : private:
482 : std::string getSymbol( const int& m ) const ;
483 : std::string getArgsForStack( const int& l, const std::string& lab ) const;
484 : public:
485 : static void registerKeywords( Keywords& keys );
486 : explicit LocalSteinhardt(const ActionOptions&);
487 : };
488 :
489 : PLUMED_REGISTER_ACTION(LocalSteinhardt,"LOCAL_Q1")
490 : PLUMED_REGISTER_ACTION(LocalSteinhardt,"LOCAL_Q3")
491 : PLUMED_REGISTER_ACTION(LocalSteinhardt,"LOCAL_Q4")
492 : PLUMED_REGISTER_ACTION(LocalSteinhardt,"LOCAL_Q6")
493 :
494 20 : void LocalSteinhardt::registerKeywords( Keywords& keys ) {
495 20 : ActionShortcut::registerKeywords( keys );
496 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");
497 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");
498 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");
499 20 : keys.add("optional","SWITCH","This keyword is used if you want to employ an alternative to the continuous swiching function defined above. "
500 : "The following provides information on the \\ref switchingfunction that are available. "
501 : "When this keyword is present you no longer need the NN, MM, D_0 and R_0 keywords.");
502 20 : keys.add("optional","MASK","the label/s for vectors that are used to determine which local steinhardt parameters to compute");
503 40 : keys.addDeprecatedFlag("LOWMEM","");
504 40 : keys.setValueDescription("vector","the values of the local steinhardt parameters for the input atoms");
505 20 : multicolvar::MultiColvarShortcuts::shortcutKeywords( keys );
506 20 : keys.needsAction("CONTACT_MATRIX");
507 20 : keys.needsAction("MATRIX_PRODUCT");
508 20 : keys.needsAction("GROUP");
509 20 : keys.needsAction("ONES");
510 20 : keys.needsAction("OUTER_PRODUCT");
511 20 : keys.needsAction("VSTACK");
512 20 : keys.needsAction("CONCATENATE");
513 20 : keys.needsAction("CUSTOM");
514 20 : keys.needsAction("TRANSPOSE");
515 40 : keys.needsAction("MATRIX_VECTOR_PRODUCT");
516 20 : keys.addFlag("USEGPU",false,"run part of this calculation on the GPU");
517 20 : }
518 :
519 98 : std::string LocalSteinhardt::getSymbol( const int& m ) const {
520 98 : if( m<0 ) {
521 : std::string num;
522 40 : Tools::convert( -1*m, num );
523 40 : return "n" + num;
524 58 : } else if( m>0 ) {
525 : std::string num;
526 49 : Tools::convert( m, num );
527 49 : return "p" + num;
528 : }
529 9 : return "0";
530 : }
531 :
532 9 : std::string LocalSteinhardt::getArgsForStack( const int& l, const std::string& sp_lab ) const {
533 : std::string numstr;
534 9 : Tools::convert( l, numstr );
535 18 : std::string data_mat = " ARG=" + sp_lab + "_sp.rm-n" + numstr + ","
536 18 : + sp_lab + "_sp.im-n" + numstr;
537 107 : for(int i=-l+1; i<=l; ++i) {
538 98 : numstr = getSymbol( i );
539 196 : data_mat += "," + sp_lab + "_sp.rm-" + numstr + ","
540 196 : + sp_lab + "_sp.im-" + numstr;
541 : }
542 9 : return data_mat;
543 : }
544 :
545 7 : LocalSteinhardt::LocalSteinhardt(const ActionOptions& ao):
546 : Action(ao),
547 7 : ActionShortcut(ao) {
548 :
549 : bool usegpu;
550 7 : parseFlag("USEGPU",usegpu);
551 14 : const std::string doUSEGPU = usegpu?" USEGPU":"";
552 :
553 : #define createLabel(name) const std::string name##Lab = getShortcutLabel()+"_"#name;
554 : bool lowmem;
555 7 : parseFlag("LOWMEM",lowmem);
556 7 : if( lowmem ) {
557 0 : warning("LOWMEM flag is deprecated and is no longer required for this action");
558 : }
559 : // Get the Q value
560 : int l;
561 14 : Tools::convert( getName().substr(7), l);
562 : // Create a vector filled with ones
563 : std::string twolplusone;
564 7 : Tools::convert( 2*(2*l+1), twolplusone );
565 7 : createLabel( uvec );
566 14 : readInputLine( uvecLab + ": ONES SIZE=" + twolplusone );
567 : // Read in species keyword
568 : std::string sp_str;
569 14 : parse("SPECIES",sp_str);
570 : std::string spa_str;
571 7 : parse("SPECIESA",spa_str);
572 7 : createLabel( dpmat );
573 7 : createLabel( cmap );
574 7 : createLabel( grp );
575 7 : if( sp_str.length()>0 ) {
576 : // Create a group with these atoms
577 12 : readInputLine( grpLab + ": GROUP ATOMS=" + sp_str );
578 6 : std::vector<std::string> sp_lab = Tools::getWords(sp_str, "\t\n ,");
579 : // This creates the stash to hold all the vectors
580 6 : createLabel( uvecs );
581 6 : createLabel( nmat );
582 6 : if( sp_lab.size()==1 ) {
583 : // The lengths of all the vectors in a vector
584 18 : readInputLine( nmatLab + ": OUTER_PRODUCT "
585 12 : "ARG=" + sp_lab[0] + "_norm," + uvecLab);
586 : // The unormalised vectors
587 12 : readInputLine( uvecsLab + ": VSTACK" + getArgsForStack( l, sp_lab[0] ) );
588 : } else {
589 0 : createLabel( mags );
590 0 : std::string len_vec = magsLab + ": CONCATENATE ARG=" + sp_lab[0] + "_norm";
591 0 : for(unsigned i=1; i<sp_lab.size(); ++i) {
592 0 : len_vec += "," + sp_lab[i] + "_norm";
593 : }
594 : // This is the vector that contains all the magnitudes
595 0 : readInputLine( len_vec );
596 0 : std::string concat_str = uvecsLab + ": CONCATENATE";
597 0 : for(unsigned i=0; i<sp_lab.size(); ++i) {
598 : std::string snum;
599 0 : Tools::convert( i+1, snum );
600 0 : concat_str += " MATRIX" + snum + "1=" + uvecsLab + snum;
601 0 : readInputLine( uvecsLab + snum + ": VSTACK" + getArgsForStack( l, sp_lab[i] ) );
602 : }
603 : // And the normalising matrix by taking the column vector of magnitudes and multiplying by the row vector of ones
604 0 : readInputLine( nmatLab + ": OUTER_PRODUCT ARG=" + magsLab + "," + uvecLab);
605 : // The unormalised vectors
606 0 : readInputLine( concat_str );
607 : }
608 : // Now normalise all the vectors by doing Hadammard "product" with normalising matrix
609 6 : createLabel( vecs );
610 18 : readInputLine( vecsLab + ": CUSTOM ARG=" + uvecsLab + ","
611 12 : + nmatLab + " FUNC=x/y PERIODIC=NO");
612 : // And transpose the matrix
613 6 : createLabel( vecsT );
614 12 : readInputLine( vecsTLab + ": TRANSPOSE ARG=" + vecsLab );
615 : std::string sw_str;
616 12 : parse("SWITCH",sw_str);
617 : std::string maskstr;
618 12 : parse("MASK",maskstr);
619 6 : if( maskstr.length()>0 ) {
620 4 : maskstr=" MASK=" + maskstr;
621 : }
622 18 : readInputLine( cmapLab + ": CONTACT_MATRIX GROUP=" + sp_str + " "
623 12 : "SWITCH={" + sw_str + "}" + maskstr + doUSEGPU);
624 : // And the matrix of dot products
625 18 : readInputLine( dpmatLab + ": MATRIX_PRODUCT ARG=" + vecsLab + ","
626 12 : + vecsTLab + " MASK=" + cmapLab + doUSEGPU);
627 7 : } else if( spa_str.length()>0 ) {
628 : // Create a group with these atoms
629 2 : readInputLine( grpLab + ": GROUP ATOMS=" + spa_str );
630 : std::string spb_str;
631 2 : parse("SPECIESB",spb_str);
632 1 : if( spb_str.length()==0 ) {
633 0 : plumed_merror("need both SPECIESA and SPECIESB in input");
634 : }
635 1 : const std::vector<std::string> sp_laba = Tools::getWords(spa_str, "\t\n ,");
636 1 : const std::vector<std::string> sp_labb = Tools::getWords(spb_str, "\t\n ,");
637 1 : createLabel(nmatA);
638 1 : createLabel(uvecsA);
639 1 : if( sp_laba.size()==1 ) {
640 : // The matrix that is used for normalising
641 3 : readInputLine( nmatALab + ": OUTER_PRODUCT "
642 2 : "ARG=" + sp_laba[0] + "_norm," + uvecLab);
643 : // The unormalised vectors
644 2 : readInputLine( uvecsALab + ": VSTACK" + getArgsForStack( l, sp_laba[0] ) );
645 : } else {
646 0 : createLabel(magsA);
647 0 : std::string len_vec = magsALab + ": CONCATENATE "
648 0 : "ARG=" + sp_laba[0] + "_norm";
649 0 : for(unsigned i=1; i<sp_laba.size(); ++i) {
650 0 : len_vec += "," + sp_laba[i] + "_norm";
651 : }
652 : // This is the vector that contains all the magnitudes
653 0 : readInputLine( len_vec );
654 0 : std::string concat_str = uvecsALab + ": CONCATENATE";
655 0 : for(unsigned i=0; i<sp_laba.size(); ++i) {
656 : std::string snum;
657 0 : Tools::convert( i+1, snum );
658 0 : concat_str += " MATRIX" + snum + "1=" + uvecsALab + snum;
659 0 : readInputLine( uvecsALab + snum + ": VSTACK" + getArgsForStack( l, sp_laba[i] ) );
660 : }
661 : // And the normalising matrix by taking the column vector of magnitudes and multiplying by the row vector of ones
662 0 : readInputLine( nmatALab + ": OUTER_PRODUCT ARG=" + magsALab + "," + uvecLab);
663 : // The unormalised vector
664 0 : readInputLine( concat_str );
665 : }
666 : // Now normalise all the vectors by doing Hadammard "product" with normalising matrix
667 1 : createLabel( vecsA );
668 3 : readInputLine( vecsALab + ": CUSTOM ARG=" + uvecsALab + ","
669 2 : + nmatALab + " FUNC=x/y PERIODIC=NO");
670 : // Now do second matrix
671 1 : createLabel(nmatB);
672 1 : createLabel(uvecsBT);
673 1 : createLabel(uvecsB);
674 1 : if( sp_labb.size()==1 ) {
675 0 : readInputLine( nmatBLab + ": OUTER_PRODUCT ARG=" + uvecLab + ","
676 0 : + sp_labb[0] + "_norm");
677 0 : readInputLine( uvecsBTLab + ": VSTACK" + getArgsForStack( l, sp_labb[0] ) );
678 0 : readInputLine( uvecsBLab + ": TRANSPOSE ARG=" + uvecsBTLab);
679 : } else {
680 1 : createLabel(magsB);
681 2 : std::string len_vec = magsBLab + ": CONCATENATE ARG=" + sp_labb[0] + "_norm";
682 2 : for(unsigned i=1; i<sp_labb.size(); ++i) {
683 2 : len_vec += "," + sp_labb[i] + "_norm";
684 : }
685 : // This is the vector that contains all the magnitudes
686 1 : readInputLine( len_vec );
687 1 : std::string concat_str = uvecsBLab + ": CONCATENATE";
688 3 : for(unsigned i=0; i<sp_labb.size(); ++i) {
689 : std::string snum;
690 2 : Tools::convert( i+1, snum );
691 4 : concat_str += " MATRIX1" + snum + "=" + uvecsBLab + snum;
692 4 : readInputLine( uvecsBTLab + snum + ": VSTACK" + getArgsForStack( l, sp_labb[i] ) );
693 4 : readInputLine( uvecsBLab + snum + ": TRANSPOSE ARG=" + uvecsBTLab + snum );
694 : }
695 : // And the normalising matrix
696 2 : readInputLine( nmatBLab + ": OUTER_PRODUCT ARG=" + uvecLab + "," + magsBLab);
697 : // The unormalised vectors
698 1 : readInputLine( concat_str );
699 : }
700 : // Now normalise all the vectors by doing Hadammard "product" with normalising matrix
701 1 : createLabel(vecsB);
702 3 : readInputLine( vecsBLab + ": CUSTOM ARG=" + uvecsBLab + ","
703 2 : + nmatBLab + " FUNC=x/y PERIODIC=NO");
704 : std::string sw_str;
705 2 : parse("SWITCH",sw_str);
706 : std::string maskstr;
707 2 : parse("MASK",maskstr);
708 1 : if( maskstr.length()>0 ) {
709 0 : maskstr=" MASK=" + maskstr;
710 : }
711 3 : readInputLine( cmapLab + ": CONTACT_MATRIX GROUPA=" + spa_str
712 2 : + " GROUPB=" + spb_str + " SWITCH={" + sw_str + "}" + maskstr + doUSEGPU);
713 3 : readInputLine( dpmatLab + ": MATRIX_PRODUCT ARG=" + vecsALab + "," + vecsBLab
714 2 : + " MASK=" + cmapLab + doUSEGPU);
715 1 : }
716 :
717 : // Now create the product matrix
718 7 : createLabel( prod );
719 14 : readInputLine( prodLab + ": CUSTOM ARG=" + cmapLab + "," + dpmatLab + " FUNC=x*y PERIODIC=NO");
720 : // Now the sum of coordination numbers times the switching functions
721 7 : ActionWithValue* av = plumed.getActionSet().selectWithLabel<ActionWithValue*>( cmapLab);
722 7 : plumed_assert( av && av->getNumberOfComponents()>0 && (av->copyOutput(0))->getRank()==2 );
723 : std::string size;
724 7 : Tools::convert( (av->copyOutput(0))->getShape()[1], size );
725 7 : createLabel( ones );
726 14 : readInputLine( onesLab + ": ONES SIZE=" + size );
727 14 : readInputLine( getShortcutLabel() + ": MATRIX_VECTOR_PRODUCT "
728 21 : "ARG=" + prodLab +"," + getShortcutLabel() +"_ones" + doUSEGPU);
729 : // And just the sum of the coordination numbers
730 14 : readInputLine( getShortcutLabel() + "_denom: MATRIX_VECTOR_PRODUCT "
731 21 : "ARG=" + cmapLab + "," + getShortcutLabel() +"_ones" + doUSEGPU);
732 : // And matheval to get the final quantity
733 7 : createLabel( av );
734 21 : readInputLine( avLab + ": CUSTOM ARG=" + getShortcutLabel() + ","
735 21 : + getShortcutLabel() + "_denom FUNC=x/y PERIODIC=NO");
736 : // And this expands everything
737 14 : multicolvar::MultiColvarShortcuts::expandFunctions( getShortcutLabel(), avLab, "", this );
738 : #undef createLabel
739 7 : }
740 :
741 : }
742 : }
|
