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 "CoordinationNumbers.h"
23 : #include "core/ActionShortcut.h"
24 : #include "multicolvar/MultiColvarShortcuts.h"
25 : #include "core/ActionWithValue.h"
26 : #include "core/PlumedMain.h"
27 : #include "core/ActionSet.h"
28 : #include "core/ActionRegister.h"
29 :
30 : #include <complex>
31 :
32 : namespace PLMD {
33 : namespace symfunc {
34 :
35 : //+PLUMEDOC MCOLVAR Q1
36 : /*
37 : Calculate 1st order Steinhardt parameters
38 :
39 : The 1st order Steinhardt parameters allow us to measure the degree to which the first coordination shell
40 : around an atom is ordered with the atoms aranged on a line. The Steinhardt parameter for atom, $i$ is complex vector whose components are
41 : calculated using the following formula:
42 :
43 : $$
44 : q_{1m}(i) = \sum_j \sigma( r_{ij} ) Y_{1m}(\mathbf{r}_{ij})
45 : $$
46 :
47 : where $Y_{1m}$ is one of the 1st order spherical harmonics so $m$ is a number that runs from $-1$ to
48 : $+1$. The function $\sigma( r_{ij} )$ is a switching function that acts on the distance between
49 : atoms $i$ and $j$. The parameters of this function should be set so that it the function is equal to one
50 : when atom $j$ is in the first coordination sphere of atom $i$ and is zero otherwise.
51 :
52 : As discussed on [this page](https://www.plumed-tutorials.org/lessons/23/001/data/Steinhardt.html), the Steinhardt parameters can
53 : be used to measure the degree of order in the system in a variety of different ways. The
54 : simplest way of measuring whether or not the coordination sphere is ordered is to simply take the norm of the above vector i.e.
55 :
56 : $$
57 : Q_1(i) = \frac{1}{\sum_j \sigma(r_{ij}) } \sqrt{ \sum_{m=-1}^1 q_{1m}(i)^{*} q_{1m}(i) }
58 : $$
59 :
60 : This norm is small when the coordination shell is disordered and larger when the coordination shell is ordered. The following input illustrates
61 : how by averaging the value of this norm over all the atoms in the system you can measure the global degree of order in the system:
62 :
63 : ```plumed
64 : q1: Q1 SPECIES=1-64 D_0=1.3 R_0=0.2
65 : q1_mean: MEAN ARG=q1 PERIODIC=NO
66 : PRINT ARG=q1_mean FILE=colvar
67 : ```
68 :
69 : In the above input the rational [switching function](LESS_THAN.md) with the parameters above. We would recommend using SWITCH syntax
70 : rather than the syntax above when giving the parameters for the switching function as you can then use any of the switching functions described
71 : 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
72 : switching function as demonstrated below:
73 :
74 : ```plumed
75 : q1: Q1 SPECIES=1-64 SWITCH={RATIONAL D_0=1.3 R_0=0.2 D_MAX=3.0}
76 : q1_mean: MEAN ARG=q1 PERIODIC=NO
77 : PRINT ARG=q1_mean FILE=colvar
78 : ```
79 :
80 : 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
81 : of the documentation for [CONTACT_MATRIX](CONTACT_MATRIX.md).
82 :
83 : ## Working with two types of atoms
84 :
85 : The command below could be used to measure the Q1 parameters that describe the arrangement of chlorine ions around the
86 : sodium atoms in sodium chloride. The imagined system here is composed of 64 NaCl formula units and the atoms are arranged in the input
87 : with the 64 Na$^+$ ions followed by the 64 Cl$-$ ions. Once again the average Q1 parameter is calculated and output to a
88 : file called colvar
89 :
90 : ```plumed
91 : q1: Q1 SPECIESA=1-64 SPECIESB=65-128 D_0=1.3 R_0=0.2
92 : q1_mean: MEAN ARG=q1 PERIODIC=NO
93 : PRINT ARG=q1_mean FILE=colvar
94 : ```
95 :
96 : If you simply want to examine the values of the Q1 parameters for each of the atoms in your system you can do so by exploiting the
97 : command [DUMPATOMS](DUMPATOMS.md) as shown in the example below. The following output file will output a file in an extended xyz format
98 : called q1.xyz for each frame of the analyzed MD trajectory. The first column in this file will contain a dummy name for each of the
99 : atoms, columns 2-4 will then contain the x, y and z positions of the atoms, column 5 will contain the value of the Q1 parameter, columns
100 : 6-8 will contain the real parts of the director of the $q_{1m}$ vector while columns 9-11 will contain the imaginary parts of this director.
101 :
102 : ```plumed
103 : q1: Q1 SPECIESA=1-64 SPECIESB=65-128 D_0=1.3 R_0=0.2 MEAN
104 : DUMPATOMS ATOMS=q1 ARG=q1 FILE=q1.xyz
105 : ```
106 :
107 : ## The MASK keyword
108 :
109 : 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
110 : input, which tells Q1 that it is safe not to calculate the Q1 parameter for some of the atoms. As illustrated below, this is useful if you are using functionality
111 : from the [volumes module](module_volumes.md) to calculate the average value of the Q1 parameter for only those atoms that lie in a certain part of the simulation box.
112 :
113 : ```plumed
114 : # Fixed virtual atom which serves as the probe volume's center (pos. in nm)
115 : center: FIXEDATOM AT=2.5,2.5,2.5
116 : # Vector in which element i is one if atom i is in sphere of interest and zero otherwise
117 : sphere: INSPHERE ATOMS=1-400 CENTER=center RADIUS={GAUSSIAN D_0=0.5 R_0=0.01 D_MAX=0.52}
118 : # Calculate the fccubic parameter of the atoms
119 : cc: Q1 ...
120 : SPECIES=1-400 MASK=sphere
121 : SWITCH={RATIONAL D_0=3.0 R_0=1.5 D_MAX=6.0}
122 : ...
123 : # Multiply fccubic parameters numbers by sphere vector
124 : prod: CUSTOM ARG=cc,sphere FUNC=x*y PERIODIC=NO
125 : # Sum of coordination numbers for atoms that are in the sphere of interest
126 : numer: SUM ARG=prod PERIODIC=NO
127 : # Number of atoms that are in sphere of interest
128 : denom: SUM ARG=sphere PERIODIC=NO
129 : # Average coordination number for atoms in sphere of interest
130 : av: CUSTOM ARG=numer,denom FUNC=x/y PERIODIC=NO
131 : # And print out final CV to a file
132 : PRINT ARG=av FILE=colvar STRIDE=1
133 : ```
134 :
135 : This input calculate the average value of the Q1 parameter for only those atoms that are within a spherical region that is centered on the point
136 : $(2.5,2.5,2.5)$.
137 :
138 : ## Deprecated syntax
139 :
140 : More information on the deprecated keywords that are given below is available in the documentation for the [DISTANCES](DISTANCES.md) command.
141 : Below is an example where these deprecated keywords are used to calculate the histogram of Q1 parameters for the 64 atoms in a box of Lennard Jones print them
142 : to a file called colvar:
143 :
144 : ```plumed
145 : q1: Q1 SPECIES=1-64 D_0=1.3 R_0=0.2 HISTOGRAM={GAUSSIAN LOWER=0.0 UPPER=1.0 NBINS=20 SMEAR=0.1}
146 : PRINT ARG=q1.* FILE=colvar
147 : ```
148 :
149 : The following example illustrates how you can use VSUM to calculate a global vector of $Q_{1m}$ values as follows:
150 :
151 : $$
152 : Q_{1m} = \sum_i \frac{q_{1m}(i)}{\sum_j \sigma(r_{ij})}
153 : $$
154 :
155 : where the sum runs over all the atoms. You can then take these $Q_{1m}$ values and compute the following norm:
156 :
157 : $$
158 : s = \sqrt{ \sum_{m=-1}^1 Q_{1m}^{*} Q_{1m} }
159 : $$
160 :
161 : The VMEAN command that is also used in the input below performs a similar operations. The only difference is that
162 : we divide the sums in the first expression above by the number of atoms.
163 :
164 : ```plumed
165 : q1: Q1 SPECIES=1-64 D_0=1.3 R_0=0.2 VMEAN VSUM
166 : PRINT ARG=q1.* FILE=colvar
167 : ```
168 :
169 : */
170 : //+ENDPLUMEDOC
171 :
172 : //+PLUMEDOC MCOLVAR Q3
173 : /*
174 : Calculate 3rd order Steinhardt parameters.
175 :
176 : The 3rd order Steinhardt parameters allow us to measure the degree to which the first coordination shell
177 : around an atom is ordered. The Steinhardt parameter for atom, $i$ is complex vector whose components are
178 : calculated using the following formula:
179 :
180 : $$
181 : q_{3m}(i) = \sum_j \sigma( r_{ij} ) Y_{3m}(\mathbf{r}_{ij})
182 : $$
183 :
184 : where $Y_{3m}$ is one of the 3rd order spherical harmonics so $m$ is a number that runs from $-3$ to
185 : $+3$. The function $\sigma( r_{ij} )$ is a switching function that acts on the distance between
186 : atoms $i$ and $j$. The parameters of this function should be set so that it the function is equal to one
187 : when atom $j$ is in the first coordination sphere of atom $i$ and is zero otherwise.
188 :
189 : As discussed on [this page](https://www.plumed-tutorials.org/lessons/23/001/data/Steinhardt.html), the Steinhardt parameters can
190 : be used to measure the degree of order in the system in a variety of different ways. The
191 : simplest way of measuring whether or not the coordination sphere is ordered is to simply take the norm of the above vector i.e.
192 :
193 : $$
194 : Q_3(i) = \frac{1}{\sum_j \sigma(r_{ij}) } \sqrt{ \sum_{m=-3}^3 q_{3m}(i)^{*} q_{3m}(i) }
195 : $$
196 :
197 : This norm is small when the coordination shell is disordered and larger when the coordination shell is ordered. The following input illustrates
198 : how by averaging the value of this norm over all the atoms in the system you can measure the global degree of order in the system:
199 :
200 : ```plumed
201 : q3: Q3 SPECIES=1-64 D_0=1.3 R_0=0.2
202 : q3_mean: MEAN ARG=q3 PERIODIC=NO
203 : PRINT ARG=q3_mean FILE=colvar
204 : ```
205 :
206 : In the above input the rational [switching function](LESS_THAN.md) with the parameters above. We would recommend using SWITCH syntax
207 : rather than the syntax above when giving the parameters for the switching function as you can then use any of the switching functions described
208 : 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
209 : switching function as demonstrated below:
210 :
211 : ```plumed
212 : q3: Q3 SPECIES=1-64 SWITCH={RATIONAL D_0=1.3 R_0=0.2 D_MAX=3.0}
213 : q3_mean: MEAN ARG=q3 PERIODIC=NO
214 : PRINT ARG=q3_mean FILE=colvar
215 : ```
216 :
217 : 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
218 : of the documentation for [CONTACT_MATRIX](CONTACT_MATRIX.md).
219 :
220 : ## Working with two types of atoms
221 :
222 : The command below could be used to measure the Q3 parameters that describe the arrangement of chlorine ions around the
223 : sodium atoms in sodium chloride. The imagined system here is composed of 64 NaCl formula units and the atoms are arranged in the input
224 : with the 64 Na$^+$ ions followed by the 64 Cl$-$ ions. Once again the average Q3 parameter is calculated and output to a
225 : file called colvar
226 :
227 : ```plumed
228 : q3: Q3 SPECIESA=1-64 SPECIESB=65-128 D_0=1.3 R_0=0.2
229 : q3_mean: MEAN ARG=q3 PERIODIC=NO
230 : PRINT ARG=q3_mean FILE=colvar
231 : ```
232 :
233 : If you simply want to examine the values of the Q3 parameters for each of the atoms in your system you can do so by exploiting the
234 : command [DUMPATOMS](DUMPATOMS.md) as shown in the example below. The following output file will output a file in an extended xyz format
235 : called q3.xyz for each frame of the analyzed MD trajectory. The first column in this file will contain a dummy name for each of the
236 : atoms, columns 2-4 will then contain the x, y and z positions of the atoms, column 5 will contain the value of the Q3 parameter, columns
237 : 6-12 will contain the real parts of the director of the $q_{3m}$ vector while columns 13-19 will contain the imaginary parts of this director.
238 :
239 : ```plumed
240 : q3: Q3 SPECIESA=1-64 SPECIESB=65-128 D_0=1.3 R_0=0.2 MEAN
241 : DUMPATOMS ATOMS=q3 ARG=q3 FILE=q3.xyz
242 : ```
243 :
244 : ## The MASK keyword
245 :
246 : 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
247 : input, which tells Q3 that it is safe not to calculate the Q3 parameter for some of the atoms. As illustrated below, this is useful if you are using functionality
248 : from the [volumes module](module_volumes.md) to calculate the average value of the Q3 parameter for only those atoms that lie in a certain part of the simulation box.
249 :
250 : ```plumed
251 : # Fixed virtual atom which serves as the probe volume's center (pos. in nm)
252 : center: FIXEDATOM AT=2.5,2.5,2.5
253 : # Vector in which element i is one if atom i is in sphere of interest and zero otherwise
254 : sphere: INSPHERE ATOMS=1-400 CENTER=center RADIUS={GAUSSIAN D_0=0.5 R_0=0.01 D_MAX=0.52}
255 : # Calculate the fccubic parameter of the atoms
256 : cc: Q3 ...
257 : SPECIES=1-400 MASK=sphere
258 : SWITCH={RATIONAL D_0=3.0 R_0=1.5 D_MAX=6.0}
259 : ...
260 : # Multiply fccubic parameters numbers by sphere vector
261 : prod: CUSTOM ARG=cc,sphere FUNC=x*y PERIODIC=NO
262 : # Sum of coordination numbers for atoms that are in the sphere of interest
263 : numer: SUM ARG=prod PERIODIC=NO
264 : # Number of atoms that are in sphere of interest
265 : denom: SUM ARG=sphere PERIODIC=NO
266 : # Average coordination number for atoms in sphere of interest
267 : av: CUSTOM ARG=numer,denom FUNC=x/y PERIODIC=NO
268 : # And print out final CV to a file
269 : PRINT ARG=av FILE=colvar STRIDE=1
270 : ```
271 :
272 : This input calculate the average value of the Q3 parameter for only those atoms that are within a spherical region that is centered on the point
273 : $(2.5,2.5,2.5)$.
274 :
275 : ## Deprecated syntax
276 :
277 : More information on the deprecated keywords that are given below is available in the documentation for the [DISTANCES](DISTANCES.md) command.
278 : Below is an example where these deprecated keywords are used to calculate the histogram of Q3 parameters for the 64 atoms in a box of Lennard Jones print them
279 : to a file called colvar:
280 :
281 : ```plumed
282 : q3: Q3 SPECIES=1-64 D_0=1.3 R_0=0.2 HISTOGRAM={GAUSSIAN LOWER=0.0 UPPER=1.0 NBINS=20 SMEAR=0.1}
283 : PRINT ARG=q3.* FILE=colvar
284 : ```
285 :
286 : The following example illustrates how you can use VSUM to calculate a global vector of $Q_{3m}$ values as follows:
287 :
288 : $$
289 : Q_{3m} = \sum_i \frac{q_{3m}(i)}{\sum_j \sigma(r_{ij})}
290 : $$
291 :
292 : where the sum runs over all the atoms. You can then take these $Q_{3m}$ values and compute the following norm:
293 :
294 : $$
295 : s = \sqrt{ \sum_{m=-3}^3 Q_{3m}^{*} Q_{3m} }
296 : $$
297 :
298 : The VMEAN command that is also used in the input below performs a similar operations. The only difference is that
299 : we divide the sums in the first expression above by the number of atoms.
300 :
301 : ```plumed
302 : q3: Q3 SPECIES=1-64 D_0=1.3 R_0=0.2 VMEAN VSUM
303 : PRINT ARG=q3.* FILE=colvar
304 : ```
305 :
306 : */
307 : //+ENDPLUMEDOC
308 :
309 : //+PLUMEDOC MCOLVAR Q4
310 : /*
311 : Calculate fourth order Steinhardt parameters.
312 :
313 : The 4th order Steinhardt parameters allow us to measure the degree to which the first coordination shell
314 : around an atom is ordered. The Steinhardt parameter for atom, $i$ is complex vector whose components are
315 : calculated using the following formula:
316 :
317 : $$
318 : q_{4m}(i) = \sum_j \sigma( r_{ij} ) Y_{4m}(\mathbf{r}_{ij})
319 : $$
320 :
321 : where $Y_{4m}$ is one of the 4th order spherical harmonics so $m$ is a number that runs from $-4$ to
322 : $+4$. The function $\sigma( r_{ij} )$ is a switching function that acts on the distance between
323 : atoms $i$ and $j$. The parameters of this function should be set so that it the function is equal to one
324 : when atom $j$ is in the first coordination sphere of atom $i$ and is zero otherwise.
325 :
326 : As discussed on [this page](https://www.plumed-tutorials.org/lessons/23/001/data/Steinhardt.html), the Steinhardt parameters can
327 : be used to measure the degree of order in the system in a variety of different ways. The
328 : simplest way of measuring whether or not the coordination sphere is ordered is to simply take the norm of the above vector i.e.
329 :
330 : $$
331 : Q_4(i) = \frac{1}{\sum_j \sigma(r_{ij}) } \sqrt{ \sum_{m=-4}^4 q_{4m}(i)^{*} q_{4m}(i) }
332 : $$
333 :
334 : This norm is small when the coordination shell is disordered and larger when the coordination shell is ordered. The following input illustrates
335 : how by averaging the value of this norm over all the atoms in the system you can measure the global degree of order in the system:
336 :
337 : ```plumed
338 : q4: Q4 SPECIES=1-64 D_0=1.3 R_0=0.2
339 : q4_mean: MEAN ARG=q4 PERIODIC=NO
340 : PRINT ARG=q4_mean FILE=colvar
341 : ```
342 :
343 : In the above input the rational [switching function](LESS_THAN.md) with the parameters above. We would recommend using SWITCH syntax
344 : rather than the syntax above when giving the parameters for the switching function as you can then use any of the switching functions described
345 : 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
346 : switching function as demonstrated below:
347 :
348 : ```plumed
349 : q4: Q4 SPECIES=1-64 SWITCH={RATIONAL D_0=1.3 R_0=0.2 D_MAX=3.0}
350 : q4_mean: MEAN ARG=q4 PERIODIC=NO
351 : PRINT ARG=q4_mean FILE=colvar
352 : ```
353 :
354 : 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
355 : of the documentation for [CONTACT_MATRIX](CONTACT_MATRIX.md).
356 :
357 : ## Working with two types of atoms
358 :
359 : The command below could be used to measure the Q4 parameters that describe the arrangement of chlorine ions around the
360 : sodium atoms in sodium chloride. The imagined system here is composed of 64 NaCl formula units and the atoms are arranged in the input
361 : with the 64 Na$^+$ ions followed by the 64 Cl$-$ ions. Once again the average Q4 parameter is calculated and output to a
362 : file called colvar
363 :
364 : ```plumed
365 : q4: Q4 SPECIESA=1-64 SPECIESB=65-128 D_0=1.3 R_0=0.2
366 : q4_mean: MEAN ARG=q4 PERIODIC=NO
367 : PRINT ARG=q4_mean FILE=colvar
368 : ```
369 :
370 : If you simply want to examine the values of the Q4 parameters for each of the atoms in your system you can do so by exploiting the
371 : command [DUMPATOMS](DUMPATOMS.md) as shown in the example below. The following output file will output a file in an extended xyz format
372 : called q4.xyz for each frame of the analyzed MD trajectory. The first column in this file will contain a dummy name for each of the
373 : atoms, columns 2-4 will then contain the x, y and z positions of the atoms, column 5 will contain the value of the Q4 parameter, columns
374 : 6-14 will contain the real parts of the director of the $q_{4m}$ vector while columns 15-23 will contain the imaginary parts of this director.
375 :
376 : ```plumed
377 : q4: Q4 SPECIESA=1-64 SPECIESB=65-128 D_0=1.3 R_0=0.2 MEAN
378 : DUMPATOMS ATOMS=q4 ARG=q4 FILE=q4.xyz
379 : ```
380 :
381 : ## The MASK keyword
382 :
383 : 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
384 : input, which tells Q4 that it is safe not to calculate the Q4 parameter for some of the atoms. As illustrated below, this is useful if you are using functionality
385 : from the [volumes module](module_volumes.md) to calculate the average value of the Q4 parameter for only those atoms that lie in a certain part of the simulation box.
386 :
387 : ```plumed
388 : # Fixed virtual atom which serves as the probe volume's center (pos. in nm)
389 : center: FIXEDATOM AT=2.5,2.5,2.5
390 : # Vector in which element i is one if atom i is in sphere of interest and zero otherwise
391 : sphere: INSPHERE ATOMS=1-400 CENTER=center RADIUS={GAUSSIAN D_0=0.5 R_0=0.01 D_MAX=0.52}
392 : # Calculate the fccubic parameter of the atoms
393 : cc: Q4 ...
394 : SPECIES=1-400 MASK=sphere
395 : SWITCH={RATIONAL D_0=3.0 R_0=1.5 D_MAX=6.0}
396 : ...
397 : # Multiply fccubic parameters numbers by sphere vector
398 : prod: CUSTOM ARG=cc,sphere FUNC=x*y PERIODIC=NO
399 : # Sum of coordination numbers for atoms that are in the sphere of interest
400 : numer: SUM ARG=prod PERIODIC=NO
401 : # Number of atoms that are in sphere of interest
402 : denom: SUM ARG=sphere PERIODIC=NO
403 : # Average coordination number for atoms in sphere of interest
404 : av: CUSTOM ARG=numer,denom FUNC=x/y PERIODIC=NO
405 : # And print out final CV to a file
406 : PRINT ARG=av FILE=colvar STRIDE=1
407 : ```
408 :
409 : This input calculate the average value of the Q4 parameter for only those atoms that are within a spherical region that is centered on the point
410 : $(2.5,2.5,2.5)$.
411 :
412 : ## Deprecated syntax
413 :
414 : More information on the deprecated keywords that are given below is available in the documentation for the [DISTANCES](DISTANCES.md) command.
415 : Below is an example where these deprecated keywords are used to calculate the histogram of Q4 parameters for the 64 atoms in a box of Lennard Jones print them
416 : to a file called colvar:
417 :
418 : ```plumed
419 : q4: Q4 SPECIES=1-64 D_0=1.3 R_0=0.2 HISTOGRAM={GAUSSIAN LOWER=0.0 UPPER=1.0 NBINS=20 SMEAR=0.1}
420 : PRINT ARG=q4.* FILE=colvar
421 : ```
422 :
423 : The following example illustrates how you can use VSUM to calculate a global vector of $Q_{4m}$ values as follows:
424 :
425 : $$
426 : Q_{4m} = \sum_i \frac{q_{4m}(i)}{\sum_j \sigma(r_{ij})}
427 : $$
428 :
429 : where the sum runs over all the atoms. You can then take these $Q_{4m}$ values and compute the following norm:
430 :
431 : $$
432 : s = \sqrt{ \sum_{m=-4}^1 Q_{4m}^{*} Q_{4m} }
433 : $$
434 :
435 : The VMEAN command that is also used in the input below performs a similar operations. The only difference is that
436 : we divide the sums in the first expression above by the number of atoms.
437 :
438 : ```plumed
439 : q4: Q4 SPECIES=1-64 D_0=1.3 R_0=0.2 VMEAN VSUM
440 : PRINT ARG=q4.* FILE=colvar
441 : ```
442 :
443 : */
444 : //+ENDPLUMEDOC
445 :
446 : //+PLUMEDOC MCOLVAR Q6
447 : /*
448 : Calculate sixth order Steinhardt parameters.
449 :
450 : The 6th order Steinhardt parameters allow us to measure the degree to which the first coordination shell
451 : around an atom is ordered. The Steinhardt parameter for atom, $i$ is complex vector whose components are
452 : calculated using the following formula:
453 :
454 : $$
455 : q_{6m}(i) = \sum_j \sigma( r_{ij} ) Y_{6m}(\mathbf{r}_{ij})
456 : $$
457 :
458 : where $Y_{6m}$ is one of the 6th order spherical harmonics so $m$ is a number that runs from $-6$ to
459 : $+6$. The function $\sigma( r_{ij} )$ is a switching function that acts on the distance between
460 : atoms $i$ and $j$. The parameters of this function should be set so that it the function is equal to one
461 : when atom $j$ is in the first coordination sphere of atom $i$ and is zero otherwise.
462 :
463 : As discussed on [this page](https://www.plumed-tutorials.org/lessons/23/001/data/Steinhardt.html), the Steinhardt parameters can
464 : be used to measure the degree of order in the system in a variety of different ways. The
465 : simplest way of measuring whether or not the coordination sphere is ordered is to simply take the norm of the above vector i.e.
466 :
467 : $$
468 : Q_6(i) = \frac{1}{\sum_j \sigma(r_{ij}) } \sqrt{ \sum_{m=-6}^6 q_{6m}(i)^{*} q_{6m}(i) }
469 : $$
470 :
471 : This norm is small when the coordination shell is disordered and larger when the coordination shell is ordered. The following input illustrates
472 : how by averaging the value of this norm over all the atoms in the system you can measure the global degree of order in the system:
473 :
474 : ```plumed
475 : q6: Q6 SPECIES=1-64 D_0=1.3 R_0=0.2
476 : q6_mean: MEAN ARG=q6 PERIODIC=NO
477 : PRINT ARG=q6_mean FILE=colvar
478 : ```
479 :
480 : In the above input the rational [switching function](LESS_THAN.md) with the parameters above. We would recommend using SWITCH syntax
481 : rather than the syntax above when giving the parameters for the switching function as you can then use any of the switching functions described
482 : 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
483 : switching function as demonstrated below:
484 :
485 : ```plumed
486 : q6: Q6 SPECIES=1-64 SWITCH={RATIONAL D_0=1.3 R_0=0.2 D_MAX=3.0}
487 : q6_mean: MEAN ARG=q6 PERIODIC=NO
488 : PRINT ARG=q6_mean FILE=colvar
489 : ```
490 :
491 : 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
492 : of the documentation for [CONTACT_MATRIX](CONTACT_MATRIX.md).
493 :
494 : ## Working with two types of atoms
495 :
496 : The command below could be used to measure the Q6 parameters that describe the arrangement of chlorine ions around the
497 : sodium atoms in sodium chloride. The imagined system here is composed of 64 NaCl formula units and the atoms are arranged in the input
498 : with the 64 Na$^+$ ions followed by the 64 Cl$-$ ions. Once again the average Q6 parameter is calculated and output to a
499 : file called colvar
500 :
501 : ```plumed
502 : q6: Q6 SPECIESA=1-64 SPECIESB=65-128 D_0=1.3 R_0=0.2
503 : q6_mean: MEAN ARG=q6 PERIODIC=NO
504 : PRINT ARG=q6_mean FILE=colvar
505 : ```
506 :
507 : If you simply want to examine the values of the Q6 parameters for each of the atoms in your system you can do so by exploiting the
508 : command [DUMPATOMS](DUMPATOMS.md) as shown in the example below. The following output file will output a file in an extended xyz format
509 : called q6.xyz for each frame of the analyzed MD trajectory. The first column in this file will contain a dummy name for each of the
510 : atoms, columns 2-4 will then contain the x, y and z positions of the atoms, column 5 will contain the value of the Q6 parameter, columns
511 : 6-18 will contain the real parts of the director of the $q_{6m}$ vector while columns 19-31 will contain the imaginary parts of this director.
512 :
513 : ```plumed
514 : q6: Q6 SPECIESA=1-64 SPECIESB=65-128 D_0=1.3 R_0=0.2 MEAN
515 : DUMPATOMS ATOMS=q6 ARG=q6 FILE=q6.xyz
516 : ```
517 :
518 : ## The MASK keyword
519 :
520 : 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
521 : input, which tells Q6 that it is safe not to calculate the Q6 parameter for some of the atoms. As illustrated below, this is useful if you are using functionality
522 : from the [volumes module](module_volumes.md) to calculate the average value of the Q6 parameter for only those atoms that lie in a certain part of the simulation box.
523 :
524 : ```plumed
525 : # Fixed virtual atom which serves as the probe volume's center (pos. in nm)
526 : center: FIXEDATOM AT=2.5,2.5,2.5
527 : # Vector in which element i is one if atom i is in sphere of interest and zero otherwise
528 : sphere: INSPHERE ATOMS=1-400 CENTER=center RADIUS={GAUSSIAN D_0=0.5 R_0=0.01 D_MAX=0.52}
529 : # Calculate the fccubic parameter of the atoms
530 : cc: Q6 ...
531 : SPECIES=1-400 MASK=sphere
532 : SWITCH={RATIONAL D_0=3.0 R_0=1.5 D_MAX=6.0}
533 : ...
534 : # Multiply fccubic parameters numbers by sphere vector
535 : prod: CUSTOM ARG=cc,sphere FUNC=x*y PERIODIC=NO
536 : # Sum of coordination numbers for atoms that are in the sphere of interest
537 : numer: SUM ARG=prod PERIODIC=NO
538 : # Number of atoms that are in sphere of interest
539 : denom: SUM ARG=sphere PERIODIC=NO
540 : # Average coordination number for atoms in sphere of interest
541 : av: CUSTOM ARG=numer,denom FUNC=x/y PERIODIC=NO
542 : # And print out final CV to a file
543 : PRINT ARG=av FILE=colvar STRIDE=1
544 : ```
545 :
546 : This input calculate the average value of the Q6 parameter for only those atoms that are within a spherical region that is centered on the point
547 : $(2.5,2.5,2.5)$.
548 :
549 : ## Deprecated syntax
550 :
551 : More information on the deprecated keywords that are given below is available in the documentation for the [DISTANCES](DISTANCES.md) command.
552 : Below is an example where these deprecated keywords are used to calculate the histogram of Q6 parameters for the 64 atoms in a box of Lennard Jones print them
553 : to a file called colvar:
554 :
555 : ```plumed
556 : q6: Q6 SPECIES=1-64 D_0=1.3 R_0=0.2 HISTOGRAM={GAUSSIAN LOWER=0.0 UPPER=1.0 NBINS=20 SMEAR=0.1}
557 : PRINT ARG=q6.* FILE=colvar
558 : ```
559 :
560 : The following example illustrates how you can use VSUM to calculate a global vector of $Q_{6m}$ values as follows:
561 :
562 : $$
563 : Q_{6m} = \sum_i \frac{q_{6m}(i)}{\sum_j \sigma(r_{ij})}
564 : $$
565 :
566 : where the sum runs over all the atoms. You can then take these $Q_{6m}$ values and compute the following norm:
567 :
568 : $$
569 : s = \sqrt{ \sum_{m=-6}^6 Q_{6m}^{*} Q_{6m} }
570 : $$
571 :
572 : The VMEAN command that is also used in the input below performs a similar operations. The only difference is that
573 : we divide the sums in the first expression above by the number of atoms.
574 :
575 : ```plumed
576 : q6: Q6 SPECIES=1-64 D_0=1.3 R_0=0.2 VMEAN VSUM
577 : PRINT ARG=q6.* FILE=colvar
578 : ```
579 :
580 : */
581 : //+ENDPLUMEDOC
582 :
583 : class Steinhardt : public ActionShortcut {
584 : private:
585 : static std::string getSymbol( int m );
586 : void createVectorNormInput( const std::string& ilab,
587 : const std::string& olab,
588 : int l,
589 : const std::string& sep,
590 : const std::string& vlab,
591 : bool usegpu = false);
592 : public:
593 : static void registerKeywords( Keywords& keys );
594 : explicit Steinhardt(const ActionOptions&);
595 : };
596 :
597 : PLUMED_REGISTER_ACTION(Steinhardt,"Q1")
598 : PLUMED_REGISTER_ACTION(Steinhardt,"Q3")
599 : PLUMED_REGISTER_ACTION(Steinhardt,"Q4")
600 : PLUMED_REGISTER_ACTION(Steinhardt,"Q6")
601 :
602 60 : void Steinhardt::registerKeywords( Keywords& keys ) {
603 60 : CoordinationNumbers::shortcutKeywords( keys );
604 120 : keys.addDeprecatedFlag("LOWMEM","");
605 60 : keys.addFlag("VMEAN",false,"calculate the norm of the mean vector.");
606 120 : keys.addOutputComponent("_vmean","VMEAN","scalar","the norm of the mean vector");
607 60 : keys.addFlag("VSUM",false,"calculate the norm of the sum of all the vectors");
608 120 : keys.addOutputComponent("_vsum","VSUM","scalar","the norm of the mean vector");
609 60 : keys.needsAction("GROUP");
610 60 : keys.needsAction("CONTACT_MATRIX");
611 60 : keys.needsAction("SPHERICAL_HARMONIC");
612 60 : keys.needsAction("ONES");
613 60 : keys.needsAction("MATRIX_VECTOR_PRODUCT");
614 60 : keys.needsAction("COMBINE");
615 60 : keys.needsAction("CUSTOM");
616 60 : keys.needsAction("MEAN");
617 60 : keys.needsAction("SUM");
618 120 : keys.setValueDescription("vector",
619 : "the norms of the vectors of spherical harmonic coefficients");
620 60 : keys.addFlag("USEGPU",false,"run part of this calculation on the GPU");
621 60 : }
622 :
623 42 : Steinhardt::Steinhardt( const ActionOptions& ao):
624 : Action(ao),
625 42 : ActionShortcut(ao) {
626 : {
627 : bool lowmem;
628 42 : parseFlag("LOWMEM",lowmem);
629 42 : if( lowmem ) {
630 0 : warning("LOWMEM flag is deprecated and is no longer required for this action");
631 : }
632 : }
633 :
634 : bool usegpu;
635 42 : parseFlag("USEGPU",usegpu);
636 84 : const std::string doUSEGPU = usegpu?" USEGPU":"";
637 :
638 : std::string sp_str;
639 84 : parse("SPECIES",sp_str);
640 : std::string specA;
641 84 : parse("SPECIESA",specA);
642 : std::string specB;
643 42 : parse("SPECIESB",specB);
644 42 : CoordinationNumbers::expandMatrix( true,
645 : getShortcutLabel(),
646 : sp_str,
647 : specA,
648 : specB,
649 : this );
650 : int l;
651 84 : std::string sph_input = getShortcutLabel() + "_sh: SPHERICAL_HARMONIC ARG="
652 126 : + getShortcutLabel() + "_mat.x,"
653 126 : + getShortcutLabel() + "_mat.y,"
654 126 : + getShortcutLabel() + "_mat.z,"
655 126 : + getShortcutLabel() + "_mat.w";
656 :
657 42 : if( getName()=="Q1" ) {
658 : sph_input +=" L=1";
659 : l=1;
660 23 : } else if( getName()=="Q3" ) {
661 : sph_input += " L=3";
662 : l=3;
663 22 : } else if( getName()=="Q4" ) {
664 : sph_input += " L=4";
665 : l=4;
666 19 : } else if( getName()=="Q6" ) {
667 : sph_input += " L=6";
668 : l=6;
669 : } else {
670 0 : plumed_merror("invalid input");
671 : }
672 42 : readInputLine( sph_input + doUSEGPU);
673 :
674 : // Input for denominator (coord)
675 42 : ActionWithValue* av = plumed.getActionSet()
676 42 : .selectWithLabel<ActionWithValue*>( getShortcutLabel() + "_mat");
677 42 : plumed_assert( av && av->getNumberOfComponents()>0 && (av->copyOutput(0))->getRank()==2 );
678 : std::string size;
679 42 : Tools::convert( (av->copyOutput(0))->getShape()[1], size );
680 84 : readInputLine( getShortcutLabel() + "_denom_ones: ONES SIZE=" + size );
681 84 : readInputLine( getShortcutLabel() + "_denom: MATRIX_VECTOR_PRODUCT "
682 126 : "ARG=" + getShortcutLabel() + "_mat.w,"
683 126 : + getShortcutLabel() + "_denom_ones"
684 84 : + doUSEGPU);
685 126 : readInputLine( getShortcutLabel() + "_sp: MATRIX_VECTOR_PRODUCT "
686 126 : "ARG=" + getShortcutLabel() + "_sh.*,"
687 126 : + getShortcutLabel() + "_denom_ones"
688 84 : + doUSEGPU);
689 :
690 : // If we are doing VMEAN determine sum of vector components
691 : std::string snum;
692 : bool do_vmean;
693 42 : parseFlag("VMEAN",do_vmean);
694 : bool do_vsum;
695 42 : parseFlag("VSUM",do_vsum);
696 42 : if( do_vmean || do_vsum ) {
697 26 : auto makeString=[&](const std::string& numstr,
698 : const char realImg)->std::string{
699 : //realImg is "r" or "i"
700 52 : return getShortcutLabel() + "_" +realImg + "mn-" + numstr + ": CUSTOM "
701 78 : "ARG=" + getShortcutLabel() + "_sp." +realImg + "m-" + numstr + ","
702 104 : + getShortcutLabel() + "_denom FUNC=x/y PERIODIC=NO";
703 1 : };
704 : // Divide all components by coordination numbers
705 14 : for(int i=-l; i<=l; ++i) {
706 13 : snum = getSymbol( i );
707 : // Real part
708 13 : readInputLine(makeString(snum,'r'));
709 : // Imaginary part
710 26 : readInputLine(makeString(snum,'i'));
711 : }
712 : }
713 :
714 42 : if( do_vmean ) {
715 26 : auto makeString=[&](const std::string& numstr,
716 : const char realImg)->std::string{
717 : //realImg is "r" or "i"
718 52 : return getShortcutLabel() + "_" +realImg + "ms-" + numstr + ": MEAN "
719 78 : "ARG=" + getShortcutLabel() + "_" +realImg + "mn-" + numstr
720 52 : + " PERIODIC=NO";
721 1 : };
722 14 : for(int i=-l; i<=l; ++i) {
723 13 : snum = getSymbol( i );
724 : // Real part
725 13 : readInputLine(makeString(snum,'r'));
726 : // Imaginary part
727 26 : readInputLine(makeString(snum,'i'));
728 : }
729 : // Now calculate the total length of the vector
730 3 : createVectorNormInput( getShortcutLabel(),
731 2 : getShortcutLabel() + "_vmean",
732 : l,
733 : "_",
734 : "ms",
735 : usegpu );
736 : }
737 42 : if( do_vsum ) {
738 0 : auto makeString=[&](const std::string& numstr,
739 : const std::string& realImg)->std::string{
740 0 : return getShortcutLabel() + "_" +realImg + "mz-" + numstr + ": SUM "
741 0 : "ARG=" + getShortcutLabel() + "_" + realImg + "mn-" + numstr
742 0 : + " PERIODIC=NO";
743 0 : };
744 0 : for(int i=-l; i<=l; ++i) {
745 0 : snum = getSymbol( i );
746 : // Real part
747 0 : readInputLine(makeString(snum,"r"));
748 : // Imaginary part
749 0 : readInputLine(makeString(snum,"i"));
750 : }
751 : // Now calculate the total length of the vector
752 0 : createVectorNormInput( getShortcutLabel(),
753 0 : getShortcutLabel() + "_vsum",
754 : l,
755 : "_",
756 : "mz",
757 : usegpu );
758 : }
759 :
760 : // Now calculate the total length of the vector
761 168 : createVectorNormInput( getShortcutLabel() + "_sp",
762 84 : getShortcutLabel() + "_norm",
763 : l,
764 : ".",
765 : "m",
766 : usegpu );
767 : // And take average
768 84 : readInputLine( getShortcutLabel() + ": CUSTOM "
769 126 : "ARG=" + getShortcutLabel() + "_norm,"
770 126 : + getShortcutLabel() + "_denom FUNC=x/y PERIODIC=NO");
771 84 : multicolvar::MultiColvarShortcuts::expandFunctions( getShortcutLabel(),
772 : getShortcutLabel(), "", this );
773 42 : }
774 :
775 43 : void Steinhardt::createVectorNormInput( const std::string& ilab,
776 : const std::string& olab,
777 : const int l,
778 : const std::string& sep,
779 : const std::string& vlab,
780 : const bool usegpu) {
781 43 : std::string norm_input = olab + "2: COMBINE PERIODIC=NO POWERS=";
782 43 : std::string snum = getSymbol( -l );
783 43 : std::string arg_inp = "";
784 :
785 86 : std::string arg_inp_real = ilab + sep + "r" + vlab + "-";
786 86 : std::string arg_inp_img = ilab + sep + "i" + vlab + "-";
787 43 : std::string comma="";
788 394 : for(int i=-l; i<=l; ++i) {
789 351 : snum = getSymbol( i );
790 702 : arg_inp += comma + arg_inp_real + snum + "," + arg_inp_img + snum;
791 702 : norm_input += comma + "2,2";
792 : comma=",";
793 : }
794 129 : readInputLine( norm_input + " ARG=" + arg_inp + (usegpu?" USEGPU":"") );
795 86 : readInputLine( olab + ": CUSTOM ARG=" + olab + "2 FUNC=sqrt(x) PERIODIC=NO");
796 43 : }
797 :
798 420 : std::string Steinhardt::getSymbol( const int m ) {
799 420 : if( m<0 ) {
800 : std::string num;
801 209 : Tools::convert( -m, num );
802 209 : return "n" + num;
803 211 : } else if( m>0 ) {
804 : std::string num;
805 166 : Tools::convert( m, num );
806 166 : return "p" + num;
807 : }
808 45 : return "0";
809 : }
810 :
811 : } // namespace symfunc
812 : } // namespace PLMD
|
