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 : \par Examples
40 :
41 : */
42 : //+ENDPLUMEDOC
43 :
44 : //+PLUMEDOC MCOLVAR Q3
45 : /*
46 : Calculate 3rd order Steinhardt parameters.
47 :
48 : The 3rd order Steinhardt parameters allow us to measure the degree to which the first coordination shell
49 : around an atom is ordered. The Steinhardt parameter for atom, \f$i\f$ is complex vector whose components are
50 : calculated using the following formula:
51 :
52 : \f[
53 : q_{3m}(i) = \frac{\sum_j \sigma( r_{ij} ) Y_{3m}(\mathbf{r}_{ij}) }{\sum_j \sigma( r_{ij} ) }
54 : \f]
55 :
56 : where \f$Y_{3m}\f$ is one of the 3rd order spherical harmonics so \f$m\f$ is a number that runs from \f$-3\f$ to
57 : \f$+3\f$. The function \f$\sigma( r_{ij} )\f$ is a \ref switchingfunction that acts on the distance between
58 : atoms \f$i\f$ and \f$j\f$. The parameters of this function should be set so that it the function is equal to one
59 : when atom \f$j\f$ is in the first coordination sphere of atom \f$i\f$ and is zero otherwise.
60 :
61 : The Steinhardt parameters can be used to measure the degree of order in the system in a variety of different ways. The
62 : simplest way of measuring whether or not the coordination sphere is ordered is to simply take the norm of the above vector i.e.
63 :
64 : \f[
65 : Q_3(i) = \sqrt{ \sum_{m=-3}^3 q_{3m}(i)^{*} q_{3m}(i) }
66 : \f]
67 :
68 : This norm is small when the coordination shell is disordered and larger when the coordination shell is ordered. Furthermore, when
69 : the keywords LESS_THAN, MIN, MAX, HISTOGRAM, MEAN and so on are used with this colvar it is the distribution of these normed quantities
70 : that is investigated.
71 :
72 : Other measures of order can be taken by averaging the components of the individual \f$q_3\f$ vectors individually or by taking dot products of
73 : the \f$q_{3}\f$ vectors on adjacent atoms. More information on these variables can be found in the documentation for \ref LOCAL_Q3,
74 : \ref LOCAL_AVERAGE and \ref NLINKS.
75 :
76 : \par Examples
77 :
78 : The following command calculates the average Q3 parameter for the 64 atoms in a box of Lennard Jones and prints this
79 : quantity to a file called colvar:
80 :
81 : \plumedfile
82 : Q3 SPECIES=1-64 D_0=1.3 R_0=0.2 MEAN LABEL=q3
83 : PRINT ARG=q3.mean FILE=colvar
84 : \endplumedfile
85 :
86 : The following command calculates the histogram of Q3 parameters for the 64 atoms in a box of Lennard Jones and prints these
87 : quantities to a file called colvar:
88 :
89 : \plumedfile
90 : 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} LABEL=q3
91 : PRINT ARG=q3.* FILE=colvar
92 : \endplumedfile
93 :
94 : The following command could be used to measure the Q3 parameters that describe the arrangement of chlorine ions around the
95 : sodium atoms in sodium chloride. The imagined system here is composed of 64 NaCl formula units and the atoms are arranged in the input
96 : with the 64 Na\f$^+\f$ ions followed by the 64 Cl\f$-\f$ ions. Once again the average Q3 parameter is calculated and output to a
97 : file called colvar
98 :
99 : \plumedfile
100 : Q3 SPECIESA=1-64 SPECIESB=65-128 D_0=1.3 R_0=0.2 MEAN LABEL=q3
101 : PRINT ARG=q3.mean FILE=colvar
102 : \endplumedfile
103 :
104 : 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
105 : command \ref DUMPATOMS as shown in the example below. The following output file will output a file in an extended xyz format
106 : 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
107 : 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
108 : 6-12 will contain the real parts of the director of the \f$q_{3m}\f$ vector while columns 12-19 will contain the imaginary parts of this director.
109 :
110 : \plumedfile
111 : q3: Q3 SPECIESA=1-64 SPECIESB=65-128 D_0=1.3 R_0=0.2 MEAN
112 : DUMPATOMS ATOMS=q3 ARG=q3_anorm FILE=q3.xyz
113 : \endplumedfile
114 :
115 : */
116 : //+ENDPLUMEDOC
117 :
118 : //+PLUMEDOC MCOLVAR Q4
119 : /*
120 : Calculate fourth order Steinhardt parameters.
121 :
122 : The fourth order Steinhardt parameters allow us to measure the degree to which the first coordination shell
123 : around an atom is ordered. The Steinhardt parameter for atom, \f$i\f$ is complex vector whose components are
124 : calculated using the following formula:
125 :
126 : \f[
127 : q_{4m}(i) = \frac{\sum_j \sigma( r_{ij} ) Y_{4m}(\mathbf{r}_{ij}) }{\sum_j \sigma( r_{ij} ) }
128 : \f]
129 :
130 : where \f$Y_{4m}\f$ is one of the fourth order spherical harmonics so \f$m\f$ is a number that runs from \f$-4\f$ to
131 : \f$+4\f$. The function \f$\sigma( r_{ij} )\f$ is a \ref switchingfunction that acts on the distance between
132 : atoms \f$i\f$ and \f$j\f$. The parameters of this function should be set so that it the function is equal to one
133 : when atom \f$j\f$ is in the first coordination sphere of atom \f$i\f$ and is zero otherwise.
134 :
135 : The Steinhardt parameters can be used to measure the degree of order in the system in a variety of different ways. The
136 : simplest way of measuring whether or not the coordination sphere is ordered is to simply take the norm of the above vector i.e.
137 :
138 : \f[
139 : Q_4(i) = \sqrt{ \sum_{m=-4}^4 q_{4m}(i)^{*} q_{4m}(i) }
140 : \f]
141 :
142 : This norm is small when the coordination shell is disordered and larger when the coordination shell is ordered. Furthermore, when
143 : the keywords LESS_THAN, MIN, MAX, HISTOGRAM, MEAN and so on are used with this colvar it is the distribution of these normed quantities
144 : that is investigated.
145 :
146 : Other measures of order can be taken by averaging the components of the individual \f$q_4\f$ vectors individually or by taking dot products of
147 : the \f$q_{4}\f$ vectors on adjacent atoms. More information on these variables can be found in the documentation for \ref LOCAL_Q4,
148 : \ref LOCAL_AVERAGE and \ref NLINKS.
149 :
150 : \par Examples
151 :
152 : The following command calculates the average Q4 parameter for the 64 atoms in a box of Lennard Jones and prints this
153 : quantity to a file called colvar:
154 :
155 : \plumedfile
156 : Q4 SPECIES=1-64 D_0=1.3 R_0=0.2 MEAN LABEL=q4
157 : PRINT ARG=q4.mean FILE=colvar
158 : \endplumedfile
159 :
160 : The following command calculates the histogram of Q4 parameters for the 64 atoms in a box of Lennard Jones and prints these
161 : quantities to a file called colvar:
162 :
163 : \plumedfile
164 : 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} LABEL=q4
165 : PRINT ARG=q4.* FILE=colvar
166 : \endplumedfile
167 :
168 : The following command could be used to measure the Q4 parameters that describe the arrangement of chlorine ions around the
169 : sodium atoms in sodium chloride. The imagined system here is composed of 64 NaCl formula units and the atoms are arranged in the input
170 : with the 64 Na\f$^+\f$ ions followed by the 64 Cl\f$-\f$ ions. Once again the average Q4 parameter is calculated and output to a
171 : file called colvar
172 :
173 : \plumedfile
174 : Q4 SPECIESA=1-64 SPECIESB=65-128 D_0=1.3 R_0=0.2 MEAN LABEL=q4
175 : PRINT ARG=q4.mean FILE=colvar
176 : \endplumedfile
177 :
178 : 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
179 : command \ref DUMPATOMS as shown in the example below. The following output file will output a file in an extended xyz format
180 : called q$.xyz for each frame of the analyzed MD trajectory. The first column in this file will contain a dummy name for each of the
181 : 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
182 : 6-15 will contain the real parts of the director of the \f$q_{6m}\f$ vector while columns 15-24 will contain the imaginary parts of this director.
183 :
184 : \plumedfile
185 : q4: Q4 SPECIESA=1-64 SPECIESB=65-128 D_0=1.3 R_0=0.2 MEAN
186 : DUMPATOMS ATOMS=q4 ARG=q4_anorm FILE=q4.xyz
187 : \endplumedfile
188 :
189 : */
190 : //+ENDPLUMEDOC
191 :
192 : //+PLUMEDOC MCOLVAR Q6
193 : /*
194 : Calculate sixth order Steinhardt parameters.
195 :
196 : The sixth order Steinhardt parameters allow us to measure the degree to which the first coordination shell
197 : around an atom is ordered. The Steinhardt parameter for atom, \f$i\f$ is complex vector whose components are
198 : calculated using the following formula:
199 :
200 : \f[
201 : q_{6m}(i) = \frac{\sum_j \sigma( r_{ij} ) Y_{6m}(\mathbf{r}_{ij}) }{\sum_j \sigma( r_{ij} ) }
202 : \f]
203 :
204 : where \f$Y_{6m}\f$ is one of the sixth order spherical harmonics so \f$m\f$ is a number that runs from \f$-6\f$ to
205 : \f$+6\f$. The function \f$\sigma( r_{ij} )\f$ is a \ref switchingfunction that acts on the distance between
206 : atoms \f$i\f$ and \f$j\f$. The parameters of this function should be set so that it the function is equal to one
207 : when atom \f$j\f$ is in the first coordination sphere of atom \f$i\f$ and is zero otherwise.
208 :
209 : The Steinhardt parameters can be used to measure the degree of order in the system in a variety of different ways. The
210 : simplest way of measuring whether or not the coordination sphere is ordered is to simply take the norm of the above vector i.e.
211 :
212 : \f[
213 : Q_6(i) = \sqrt{ \sum_{m=-6}^6 q_{6m}(i)^{*} q_{6m}(i) }
214 : \f]
215 :
216 : This norm is small when the coordination shell is disordered and larger when the coordination shell is ordered. Furthermore, when
217 : the keywords LESS_THAN, MIN, MAX, HISTOGRAM, MEAN and so on are used with this colvar it is the distribution of these normed quantities
218 : that is investigated.
219 :
220 : Other measures of order can be taken by averaging the components of the individual \f$q_6\f$ vectors individually or by taking dot products of
221 : the \f$q_{6}\f$ vectors on adjacent atoms. More information on these variables can be found in the documentation for \ref LOCAL_Q6,
222 : \ref LOCAL_AVERAGE and \ref NLINKS.
223 :
224 : \par Examples
225 :
226 : The following command calculates the average Q6 parameter for the 64 atoms in a box of Lennard Jones and prints this
227 : quantity to a file called colvar:
228 :
229 : \plumedfile
230 : Q6 SPECIES=1-64 D_0=1.3 R_0=0.2 MEAN LABEL=q6
231 : PRINT ARG=q6.mean FILE=colvar
232 : \endplumedfile
233 :
234 : The following command calculates the histogram of Q6 parameters for the 64 atoms in a box of Lennard Jones and prints these
235 : quantities to a file called colvar:
236 :
237 : \plumedfile
238 : 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} LABEL=q6
239 : PRINT ARG=q6.* FILE=colvar
240 : \endplumedfile
241 :
242 : The following command could be used to measure the Q6 parameters that describe the arrangement of chlorine ions around the
243 : sodium atoms in sodium chloride. The imagined system here is composed of 64 NaCl formula units and the atoms are arranged in the input
244 : with the 64 Na\f$^+\f$ ions followed by the 64 Cl\f$-\f$ ions. Once again the average Q6 parameter is calculated and output to a
245 : file called colvar
246 :
247 : \plumedfile
248 : Q6 SPECIESA=1-64 SPECIESB=65-128 D_0=1.3 R_0=0.2 MEAN LABEL=q6
249 : PRINT ARG=q6.mean FILE=colvar
250 : \endplumedfile
251 :
252 : 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
253 : command \ref DUMPATOMS as shown in the example below. The following output file will output a file in an extended xyz format
254 : 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
255 : 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
256 : 6-19 will contain the real parts of the director of the \f$q_{6m}\f$ vector while columns 20-33 will contain the imaginary parts of this director.
257 :
258 : \plumedfile
259 : q6: Q6 SPECIESA=1-64 SPECIESB=65-128 D_0=1.3 R_0=0.2 MEAN
260 : DUMPATOMS ARG=q6_anorm ATOMS=q6 FILE=q6.xyz
261 : \endplumedfile
262 :
263 : */
264 : //+ENDPLUMEDOC
265 :
266 : class Steinhardt : public ActionShortcut {
267 : private:
268 : std::string getSymbol( const int& m ) const ;
269 : void createVectorNormInput( const std::string& ilab, const std::string& olab, const int& l, const std::string& sep, const std::string& vlab );
270 : public:
271 : static void registerKeywords( Keywords& keys );
272 : explicit Steinhardt(const ActionOptions&);
273 : };
274 :
275 : PLUMED_REGISTER_ACTION(Steinhardt,"Q1")
276 : PLUMED_REGISTER_ACTION(Steinhardt,"Q3")
277 : PLUMED_REGISTER_ACTION(Steinhardt,"Q4")
278 : PLUMED_REGISTER_ACTION(Steinhardt,"Q6")
279 :
280 60 : void Steinhardt::registerKeywords( Keywords& keys ) {
281 60 : CoordinationNumbers::shortcutKeywords( keys );
282 120 : keys.addFlag("LOWMEM",false,"this flag does nothing and is present only to ensure back-compatibility");
283 120 : keys.addFlag("VMEAN",false,"calculate the norm of the mean vector.");
284 120 : keys.addOutputComponent("_vmean","VMEAN","the norm of the mean vector");
285 120 : keys.addFlag("VSUM",false,"calculate the norm of the sum of all the vectors");
286 120 : keys.addOutputComponent("_vsum","VSUM","the norm of the mean vector");
287 60 : keys.needsAction("GROUP");
288 60 : keys.needsAction("CONTACT_MATRIX");
289 60 : keys.needsAction("SPHERICAL_HARMONIC");
290 60 : keys.needsAction("ONES");
291 60 : keys.needsAction("MATRIX_VECTOR_PRODUCT");
292 60 : keys.needsAction("COMBINE");
293 60 : keys.needsAction("CUSTOM");
294 60 : keys.needsAction("MEAN");
295 60 : keys.needsAction("SUM");
296 60 : }
297 :
298 42 : Steinhardt::Steinhardt( const ActionOptions& ao):
299 : Action(ao),
300 42 : ActionShortcut(ao) {
301 : bool lowmem;
302 42 : parseFlag("LOWMEM",lowmem);
303 42 : if( lowmem ) {
304 0 : warning("LOWMEM flag is deprecated and is no longer required for this action");
305 : }
306 : std::string sp_str, specA, specB;
307 42 : parse("SPECIES",sp_str);
308 42 : parse("SPECIESA",specA);
309 42 : parse("SPECIESB",specB);
310 42 : CoordinationNumbers::expandMatrix( true, getShortcutLabel(), sp_str, specA, specB, this );
311 : int l;
312 84 : std::string sph_input = getShortcutLabel() + "_sh: SPHERICAL_HARMONIC ARG=" + getShortcutLabel() + "_mat.x," + getShortcutLabel() + "_mat.y," + getShortcutLabel() + "_mat.z," + getShortcutLabel() + "_mat.w";
313 :
314 42 : if( getName()=="Q1" ) {
315 : sph_input +=" L=1";
316 19 : l=1;
317 23 : } else if( getName()=="Q3" ) {
318 : sph_input += " L=3";
319 1 : l=3;
320 22 : } else if( getName()=="Q4" ) {
321 : sph_input += " L=4";
322 3 : l=4;
323 19 : } else if( getName()=="Q6" ) {
324 : sph_input += " L=6";
325 19 : l=6;
326 : } else {
327 0 : plumed_merror("invalid input");
328 : }
329 42 : readInputLine( sph_input );
330 :
331 : // Input for denominator (coord)
332 42 : ActionWithValue* av = plumed.getActionSet().selectWithLabel<ActionWithValue*>( getShortcutLabel() + "_mat");
333 42 : plumed_assert( av && av->getNumberOfComponents()>0 && (av->copyOutput(0))->getRank()==2 );
334 : std::string size;
335 42 : Tools::convert( (av->copyOutput(0))->getShape()[1], size );
336 84 : readInputLine( getShortcutLabel() + "_denom_ones: ONES SIZE=" + size );
337 84 : readInputLine( getShortcutLabel() + "_denom: MATRIX_VECTOR_PRODUCT ARG=" + getShortcutLabel() + "_mat.w," + getShortcutLabel() + "_denom_ones" );
338 84 : readInputLine( getShortcutLabel() + "_sp: MATRIX_VECTOR_PRODUCT ARG=" + getShortcutLabel() + "_sh.*," + getShortcutLabel() + "_denom_ones");
339 :
340 : // If we are doing VMEAN determine sum of vector components
341 : std::string snum;
342 : bool do_vmean;
343 42 : parseFlag("VMEAN",do_vmean);
344 : bool do_vsum;
345 42 : parseFlag("VSUM",do_vsum);
346 42 : if( do_vmean || do_vsum ) {
347 : // Divide all components by coordination numbers
348 14 : for(int i=-l; i<=l; ++i) {
349 13 : snum = getSymbol( i );
350 : // Real part
351 26 : readInputLine( getShortcutLabel() + "_rmn-" + snum + ": CUSTOM ARG=" + getShortcutLabel() + "_sp.rm-" + snum + "," + getShortcutLabel() + "_denom FUNC=x/y PERIODIC=NO");
352 : // Imaginary part
353 26 : readInputLine( getShortcutLabel() + "_imn-" + snum + ": CUSTOM ARG=" + getShortcutLabel() + "_sp.im-" + snum + "," + getShortcutLabel() + "_denom FUNC=x/y PERIODIC=NO");
354 : }
355 : }
356 :
357 42 : if( do_vmean ) {
358 14 : for(int i=-l; i<=l; ++i) {
359 13 : snum = getSymbol( i );
360 : // Real part
361 26 : readInputLine( getShortcutLabel() + "_rms-" + snum + ": MEAN ARG=" + getShortcutLabel() + "_rmn-" + snum + " PERIODIC=NO");
362 : // Imaginary part
363 26 : readInputLine( getShortcutLabel() + "_ims-" + snum + ": MEAN ARG=" + getShortcutLabel() + "_imn-" + snum + " PERIODIC=NO");
364 : }
365 : // Now calculate the total length of the vector
366 2 : createVectorNormInput( getShortcutLabel(), getShortcutLabel() + "_vmean", l, "_", "ms" );
367 : }
368 42 : if( do_vsum ) {
369 0 : for(int i=-l; i<=l; ++i) {
370 0 : snum = getSymbol( i );
371 : // Real part
372 0 : readInputLine( getShortcutLabel() + "_rmz-" + snum + ": SUM ARG=" + getShortcutLabel() + "_rmn-" + snum + " PERIODIC=NO");
373 : // Imaginary part
374 0 : readInputLine( getShortcutLabel() + "_imz-" + snum + ": SUM ARG=" + getShortcutLabel() + "_imn-" + snum + " PERIODIC=NO");
375 : }
376 : // Now calculate the total length of the vector
377 0 : createVectorNormInput( getShortcutLabel(), getShortcutLabel() + "_vsum", l, "_", "mz" );
378 : }
379 :
380 : // Now calculate the total length of the vector
381 84 : createVectorNormInput( getShortcutLabel() + "_sp", getShortcutLabel() + "_norm", l, ".", "m" );
382 : // And take average
383 84 : readInputLine( getShortcutLabel() + ": CUSTOM ARG=" + getShortcutLabel() + "_norm," + getShortcutLabel() + "_denom FUNC=x/y PERIODIC=NO");
384 84 : multicolvar::MultiColvarShortcuts::expandFunctions( getShortcutLabel(), getShortcutLabel(), "", this );
385 42 : }
386 :
387 43 : void Steinhardt::createVectorNormInput( const std::string& ilab, const std::string& olab, const int& l, const std::string& sep, const std::string& vlab ) {
388 43 : std::string arg_inp, norm_input = olab + "2: COMBINE PERIODIC=NO POWERS=2,2";
389 43 : std::string snum = getSymbol( -l );
390 86 : arg_inp = " ARG=" + ilab + sep + "r" + vlab + "-" + snum +"," + ilab + sep + "i" + vlab + "-" + snum;
391 351 : for(int i=-l+1; i<=l; ++i) {
392 308 : snum = getSymbol( i );
393 616 : arg_inp += "," + ilab + sep + "r" + vlab + "-" + snum + "," + ilab + sep + "i" + vlab + "-" + snum;
394 : norm_input += ",2,2";
395 : }
396 43 : readInputLine( norm_input + arg_inp );
397 86 : readInputLine( olab + ": CUSTOM ARG=" + olab + "2 FUNC=sqrt(x) PERIODIC=NO");
398 43 : }
399 :
400 377 : std::string Steinhardt::getSymbol( const int& m ) const {
401 377 : if( m<0 ) {
402 : std::string num;
403 166 : Tools::convert( -1*m, num );
404 166 : return "n" + num;
405 211 : } else if( m>0 ) {
406 : std::string num;
407 166 : Tools::convert( m, num );
408 166 : return "p" + num;
409 : }
410 45 : return "0";
411 : }
412 :
413 : }
414 : }
415 :
|
