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 SWITCH={RATIONAL 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 : 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
63 : 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
64 : expression above are then calculated as:
65 :
66 : $$
67 : \left(
68 : \begin{matrix}
69 : x_{ij} \\
70 : y_{ij} \\
71 : z_{ij}
72 : \end{matrix}
73 : \right) =
74 : \left(
75 : \begin{matrix}
76 : \cos(\psi)\cos(\phi) - \cos(\theta)\sin(\phi)\sin(\psi) & \cos(\psi)*\sin(\phi)+\cos(\theta)*\cos(\phi)*\sin(\psi) & \sin(\psi)*sin(\theta) \\
77 : -\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) \\
78 : \sin(\theta)*\sin(\phi) & \sin(\theta)*\cos(\phi) & \cos(\theta)
79 : \end{matrix}
80 : \left(
81 : \begin{matrix}
82 : x_{ij}' \\
83 : y_{ij}' \\
84 : z_{ij}'
85 : \end{matrix}
86 : \right)
87 : $$
88 :
89 : $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
90 : [rotation matrix](https://en.wikipedia.org/wiki/Rotation_matrix) that converts the lab frame to some frame of interest.
91 :
92 :
93 : */
94 : //+ENDPLUMEDOC
95 :
96 : //+PLUMEDOC MCOLVAR COORDINATION_SHELL_AVERAGE
97 : /*
98 : Calculate an arbitrary function of all the bond vectors in the first coordination sphere of an atom and take an average
99 :
100 : 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.
101 : In other words, this action allows you to compute the following:
102 :
103 : $$
104 : 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})}
105 : $$
106 :
107 : 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.
108 : $\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
109 : of atom $i$.
110 :
111 : Below you can see a simple example that shows how this action can be used in practise.
112 :
113 : ```plumed
114 : cv: COORDINATION_SHELL_AVERAGE SPECIES=1-64 SWITCH={RATIONAL 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
115 : PRINT ARG=cv FILE=colvar
116 : ```
117 :
118 : 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.
119 : 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
120 : they can be excluded from your function.
121 :
122 : Notice that you can you can rotate the bond vectors before computing the
123 : arbitrary function $f$ in the above expression as is discussed in the documentation for [COORDINATION_SHELL_FUNCTION](COORDINATION_SHELL_FUNCTION.md)
124 :
125 : */
126 : //+ENDPLUMEDOC
127 :
128 : //+PLUMEDOC MCOLVAR SIMPLECUBIC
129 : /*
130 : Calculate whether or not the coordination spheres of atoms are arranged as they would be in a simple cubic structure.
131 :
132 : This shortcut is an example of a [COORDINATION_SHELL_AVERAGE](COORDINATION_SHELL_AVERAGE.md),
133 : 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
134 : the following quantity:
135 :
136 : $$
137 : 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}) }
138 : $$
139 :
140 : In this expression $x_{ij}$, $y_{ij}$ and $z_{ij}$ are the $x$, $y$ and $z$ components of the vector connecting atom $i$ to
141 : 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
142 : over all of the other atoms in the system ensures that we are calculating an average
143 : of the function of $x_{ij}$, $y_{ij}$ and $z_{ij}$ for the atoms in the first coordination sphere around atom $i$.
144 : This quantity is once again a multicolvar so you can compute it for multiple atoms using a single PLUMED action and then compute
145 : 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
146 : so on. Notice also that you can rotate the reference frame if you are using a non-standard unit cell.
147 :
148 : The following input tells plumed to calculate the simple cubic parameter for the atoms 1-100 with themselves.
149 : The mean value is then calculated.
150 :
151 : ```plumed
152 : SIMPLECUBIC SPECIES=1-100 R_0=1.0 MEAN
153 : ```
154 :
155 : The following input tells plumed to look at the ways atoms 1-100 are within 3.0 are arranged about atoms
156 : from 101-110. The number of simple cubic parameters that are greater than 0.8 is then output
157 :
158 : ```plumed
159 : 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}
160 : ```
161 :
162 : Notice that you can you can rotate the bond vectors before computing the
163 : function in the above expression as is discussed in the documentation for [COORDINATION_SHELL_FUNCTION](COORDINATION_SHELL_FUNCTION.md)
164 :
165 : */
166 : //+ENDPLUMEDOC
167 :
168 : //+PLUMEDOC MCOLVAR TETRAHEDRAL
169 : /*
170 : Calculate the degree to which the environment about ions has a tetrahedral order.
171 :
172 : This shortcut is an example of a [COORDINATION_SHELL_AVERAGE](COORDINATION_SHELL_AVERAGE.md),
173 : which we can use to measure the degree to which the atoms in the first coordination shell around any atom, $i$ is
174 : is arranged like a tetrahedron. We perform this comparison by evaluating the following function.
175 :
176 : $$
177 : 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} +
178 : \frac{(x_{ij} - y_{ij} - z_{ij})^3}{r_{ij}^3} +
179 : \frac{(-x_{ij} + y_{ij} - z_{ij})^3}{r_{ij}^3} +
180 : \frac{(-x_{ij} - y_{ij} + z_{ij})^3}{r_{ij}^3} \right]
181 : $$
182 :
183 : Here $r_{ij}$ is the magnitude of the vector connecting atom $i$ to atom $j$ and $x_{ij}$, $y_{ij}$ and $z_{ij}$
184 : are its three components. The function $\sigma( r_{ij} )$ is a switching function that acts on the distance between
185 : atoms $i$ and $j$. The parameters of this function should be set so that the function is equal to one
186 : when atom $j$ is in the first coordination sphere of atom $i$ and is zero otherwise.
187 :
188 : The following command calculates the average value of the TETRAHEDRAL parameter for a set of 64 atoms all of the same type
189 : and outputs this quantity to a file called colvar.
190 :
191 : ```plumed
192 : tt: TETRAHEDRAL SPECIES=1-64 SWITCH={RATIONAL D_0=1.3 R_0=0.2} MEAN
193 : PRINT ARG=tt.mean FILE=colvar
194 : ```
195 :
196 : The following command calculates the number of TETRAHEDRAL parameters that are greater than 0.8 in a set of 10 atoms.
197 : In this calculation it is assumed that there are two atom types A and B and that the first coordination sphere of the
198 : 10 atoms of type A contains atoms of type B. The formula above is thus calculated for ten different A atoms and within
199 : it the sum over $j$ runs over 40 atoms of type B that could be in the first coordination sphere.
200 :
201 : ```plumed
202 : 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}
203 : PRINT ARG=tt.* FILE=colvar
204 : ```
205 :
206 : Notice that you can you can rotate the bond vectors before computing the
207 : function in the above expression as is discussed in the documentation for [COORDINATION_SHELL_FUNCTION](COORDINATION_SHELL_FUNCTION.md)
208 :
209 : */
210 : //+ENDPLUMEDOC
211 :
212 : class CoordShellVectorFunction : public ActionShortcut {
213 : public:
214 : static void registerKeywords(Keywords& keys);
215 : explicit CoordShellVectorFunction(const ActionOptions&);
216 : };
217 :
218 : PLUMED_REGISTER_ACTION(CoordShellVectorFunction,"FCCUBIC")
219 : PLUMED_REGISTER_ACTION(CoordShellVectorFunction,"TETRAHEDRAL")
220 : PLUMED_REGISTER_ACTION(CoordShellVectorFunction,"SIMPLECUBIC")
221 : PLUMED_REGISTER_ACTION(CoordShellVectorFunction,"COORDINATION_SHELL_FUNCTION")
222 : PLUMED_REGISTER_ACTION(CoordShellVectorFunction,"COORDINATION_SHELL_AVERAGE")
223 :
224 23 : void CoordShellVectorFunction::registerKeywords( Keywords& keys ) {
225 23 : CoordinationNumbers::shortcutKeywords( keys );
226 23 : keys.add("compulsory","FUNCTION","the function of the bond vectors that you would like to evaluate");
227 23 : keys.add("compulsory","PHI","0.0","The Euler rotational angle phi");
228 23 : keys.add("compulsory","THETA","0.0","The Euler rotational angle theta");
229 23 : keys.add("compulsory","PSI","0.0","The Euler rotational angle psi");
230 23 : keys.add("compulsory","ALPHA","3.0","The alpha parameter of the angular function that is used for FCCUBIC");
231 23 : keys.addFlag("LOWMEM",false,"this flag does nothing and is present only to ensure back-compatibility");
232 46 : keys.setValueDescription("vector","the symmetry function for each of the specified atoms");
233 23 : keys.needsAction("CONTACT_MATRIX");
234 23 : keys.needsAction("FCCUBIC_FUNC");
235 23 : keys.needsAction("CUSTOM");
236 23 : keys.needsAction("ONES");
237 23 : keys.needsAction("MATRIX_VECTOR_PRODUCT");
238 23 : keys.addDOI("10.1103/PhysRevB.81.125416");
239 23 : keys.addDOI("10.1103/PhysRevB.92.180102");
240 23 : keys.addDOI("10.1063/1.4997180");
241 23 : keys.addDOI("10.1063/1.5134461");
242 23 : }
243 :
244 9 : CoordShellVectorFunction::CoordShellVectorFunction(const ActionOptions& ao):
245 : Action(ao),
246 9 : ActionShortcut(ao) {
247 : std::string matlab, sp_str, specA, specB;
248 : bool lowmem;
249 9 : parseFlag("LOWMEM",lowmem);
250 9 : if( lowmem ) {
251 0 : warning("LOWMEM flag is deprecated and is no longer required for this action");
252 : }
253 9 : parse("SPECIES",sp_str);
254 9 : parse("SPECIESA",specA);
255 18 : parse("SPECIESB",specB);
256 9 : if( sp_str.length()>0 || specA.length()>0 ) {
257 9 : matlab = getShortcutLabel() + "_mat";
258 9 : CoordinationNumbers::expandMatrix( true, getShortcutLabel(), sp_str, specA, specB, this );
259 : } else {
260 0 : error("found no input atoms use SPECIES/SPECIESA");
261 : }
262 : double phi, theta, psi;
263 9 : parse("PHI",phi);
264 9 : parse("THETA",theta);
265 9 : parse("PSI",psi);
266 9 : std::vector<std::string> rotelements(9);
267 9 : std::string xvec = matlab + ".x", yvec = matlab + ".y", zvec = matlab + ".z";
268 9 : if( phi!=0 || theta!=0 || psi!=0 ) {
269 0 : Tools::convert( std::cos(psi)*std::cos(phi)-std::cos(theta)*std::sin(phi)*std::sin(psi), rotelements[0] );
270 0 : Tools::convert( std::cos(psi)*std::sin(phi)+std::cos(theta)*std::cos(phi)*std::sin(psi), rotelements[1] );
271 0 : Tools::convert( std::sin(psi)*std::sin(theta), rotelements[2] );
272 :
273 0 : Tools::convert( -std::sin(psi)*std::cos(phi)-std::cos(theta)*std::sin(phi)*std::cos(psi), rotelements[3] );
274 0 : Tools::convert( -std::sin(psi)*std::sin(phi)+std::cos(theta)*std::cos(phi)*std::cos(psi), rotelements[4] );
275 0 : Tools::convert( std::cos(psi)*std::sin(theta), rotelements[5] );
276 :
277 0 : Tools::convert( std::sin(theta)*std::sin(phi), rotelements[6] );
278 0 : Tools::convert( -std::sin(theta)*std::cos(phi), rotelements[7] );
279 0 : Tools::convert( std::cos(theta), rotelements[8] );
280 0 : readInputLine( getShortcutLabel() + "_xrot: CUSTOM ARG=" + matlab + ".x," + matlab + ".y," + matlab + ".z FUNC=" + rotelements[0] + "*x+" + rotelements[1] + "*y+" + rotelements[2] + "*z PERIODIC=NO");
281 0 : readInputLine( getShortcutLabel() + "_yrot: CUSTOM ARG=" + matlab + ".x," + matlab + ".y," + matlab + ".z FUNC=" + rotelements[3] + "*x+" + rotelements[4] + "*y+" + rotelements[5] + "*z PERIODIC=NO");
282 0 : readInputLine( getShortcutLabel() + "_zrot: CUSTOM ARG=" + matlab + ".x," + matlab + ".y," + matlab + ".z FUNC=" + rotelements[6] + "*x+" + rotelements[7] + "*y+" + rotelements[8] + "*z PERIODIC=NO");
283 : }
284 : // Calculate FCC cubic function from bond vectors
285 9 : if( getName()=="FCCUBIC" ) {
286 : std::string alpha;
287 5 : parse("ALPHA",alpha);
288 10 : readInputLine( getShortcutLabel() + "_vfunc: FCCUBIC_FUNC ARG=" + xvec + "," + yvec + "," + zvec+ " ALPHA=" + alpha);
289 4 : } else if( getName()=="TETRAHEDRAL" ) {
290 2 : readInputLine( getShortcutLabel() + "_r: CUSTOM ARG=" + xvec + "," + yvec + "," + zvec + " PERIODIC=NO FUNC=sqrt(x*x+y*y+z*z)");
291 3 : readInputLine( getShortcutLabel() + "_vfunc: CUSTOM ARG=" + xvec + "," + yvec + "," + zvec + "," + getShortcutLabel() + "_r"
292 2 : + " VAR=x,y,z,r PERIODIC=NO FUNC=((x+y+z)/r)^3+((x-y-z)/r)^3+((-x+y-z)/r)^3+((-x-y+z)/r)^3" );
293 3 : } else if( getName()=="SIMPLECUBIC" ) {
294 2 : readInputLine( getShortcutLabel() + "_r: CUSTOM ARG=" + xvec + "," + yvec + "," + zvec + " PERIODIC=NO FUNC=sqrt(x*x+y*y+z*z)");
295 3 : readInputLine( getShortcutLabel() + "_vfunc: CUSTOM ARG=" + xvec + "," + yvec + "," + zvec + "," + getShortcutLabel() + "_r"
296 2 : + " VAR=x,y,z,r PERIODIC=NO FUNC=(x^4+y^4+z^4)/(r^4)" );
297 : } else {
298 : std::string myfunc;
299 2 : parse("FUNCTION",myfunc);
300 2 : if( myfunc.find("r")!=std::string::npos ) {
301 4 : readInputLine( getShortcutLabel() + "_r: CUSTOM ARG=" + xvec + "," + yvec + "," + zvec + " PERIODIC=NO FUNC=sqrt(x*x+y*y+z*z)");
302 4 : readInputLine( getShortcutLabel() + "_vfunc: CUSTOM ARG=" + xvec + "," + yvec + "," + zvec + "," + getShortcutLabel() + "_r VAR=x,y,z,r PERIODIC=NO FUNC=" + myfunc );
303 : } else {
304 0 : readInputLine( getShortcutLabel() + "_vfunc: CUSTOM ARG=" + xvec + "," + yvec + "," + zvec + " PERIODIC=NO FUNC=" + myfunc );
305 : }
306 : }
307 : // Hadamard product of function above and weights
308 18 : readInputLine( getShortcutLabel() + "_wvfunc: CUSTOM ARG=" + getShortcutLabel() + "_vfunc," + matlab + ".w FUNC=x*y PERIODIC=NO");
309 : // And coordination numbers
310 9 : ActionWithValue* av = plumed.getActionSet().selectWithLabel<ActionWithValue*>( getShortcutLabel() + "_mat");
311 9 : plumed_assert( av && av->getNumberOfComponents()>0 && (av->copyOutput(0))->getRank()==2 );
312 : std::string size;
313 9 : Tools::convert( (av->copyOutput(0))->getShape()[1], size );
314 18 : readInputLine( getShortcutLabel() + "_ones: ONES SIZE=" + size );
315 18 : readInputLine( getShortcutLabel() + ": MATRIX_VECTOR_PRODUCT ARG=" + getShortcutLabel() + "_wvfunc," + getShortcutLabel() + "_ones");
316 9 : std::string olab=getShortcutLabel();
317 9 : if( getName()!="COORDINATION_SHELL_FUNCTION" ) {
318 7 : olab = getShortcutLabel() + "_n";
319 : // Calculate coordination numbers for denominator
320 14 : readInputLine( getShortcutLabel() + "_denom: MATRIX_VECTOR_PRODUCT ARG=" + matlab + ".w," + getShortcutLabel() + "_ones");
321 : // And normalise
322 14 : readInputLine( getShortcutLabel() + "_n: CUSTOM ARG=" + getShortcutLabel() + "," + getShortcutLabel() + "_denom FUNC=x/y PERIODIC=NO");
323 : }
324 : // And expand the functions
325 : std::map<std::string,std::string> keymap;
326 9 : multicolvar::MultiColvarShortcuts::readShortcutKeywords( keymap, this );
327 18 : multicolvar::MultiColvarShortcuts::expandFunctions( getShortcutLabel(), olab, "", keymap, this );
328 18 : }
329 :
330 : }
331 : }
332 :
|
