Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2014-2020 The plumed team
3 : (see the PEOPLE file at the root of the distribution for a list of names)
4 :
5 : See http://www.plumed.org for more information.
6 :
7 : This file is part of plumed, version 2.
8 :
9 : plumed is free software: you can redistribute it and/or modify
10 : it under the terms of the GNU Lesser General Public License as published by
11 : the Free Software Foundation, either version 3 of the License, or
12 : (at your option) any later version.
13 :
14 : plumed is distributed in the hope that it will be useful,
15 : but WITHOUT ANY WARRANTY; without even the implied warranty of
16 : MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
17 : GNU Lesser General Public License for more details.
18 :
19 : You should have received a copy of the GNU Lesser General Public License
20 : along with plumed. If not, see <http://www.gnu.org/licenses/>.
21 : +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
22 : #include "CoordinationNumbers.h"
23 : #include "multicolvar/MultiColvarShortcuts.h"
24 : #include "core/ActionWithValue.h"
25 : #include "core/ActionShortcut.h"
26 : #include "core/ActionRegister.h"
27 : #include "core/PlumedMain.h"
28 : #include "core/ActionSet.h"
29 :
30 : #include <string>
31 : #include <cmath>
32 :
33 : namespace PLMD {
34 : namespace symfunc {
35 :
36 : //+PLUMEDOC MCOLVAR COORDINATION_SHELL_FUNCTION
37 : /*
38 : Calculate an arbitrary function of all the bond vectors in the first coordination sphere of an atom
39 :
40 : This shortcut allows you to calculate the sum for an arbitrary function of the bond vectors that connect an atom to each of its neighbours.
41 : In other words, this action allows you to compute the following:
42 :
43 : $$
44 : s_i = \sum_{i \ne j} \sigma(r_{ij}) f(x_{ij}, y_{ij}, z_{ij}, r_{ij}) )
45 : $$
46 :
47 : In this expression, $x_{ij}, y_{ij}, z_{ij}$ are the components of the vector connecting atoms $i$ and $j$ and $r_{ij}$ is the magnitude of this vector.
48 : $\sigma(r_{ij})$ is then a switching function that ensures that the aritrary function $f$ is only evaluated for if atom $j$ is within a certain cutoff distance
49 : of atom $i$.
50 :
51 : Below you can see a simple example that shows how this action can be used in practise.
52 :
53 : ```plumed
54 : cv: COORDINATION_SHELL_FUNCTION SPECIES=1-64 D_0=3.0 R_0=1.5 FUNCTION=((x+y+z)/r)^3+((x-y-z)/r)^3+((-x+y-z)/r)^3+((-x-y+z)/r)^3
55 : PRINT ARG=cv FILE=colvar
56 : ```
57 :
58 : The above input calculates 64 $s_i$ values - one $s_i$ values for each of the atoms specified using the SPECIES keyword. These 64 numbers are then output to a file.
59 : The function of x, y, z and r to be evaluated is specified using the FUNCTION keyword. Obviously, if your function does not depend on all four of these variables
60 : they can be excluded from your function.
61 :
62 : In the input above we use a rational [switching function](LESS_THAN.md) with the parameters above. We would recommend using SWITCH syntax
63 : rather than the syntax above when giving the parameters for the switching function as you can then use any of the switching functions described
64 : 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
65 : switching function as demonstrated below:
66 :
67 : ```plumed
68 : cv: COORDINATION_SHELL_FUNCTION ...
69 : SPECIES=1-64 SWITCH={RATIONAL D_0=3.0 R_0=1.5 D_MAX=7.0}
70 : FUNCTION=((x+y+z)/r)^3+((x-y-z)/r)^3+((-x+y-z)/r)^3+((-x-y+z)/r)^3
71 : ...
72 : PRINT ARG=cv FILE=colvar
73 : ```
74 :
75 : 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
76 : of the documentation for [CONTACT_MATRIX](CONTACT_MATRIX.md).
77 :
78 : ## Rotating the reference frame
79 :
80 : As discussed in [this paper](https://journals.aps.org/prb/abstract/10.1103/PhysRevB.92.180102) it is sometimes useful to rotate the bond vectors before computing the
81 : arbitrary function $f$ in the above expression. To perform such rotations you use the PHI, THETA and PSI keywords. The $x_{ij}, y_{ij}$ and $z_{ij}$ values that enter $f$ in the
82 : expression above are then calculated as:
83 :
84 : $$
85 : \left(
86 : \begin{matrix}
87 : x_{ij} \\
88 : y_{ij} \\
89 : z_{ij}
90 : \end{matrix}
91 : \right) =
92 : \left(
93 : \begin{matrix}
94 : \cos(\psi)\cos(\phi) - \cos(\theta)\sin(\phi)\sin(\psi) & \cos(\psi)*\sin(\phi)+\cos(\theta)*\cos(\phi)*\sin(\psi) & \sin(\psi)*sin(\theta) \\
95 : -\sin(\psi)*\cos(\phi)-\cos(\theta)*\sin(\phi)*\cos(\psi) & -\sin(\psi)*\sin(\phi)+\cos(\theta)*\cos(\phi)*std::cos(\psi), & \cos(\psi)*\sin(\theta) \\
96 : \sin(\theta)*\sin(\phi) & \sin(\theta)*\cos(\phi) & \cos(\theta)
97 : \end{matrix}
98 : \right)
99 : \left(
100 : \begin{matrix}
101 : x_{ij}' \\
102 : y_{ij}' \\
103 : z_{ij}'
104 : \end{matrix}
105 : \right)
106 : $$
107 :
108 : $x_{ij}', y_{ij}'$ and $z_{ij}'$ in this expression are the bond vectors that are calculated in the lab frame. The matrix in the above expression is thus just a
109 : [rotation matrix](https://en.wikipedia.org/wiki/Rotation_matrix) that converts the lab frame to some frame of interest.
110 :
111 : ## Working with two types of atom
112 :
113 : If you would like to calculate a function of the bonds connecting the atoms in GROUPA to the atoms in GROUPB you can use an input like the one
114 : shown below:
115 :
116 : ```plumed
117 : d: COORDINATION_SHELL_FUNCTION ...
118 : SPECIESA=1-64 SPECIESB=65-200
119 : FUNCTION=((x+y+z)/r)^3+((x-y-z)/r)^3+((-x+y-z)/r)^3+((-x-y+z)/r)^3
120 : SWITCH={RATIONAL D_0=3.0 R_0=1.5 D_MAX=6.0}
121 : ...
122 : lt: MORE_THAN ARG=d SWITCH={RATIONAL R_0=0.5}
123 : s: SUM ARG=lt PERIODIC=NO
124 : PRINT ARG=s FILE=colv
125 : ```
126 :
127 : ## The MASK keyword
128 :
129 : 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
130 : 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
131 : from the [volumes module](module_volumes.md) to calculate the average value of the TETRAHEDRAL parameter for only those atoms that lie in a certain part of the simulation box.
132 :
133 : ```plumed
134 : # Fixed virtual atom which serves as the probe volume's center (pos. in nm)
135 : center: FIXEDATOM AT=2.5,2.5,2.5
136 : # Vector in which element i is one if atom i is in sphere of interest and zero otherwise
137 : sphere: INSPHERE ATOMS=1-400 CENTER=center RADIUS={GAUSSIAN D_0=0.5 R_0=0.01 D_MAX=0.52}
138 : # Calculate the tetrahedral parameter of the atoms
139 : d: COORDINATION_SHELL_FUNCTION ...
140 : SPECIES=1-400 MASK=sphere
141 : FUNCTION=((x+y+z)/r)^3+((x-y-z)/r)^3+((-x+y-z)/r)^3+((-x-y+z)/r)^3
142 : SWITCH={RATIONAL D_0=3.0 R_0=1.5 D_MAX=6.0}
143 : ...
144 : # Multiply fccubic parameters numbers by sphere vector
145 : prod: CUSTOM ARG=d,sphere FUNC=x*y PERIODIC=NO
146 : # Sum of coordination numbers for atoms that are in the sphere of interest
147 : numer: SUM ARG=prod PERIODIC=NO
148 : # Number of atoms that are in sphere of interest
149 : denom: SUM ARG=sphere PERIODIC=NO
150 : # Average coordination number for atoms in sphere of interest
151 : av: CUSTOM ARG=numer,denom FUNC=x/y PERIODIC=NO
152 : # And print out final CV to a file
153 : PRINT ARG=av FILE=colvar STRIDE=1
154 : ```
155 :
156 : This input calculate the average value of the function for only those atoms that are within a spherical region that is centered on the point
157 : $(2.5,2.5,2.5)$.
158 :
159 : ## Deprecated syntax
160 :
161 : More information on the deprecated keywords that are given below is available in the documentation for the [DISTANCES](DISTANCES.md) command.
162 :
163 : */
164 : //+ENDPLUMEDOC
165 :
166 : //+PLUMEDOC MCOLVAR COORDINATION_SHELL_AVERAGE
167 : /*
168 : Calculate an arbitrary function of all the bond vectors in the first coordination sphere of an atom and take an average
169 :
170 : This shortcut allows you to calculate the average for an arbitrary function of the bond vectors that connect an atom to each of its neighbours.
171 : In other words, this action allows you to compute the following:
172 :
173 : $$
174 : s_i = \frac{\sum_{i \ne j} \sigma(r_{ij}) f(x_{ij}, y_{ij}, z_{ij}, r_{ij}) )}{\sum_{i \ne j} \sigma(r_{ij})}
175 : $$
176 :
177 : In this expression, $x_{ij}, y_{ij}, z_{ij}$ are the components of the vector connecting atoms $i$ and $j$ and $r_{ij}$ is the magnitude of this vector.
178 : $\sigma(r_{ij})$ is then a switching function that ensures that the aritrary function $f$ is only evaluated for if atom $j$ is within a certain cutoff distance
179 : of atom $i$.
180 :
181 : Below you can see a simple example that shows how this action can be used in practise.
182 :
183 : ```plumed
184 : cv: COORDINATION_SHELL_AVERAGE SPECIES=1-64 D_0=3.0 R_0=1.5 FUNCTION=((x+y+z)/r)^3+((x-y-z)/r)^3+((-x+y-z)/r)^3+((-x-y+z)/r)^3
185 : PRINT ARG=cv FILE=colvar
186 : ```
187 :
188 : The above input calculates 64 $s_i$ values - one $s_i$ values for each of the atoms specified using the SPECIES keyword. These 64 numbers are then output to a file.
189 : The function of x, y, z and r to be evaluated is specified using the FUNCTION keyword. Obviously, if your function does not depend on all four of these variables
190 : they can be excluded from your function.
191 :
192 : In the input above we use a rational [switching function](LESS_THAN.md) with the parameters above. We would recommend using SWITCH syntax
193 : rather than the syntax above when giving the parameters for the switching function as you can then use any of the switching functions described
194 : 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
195 : switching function as demonstrated below:
196 :
197 : ```plumed
198 : cv: COORDINATION_SHELL_AVERAGE ...
199 : SPECIES=1-64 SWITCH={RATIONAL D_0=3.0 R_0=1.5 D_MAX=7.0}
200 : FUNCTION=((x+y+z)/r)^3+((x-y-z)/r)^3+((-x+y-z)/r)^3+((-x-y+z)/r)^3
201 : ...
202 : PRINT ARG=cv FILE=colvar
203 : ```
204 :
205 : 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
206 : of the documentation for [CONTACT_MATRIX](CONTACT_MATRIX.md).
207 :
208 : Notice that you can you can rotate the bond vectors before computing the
209 : arbitrary function $f$ in the above expression as is discussed in the documentation for [COORDINATION_SHELL_FUNCTION](COORDINATION_SHELL_FUNCTION.md)
210 :
211 : ## Working with two types of atom
212 :
213 : If you would like to calculate a function of the bonds connecting the atoms in GROUPA to the atoms in GROUPB you can use an input like the one
214 : shown below:
215 :
216 : ```plumed
217 : d: COORDINATION_SHELL_AVERAGE ...
218 : SPECIESA=1-64 SPECIESB=65-200
219 : FUNCTION=((x+y+z)/r)^3+((x-y-z)/r)^3+((-x+y-z)/r)^3+((-x-y+z)/r)^3
220 : SWITCH={RATIONAL D_0=3.0 R_0=1.5 D_MAX=6.0}
221 : ...
222 : lt: MORE_THAN ARG=d SWITCH={RATIONAL R_0=0.5}
223 : s: SUM ARG=lt PERIODIC=NO
224 : PRINT ARG=s FILE=colv
225 : ```
226 :
227 : ## The MASK keyword
228 :
229 : 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
230 : 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
231 : from the [volumes module](module_volumes.md) to calculate the average value of the TETRAHEDRAL parameter for only those atoms that lie in a certain part of the simulation box.
232 :
233 : ```plumed
234 : # Fixed virtual atom which serves as the probe volume's center (pos. in nm)
235 : center: FIXEDATOM AT=2.5,2.5,2.5
236 : # Vector in which element i is one if atom i is in sphere of interest and zero otherwise
237 : sphere: INSPHERE ATOMS=1-400 CENTER=center RADIUS={GAUSSIAN D_0=0.5 R_0=0.01 D_MAX=0.52}
238 : # Calculate the function for the atoms
239 : d: COORDINATION_SHELL_AVERAGE ...
240 : SPECIES=1-400 MASK=sphere
241 : FUNCTION=((x+y+z)/r)^3+((x-y-z)/r)^3+((-x+y-z)/r)^3+((-x-y+z)/r)^3
242 : SWITCH={RATIONAL D_0=3.0 R_0=1.5 D_MAX=6.0}
243 : ...
244 : # Multiply fccubic parameters numbers by sphere vector
245 : prod: CUSTOM ARG=d,sphere FUNC=x*y PERIODIC=NO
246 : # Sum of coordination numbers for atoms that are in the sphere of interest
247 : numer: SUM ARG=prod PERIODIC=NO
248 : # Number of atoms that are in sphere of interest
249 : denom: SUM ARG=sphere PERIODIC=NO
250 : # Average coordination number for atoms in sphere of interest
251 : av: CUSTOM ARG=numer,denom FUNC=x/y PERIODIC=NO
252 : # And print out final CV to a file
253 : PRINT ARG=av FILE=colvar STRIDE=1
254 : ```
255 :
256 : This input calculate the average value of the function for only those atoms that are within a spherical region that is centered on the point
257 : $(2.5,2.5,2.5)$.
258 :
259 : ## Deprecated syntax
260 :
261 : More information on the deprecated keywords that are given below is available in the documentation for the [DISTANCES](DISTANCES.md) command.
262 :
263 : */
264 : //+ENDPLUMEDOC
265 :
266 : //+PLUMEDOC MCOLVAR SIMPLECUBIC
267 : /*
268 : Calculate whether or not the coordination spheres of atoms are arranged as they would be in a simple cubic structure.
269 :
270 : This shortcut is an example of a [COORDINATION_SHELL_AVERAGE](COORDINATION_SHELL_AVERAGE.md),
271 : which we can use to measure how similar the environment around atom $i$ is to a simple cubic structure. We perform this comparison by evaluating
272 : the following quantity:
273 :
274 : $$
275 : s_i = \frac{ \sum_{i \ne j} \sigma(r_{ij}) \left[ \frac{ x_{ij}^4 + y_{ij}^4 + z_{ij}^4 }{r_{ij}^4} \right] }{ \sum_{i \ne j} \sigma(r_{ij}) }
276 : $$
277 :
278 : In this expression $x_{ij}$, $y_{ij}$ and $z_{ij}$ are the $x$, $y$ and $z$ components of the vector connecting atom $i$ to
279 : atom $j$ and $r_{ij}$ is the magnitude of this vector. $\sigma(r_{ij})$ is a switching function that acts on the distance between atom $i$ and atom $j$ and its inclusion in the numerator and the denominator of the above expression as well as the fact that we are summing
280 : over all of the other atoms in the system ensures that we are calculating an average
281 : of the function of $x_{ij}$, $y_{ij}$ and $z_{ij}$ for the atoms in the first coordination sphere around atom $i$.
282 : This quantity is once again a multicolvar so you can compute it for multiple atoms using a single PLUMED action and then compute
283 : the average value for the atoms in your system, the number of atoms that have an $s_i$ value that is more that some target and
284 : so on. Notice also that you can rotate the reference frame if you are using a non-standard unit cell.
285 :
286 : The following input tells plumed to calculate the simple cubic parameter for the atoms 1-100 with themselves.
287 : The mean value is then calculated.
288 :
289 : ```plumed
290 : sc: SIMPLECUBIC SPECIES=1-100 R_0=1.0
291 : sc_mean: MEAN ARG=sc PERIODIC=NO
292 : PRINT ARG=sc_mean FILE=colvar
293 : ```
294 :
295 : In the input above we use a rational [switching function](LESS_THAN.md) with the parameters above. We would recommend using SWITCH syntax
296 : rather than the syntax above when giving the parameters for the switching function as you can then use any of the switching functions described
297 : 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
298 : switching function as demonstrated below:
299 :
300 : ```plumed
301 : d: SIMPLECUBIC SPECIES=1-64 SWITCH={RATIONAL D_0=3.0 R_0=1.5 D_MAX=6.0}
302 : dm: MEAN ARG=d PERIODIC=NO
303 : PRINT ARG=dm FILE=colv
304 : ```
305 :
306 : 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
307 : of the documentation for [CONTACT_MATRIX](CONTACT_MATRIX.md).
308 :
309 : ## Working with two types of atom
310 :
311 : If you would like to calculate whether the atoms in GROUPB are arranged around the atoms in GROUPA as they in an tetrahedral structure you use an input like the one
312 : shown below:
313 :
314 : ```plumed
315 : d: SIMPLECUBIC SPECIESA=1-64 SPECIESB=65-200 SWITCH={RATIONAL D_0=3.0 R_0=1.5 D_MAX=6.0}
316 : lt: MORE_THAN ARG=d SWITCH={RATIONAL R_0=0.5}
317 : s: SUM ARG=lt PERIODIC=NO
318 : PRINT ARG=s FILE=colv
319 : ```
320 :
321 : ## The MASK keyword
322 :
323 : 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
324 : input, which tells SIMPLECUBIC that it is safe not to calculate the SIMPLECUBIC parameter for some of the atoms. As illustrated below, this is useful if you are using functionality
325 : from the [volumes module](module_volumes.md) to calculate the average value of the SIMPLECUBIC parameter for only those atoms that lie in a certain part of the simulation box.
326 :
327 : ```plumed
328 : # Fixed virtual atom which serves as the probe volume's center (pos. in nm)
329 : center: FIXEDATOM AT=2.5,2.5,2.5
330 : # Vector in which element i is one if atom i is in sphere of interest and zero otherwise
331 : sphere: INSPHERE ATOMS=1-400 CENTER=center RADIUS={GAUSSIAN D_0=0.5 R_0=0.01 D_MAX=0.52}
332 : # Calculate the tetrahedral parameter of the atoms
333 : cc: SIMPLECUBIC ...
334 : SPECIES=1-400 MASK=sphere
335 : SWITCH={RATIONAL D_0=3.0 R_0=1.5 D_MAX=6.0}
336 : ...
337 : # Multiply fccubic parameters numbers by sphere vector
338 : prod: CUSTOM ARG=cc,sphere FUNC=x*y PERIODIC=NO
339 : # Sum of coordination numbers for atoms that are in the sphere of interest
340 : numer: SUM ARG=prod PERIODIC=NO
341 : # Number of atoms that are in sphere of interest
342 : denom: SUM ARG=sphere PERIODIC=NO
343 : # Average coordination number for atoms in sphere of interest
344 : av: CUSTOM ARG=numer,denom FUNC=x/y PERIODIC=NO
345 : # And print out final CV to a file
346 : PRINT ARG=av FILE=colvar STRIDE=1
347 : ```
348 :
349 : This input calculate the average value of the SIMPLECUBIC parameter for only those atoms that are within a spherical region that is centered on the point
350 : $(2.5,2.5,2.5)$.
351 :
352 : ## Deprecated syntax
353 :
354 : More information on the deprecated keywords that are given below is available in the documentation for the [DISTANCES](DISTANCES.md) command.
355 :
356 : The following command illustrates how this deprecated syntax works. It tells plumed to look at the ways atoms 1-100 are within 3.0 are arranged about atoms
357 : from 101-110. The number of simple cubic parameters that are greater than 0.8 is then output
358 :
359 : ```plumed
360 : SIMPLECUBIC SPECIESA=101-110 SPECIESB=1-100 R_0=3.0 MORE_THAN={RATIONAL R_0=0.8 NN=6 MM=12 D_0=0}
361 : ```
362 :
363 : Notice that you can you can rotate the bond vectors before computing the
364 : function in the above expression as is discussed in the documentation for [COORDINATION_SHELL_FUNCTION](COORDINATION_SHELL_FUNCTION.md)
365 :
366 : */
367 : //+ENDPLUMEDOC
368 :
369 : //+PLUMEDOC MCOLVAR TETRAHEDRAL
370 : /*
371 : Calculate the degree to which the environment about ions has a tetrahedral order.
372 :
373 : This shortcut is an example of a [COORDINATION_SHELL_AVERAGE](COORDINATION_SHELL_AVERAGE.md),
374 : which we can use to measure the degree to which the atoms in the first coordination shell around any atom, $i$ is
375 : is arranged like a tetrahedron. We perform this comparison by evaluating the following function.
376 :
377 : $$
378 : s(i) = \frac{1}{\sum_j \sigma( r_{ij} )} \sum_j \sigma( r_{ij} )\left[ \frac{(x_{ij} + y_{ij} + z_{ij})^3}{r_{ij}^3} +
379 : \frac{(x_{ij} - y_{ij} - z_{ij})^3}{r_{ij}^3} +
380 : \frac{(-x_{ij} + y_{ij} - z_{ij})^3}{r_{ij}^3} +
381 : \frac{(-x_{ij} - y_{ij} + z_{ij})^3}{r_{ij}^3} \right]
382 : $$
383 :
384 : Here $r_{ij}$ is the magnitude of the vector connecting atom $i$ to atom $j$ and $x_{ij}$, $y_{ij}$ and $z_{ij}$
385 : are its three components. The function $\sigma( r_{ij} )$ is a switching function that acts on the distance between
386 : atoms $i$ and $j$. The parameters of this function should be set so that the function is equal to one
387 : when atom $j$ is in the first coordination sphere of atom $i$ and is zero otherwise.
388 :
389 : The following command calculates the average value of the TETRAHEDRAL parameter for a set of 64 atoms all of the same type
390 : and outputs this quantity to a file called colvar.
391 :
392 : ```plumed
393 : tt: TETRAHEDRAL SPECIES=1-64 D_0=3.0 R_0=1.5 NN=6 MM=12
394 : tt_mean: MEAN ARG=tt PERIODIC=NO
395 : PRINT ARG=tt_mean FILE=colvar
396 : ```
397 :
398 : In the input above we use a rational [switching function](LESS_THAN.md) with the parameters above. We would recommend using SWITCH syntax
399 : rather than the syntax above when giving the parameters for the switching function as you can then use any of the switching functions described
400 : 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
401 : switching function as demonstrated below:
402 :
403 : ```plumed
404 : d: TETRAHEDRAL SPECIES=1-64 SWITCH={RATIONAL D_0=3.0 R_0=1.5 D_MAX=6.0}
405 : dm: MEAN ARG=d PERIODIC=NO
406 : PRINT ARG=dm FILE=colv
407 : ```
408 :
409 : 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
410 : of the documentation for [CONTACT_MATRIX](CONTACT_MATRIX.md).
411 :
412 : ## Working with two types of atom
413 :
414 : If you would like to calculate whether the atoms in GROUPB are arranged around the atoms in GROUPA as they in an tetrahedral structure you use an input like the one
415 : shown below:
416 :
417 : ```plumed
418 : d: TETRAHEDRAL SPECIESA=1-64 SPECIESB=65-200 SWITCH={RATIONAL D_0=3.0 R_0=1.5 D_MAX=6.0}
419 : lt: MORE_THAN ARG=d SWITCH={RATIONAL R_0=0.5}
420 : s: SUM ARG=lt PERIODIC=NO
421 : PRINT ARG=s FILE=colv
422 : ```
423 :
424 : ## The MASK keyword
425 :
426 : 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
427 : input, which tells TETRAHEDRAL that it is safe not to calculate the TETRAHEDRAL parameter for some of the atoms. As illustrated below, this is useful if you are using functionality
428 : from the [volumes module](module_volumes.md) to calculate the average value of the TETRAHEDRAL parameter for only those atoms that lie in a certain part of the simulation box.
429 :
430 : ```plumed
431 : # Fixed virtual atom which serves as the probe volume's center (pos. in nm)
432 : center: FIXEDATOM AT=2.5,2.5,2.5
433 : # Vector in which element i is one if atom i is in sphere of interest and zero otherwise
434 : sphere: INSPHERE ATOMS=1-400 CENTER=center RADIUS={GAUSSIAN D_0=0.5 R_0=0.01 D_MAX=0.52}
435 : # Calculate the tetrahedral parameter of the atoms
436 : cc: TETRAHEDRAL ...
437 : SPECIES=1-400 MASK=sphere
438 : SWITCH={RATIONAL D_0=3.0 R_0=1.5 D_MAX=6.0}
439 : ...
440 : # Multiply fccubic parameters numbers by sphere vector
441 : prod: CUSTOM ARG=cc,sphere FUNC=x*y PERIODIC=NO
442 : # Sum of coordination numbers for atoms that are in the sphere of interest
443 : numer: SUM ARG=prod PERIODIC=NO
444 : # Number of atoms that are in sphere of interest
445 : denom: SUM ARG=sphere PERIODIC=NO
446 : # Average coordination number for atoms in sphere of interest
447 : av: CUSTOM ARG=numer,denom FUNC=x/y PERIODIC=NO
448 : # And print out final CV to a file
449 : PRINT ARG=av FILE=colvar STRIDE=1
450 : ```
451 :
452 : This input calculate the average value of the TETRAHEDRAL parameter for only those atoms that are within a spherical region that is centered on the point
453 : $(2.5,2.5,2.5)$.
454 :
455 : ## Deprecated syntax
456 :
457 : More information on the deprecated keywords that are given below is available in the documentation for the [DISTANCES](DISTANCES.md) command.
458 :
459 : The following command illustrates how this deprecated syntax works. It calculates the number of TETRAHEDRAL parameters that are greater than 0.8 in a set of 10 atoms.
460 : In this calculation it is assumed that there are two atom types A and B and that the first coordination sphere of the
461 : 10 atoms of type A contains atoms of type B. The formula above is thus calculated for ten different A atoms and within
462 : it the sum over $j$ runs over 40 atoms of type B that could be in the first coordination sphere.
463 :
464 : ```plumed
465 : tt: TETRAHEDRAL SPECIESA=1-10 SPECIESB=11-40 SWITCH={RATIONAL D_0=1.3 R_0=0.2} MORE_THAN={RATIONAL R_0=0.8}
466 : PRINT ARG=tt.* FILE=colvar
467 : ```
468 :
469 : Notice that you can you can rotate the bond vectors before computing the
470 : function in the above expression as is discussed in the documentation for [COORDINATION_SHELL_FUNCTION](COORDINATION_SHELL_FUNCTION.md)
471 :
472 : */
473 : //+ENDPLUMEDOC
474 :
475 : class CoordShellVectorFunction : public ActionShortcut {
476 : public:
477 : static void registerKeywords(Keywords& keys);
478 : explicit CoordShellVectorFunction(const ActionOptions&);
479 : };
480 :
481 : PLUMED_REGISTER_ACTION(CoordShellVectorFunction,"FCCUBIC")
482 : PLUMED_REGISTER_ACTION(CoordShellVectorFunction,"TETRAHEDRAL")
483 : PLUMED_REGISTER_ACTION(CoordShellVectorFunction,"SIMPLECUBIC")
484 : PLUMED_REGISTER_ACTION(CoordShellVectorFunction,"COORDINATION_SHELL_FUNCTION")
485 : PLUMED_REGISTER_ACTION(CoordShellVectorFunction,"COORDINATION_SHELL_AVERAGE")
486 :
487 33 : void CoordShellVectorFunction::registerKeywords( Keywords& keys ) {
488 33 : CoordinationNumbers::shortcutKeywords( keys );
489 60 : if( keys.getDisplayName()=="COORDINATION_SHELL_FUNCTION" || keys.getDisplayName()=="COORDINATION_SHELL_AVERAGE" ) {
490 10 : keys.add("compulsory","FUNCTION","the function of the bond vectors that you would like to evaluate");
491 : }
492 33 : keys.add("compulsory","PHI","0.0","The Euler rotational angle phi");
493 33 : keys.add("compulsory","THETA","0.0","The Euler rotational angle theta");
494 33 : keys.add("compulsory","PSI","0.0","The Euler rotational angle psi");
495 33 : if( keys.getDisplayName()=="FCCUBIC" ) {
496 12 : keys.add("compulsory","ALPHA","3.0","The alpha parameter of the angular function that is used for FCCUBIC");
497 : }
498 66 : keys.addDeprecatedFlag("LOWMEM","");
499 66 : keys.setValueDescription("vector","the symmetry function for each of the specified atoms");
500 33 : keys.needsAction("CONTACT_MATRIX");
501 33 : keys.needsAction("FCCUBIC_FUNC");
502 33 : keys.needsAction("CUSTOM");
503 33 : keys.needsAction("ONES");
504 66 : keys.needsAction("MATRIX_VECTOR_PRODUCT");
505 33 : keys.addFlag("USEGPU",false,"run part of this calculation on the GPU");
506 33 : keys.addDOI("10.1103/PhysRevB.81.125416");
507 33 : keys.addDOI("10.1103/PhysRevB.92.180102");
508 33 : keys.addDOI("10.1063/1.4997180");
509 33 : keys.addDOI("10.1063/1.5134461");
510 33 : }
511 :
512 9 : CoordShellVectorFunction::CoordShellVectorFunction(const ActionOptions& ao):
513 : Action(ao),
514 9 : ActionShortcut(ao) {
515 : #define createLabel(name) const std::string name##Lab = getShortcutLabel()+"_"#name;
516 : bool lowmem;
517 9 : parseFlag("LOWMEM",lowmem);
518 9 : if( lowmem ) {
519 0 : warning("LOWMEM flag is deprecated and is no longer required for this action");
520 : }
521 : bool usegpu;
522 9 : parseFlag("USEGPU",usegpu);
523 18 : const std::string doUSEGPU = usegpu?" USEGPU":"";
524 : std::string sp_str;
525 18 : parse("SPECIES",sp_str);
526 : std::string specA;
527 18 : parse("SPECIESA",specA);
528 : std::string specB;
529 9 : parse("SPECIESB",specB);
530 9 : createLabel(mat);
531 9 : if( sp_str.length()>0 || specA.length()>0 ) {
532 9 : CoordinationNumbers::expandMatrix( true, getShortcutLabel(), sp_str, specA, specB, this );
533 : } else {
534 0 : error("found no input atoms use SPECIES/SPECIESA");
535 : }
536 : double phi;
537 9 : parse("PHI",phi);
538 : double theta;
539 9 : parse("THETA",theta);
540 : double psi;
541 9 : parse("PSI",psi);
542 9 : std::vector<std::string> rotelements(9);
543 9 : std::string xvec = matLab + ".x", yvec = matLab + ".y", zvec = matLab + ".z";
544 9 : if( phi!=0 || theta!=0 || psi!=0 ) {
545 0 : Tools::convert( std::cos(psi)*std::cos(phi)
546 0 : - std::cos(theta)*std::sin(phi)*std::sin(psi),
547 : rotelements[0] );
548 0 : Tools::convert( std::cos(psi)*std::sin(phi)
549 0 : + std::cos(theta)*std::cos(phi)*std::sin(psi),
550 : rotelements[1] );
551 0 : Tools::convert( std::sin(psi)*std::sin(theta), rotelements[2] );
552 :
553 0 : Tools::convert( -std::sin(psi)*std::cos(phi)
554 0 : - std::cos(theta)*std::sin(phi)*std::cos(psi),
555 : rotelements[3] );
556 0 : Tools::convert( -std::sin(psi)*std::sin(phi)
557 0 : + std::cos(theta)*std::cos(phi)*std::cos(psi),
558 : rotelements[4] );
559 0 : Tools::convert( std::cos(psi)*std::sin(theta), rotelements[5] );
560 :
561 0 : Tools::convert( std::sin(theta)*std::sin(phi), rotelements[6] );
562 0 : Tools::convert( -std::sin(theta)*std::cos(phi), rotelements[7] );
563 0 : Tools::convert( std::cos(theta), rotelements[8] );
564 0 : createLabel(xrot);
565 0 : createLabel(yrot);
566 0 : createLabel(zrot);
567 0 : readInputLine( xrotLab + ": CUSTOM ARG=" + matLab + ".x," + matLab + ".y," + matLab + ".z "
568 0 : "FUNC=" + rotelements[0] + "*x+" + rotelements[1] + "*y+" + rotelements[2] + "*z "
569 : "PERIODIC=NO");
570 0 : readInputLine( yrotLab + ": CUSTOM ARG=" + matLab + ".x," + matLab + ".y," + matLab + ".z "
571 0 : "FUNC=" + rotelements[3] + "*x+" + rotelements[4] + "*y+" + rotelements[5] + "*z "
572 : "PERIODIC=NO");
573 0 : readInputLine( zrotLab + ": CUSTOM ARG=" + matLab + ".x," + matLab + ".y," + matLab + ".z "
574 0 : "FUNC=" + rotelements[6] + "*x+" + rotelements[7] + "*y+" + rotelements[8] + "*z "
575 : "PERIODIC=NO");
576 : }
577 : // Calculate FCC cubic function from bond vectors
578 9 : createLabel(vfunc)
579 9 : createLabel(r);
580 9 : if( getName()=="FCCUBIC" ) {
581 : std::string alpha;
582 5 : parse("ALPHA",alpha);
583 5 : readInputLine( vfuncLab + ": FCCUBIC_FUNC "
584 25 : "ARG=" + xvec + "," + yvec + "," + zvec+ " ALPHA=" + alpha);
585 4 : } else if( getName()=="TETRAHEDRAL" ) {
586 4 : readInputLine( rLab + ": CUSTOM ARG=" + xvec + "," + yvec + "," + zvec + " "
587 : "PERIODIC=NO FUNC=sqrt(x*x+y*y+z*z)");
588 1 : readInputLine( vfuncLab + ": CUSTOM "
589 4 : "ARG=" + xvec + "," + yvec + "," + zvec + "," + rLab
590 1 : + " VAR=x,y,z,r PERIODIC=NO "
591 : "FUNC=((x+y+z)/r)^3+((x-y-z)/r)^3+((-x+y-z)/r)^3+((-x-y+z)/r)^3" );
592 3 : } else if( getName()=="SIMPLECUBIC" ) {
593 4 : readInputLine( rLab + ": CUSTOM ARG=" + xvec + "," + yvec + "," + zvec
594 1 : + " PERIODIC=NO FUNC=sqrt(x*x+y*y+z*z)");
595 1 : readInputLine( vfuncLab + ": CUSTOM "
596 4 : "ARG=" + xvec + "," + yvec + "," + zvec + "," + rLab
597 1 : + " VAR=x,y,z,r PERIODIC=NO FUNC=(x^4+y^4+z^4)/(r^4)" );
598 : } else {
599 : std::string myfunc;
600 2 : parse("FUNCTION",myfunc);
601 2 : if( myfunc.find("r")!=std::string::npos ) {
602 8 : readInputLine( rLab + ": CUSTOM ARG=" + xvec + "," + yvec + "," + zvec
603 2 : + " PERIODIC=NO FUNC=sqrt(x*x+y*y+z*z)");
604 2 : readInputLine( vfuncLab + ": CUSTOM "
605 8 : "ARG=" + xvec + "," + yvec + "," + zvec + "," + rLab
606 4 : + " VAR=x,y,z,r PERIODIC=NO FUNC=" + myfunc );
607 : } else {
608 0 : readInputLine( vfuncLab + ": CUSTOM ARG=" + xvec + "," + yvec + "," + zvec
609 0 : + " PERIODIC=NO FUNC=" + myfunc );
610 : }
611 : }
612 : // Hadamard product of function above and weights
613 9 : createLabel(wvfunc);
614 27 : readInputLine( wvfuncLab + ": CUSTOM ARG=" + vfuncLab + "," + matLab + ".w "
615 : "FUNC=x*y PERIODIC=NO");
616 : // And coordination numbers
617 9 : ActionWithValue* av = plumed.getActionSet().selectWithLabel<ActionWithValue*>( matLab );
618 9 : plumed_assert( av && av->getNumberOfComponents()>0 && (av->copyOutput(0))->getRank()==2 );
619 : std::string size;
620 9 : Tools::convert( (av->copyOutput(0))->getShape()[1], size );
621 9 : createLabel(ones);
622 9 : readInputLine( onesLab + ": ONES SIZE=" + size );
623 18 : readInputLine( getShortcutLabel() + ": MATRIX_VECTOR_PRODUCT "
624 18 : "ARG=" + wvfuncLab + "," + onesLab + doUSEGPU);
625 9 : std::string olab=getShortcutLabel();
626 9 : if( getName()!="COORDINATION_SHELL_FUNCTION" ) {
627 7 : olab = getShortcutLabel() + "_n";
628 : // Calculate coordination numbers for denominator
629 7 : createLabel(denom);
630 7 : readInputLine( denomLab + ": MATRIX_VECTOR_PRODUCT "
631 14 : "ARG=" + matLab + ".w," + onesLab + doUSEGPU);
632 : // And normalise
633 7 : readInputLine( olab + ": CUSTOM "
634 21 : "ARG=" + getShortcutLabel() + "," + denomLab + " FUNC=x/y PERIODIC=NO");
635 : }
636 : // And expand the functions
637 : std::map<std::string,std::string> keymap;
638 9 : multicolvar::MultiColvarShortcuts::readShortcutKeywords( keymap, this );
639 18 : multicolvar::MultiColvarShortcuts::expandFunctions( getShortcutLabel(),
640 : olab, "", keymap, this );
641 : #undef createLabel
642 18 : }
643 :
644 : }
645 : }
646 :
|
