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 240 : keys.needsAction("GROUP"); keys.needsAction("CONTACT_MATRIX"); keys.needsAction("SPHERICAL_HARMONIC"); keys.needsAction("ONES");
288 300 : keys.needsAction("MATRIX_VECTOR_PRODUCT"); keys.needsAction("COMBINE"); keys.needsAction("CUSTOM"); keys.needsAction("MEAN"); keys.needsAction("SUM");
289 60 : }
290 :
291 42 : Steinhardt::Steinhardt( const ActionOptions& ao):
292 : Action(ao),
293 42 : ActionShortcut(ao)
294 : {
295 42 : bool lowmem; parseFlag("LOWMEM",lowmem);
296 42 : if( lowmem ) warning("LOWMEM flag is deprecated and is no longer required for this action");
297 126 : std::string sp_str, specA, specB; parse("SPECIES",sp_str); parse("SPECIESA",specA); parse("SPECIESB",specB);
298 42 : CoordinationNumbers::expandMatrix( true, getShortcutLabel(), sp_str, specA, specB, this ); int l;
299 84 : std::string sph_input = getShortcutLabel() + "_sh: SPHERICAL_HARMONIC ARG=" + getShortcutLabel() + "_mat.x," + getShortcutLabel() + "_mat.y," + getShortcutLabel() + "_mat.z," + getShortcutLabel() + "_mat.w";
300 :
301 42 : if( getName()=="Q1" ) {
302 19 : sph_input +=" L=1"; l=1;
303 23 : } else if( getName()=="Q3" ) {
304 1 : sph_input += " L=3"; l=3;
305 22 : } else if( getName()=="Q4" ) {
306 3 : sph_input += " L=4"; l=4;
307 19 : } else if( getName()=="Q6" ) {
308 19 : sph_input += " L=6"; l=6;
309 : } else {
310 0 : plumed_merror("invalid input");
311 : }
312 42 : readInputLine( sph_input );
313 :
314 : // Input for denominator (coord)
315 42 : ActionWithValue* av = plumed.getActionSet().selectWithLabel<ActionWithValue*>( getShortcutLabel() + "_mat");
316 42 : plumed_assert( av && av->getNumberOfComponents()>0 && (av->copyOutput(0))->getRank()==2 );
317 42 : std::string size; Tools::convert( (av->copyOutput(0))->getShape()[1], size );
318 84 : readInputLine( getShortcutLabel() + "_denom_ones: ONES SIZE=" + size );
319 84 : readInputLine( getShortcutLabel() + "_denom: MATRIX_VECTOR_PRODUCT ARG=" + getShortcutLabel() + "_mat.w," + getShortcutLabel() + "_denom_ones" );
320 84 : readInputLine( getShortcutLabel() + "_sp: MATRIX_VECTOR_PRODUCT ARG=" + getShortcutLabel() + "_sh.*," + getShortcutLabel() + "_denom_ones");
321 :
322 : // If we are doing VMEAN determine sum of vector components
323 : std::string snum;
324 42 : bool do_vmean; parseFlag("VMEAN",do_vmean);
325 42 : bool do_vsum; parseFlag("VSUM",do_vsum);
326 42 : if( do_vmean || do_vsum ) {
327 : // Divide all components by coordination numbers
328 14 : for(int i=-l; i<=l; ++i) {
329 13 : snum = getSymbol( i );
330 : // Real part
331 26 : readInputLine( getShortcutLabel() + "_rmn-" + snum + ": CUSTOM ARG=" + getShortcutLabel() + "_sp.rm-" + snum + "," + getShortcutLabel() + "_denom FUNC=x/y PERIODIC=NO");
332 : // Imaginary part
333 26 : readInputLine( getShortcutLabel() + "_imn-" + snum + ": CUSTOM ARG=" + getShortcutLabel() + "_sp.im-" + snum + "," + getShortcutLabel() + "_denom FUNC=x/y PERIODIC=NO");
334 : }
335 : }
336 :
337 42 : if( do_vmean ) {
338 14 : for(int i=-l; i<=l; ++i) {
339 13 : snum = getSymbol( i );
340 : // Real part
341 26 : readInputLine( getShortcutLabel() + "_rms-" + snum + ": MEAN ARG=" + getShortcutLabel() + "_rmn-" + snum + " PERIODIC=NO");
342 : // Imaginary part
343 26 : readInputLine( getShortcutLabel() + "_ims-" + snum + ": MEAN ARG=" + getShortcutLabel() + "_imn-" + snum + " PERIODIC=NO");
344 : }
345 : // Now calculate the total length of the vector
346 2 : createVectorNormInput( getShortcutLabel(), getShortcutLabel() + "_vmean", l, "_", "ms" );
347 : }
348 42 : if( do_vsum ) {
349 0 : for(int i=-l; i<=l; ++i) {
350 0 : snum = getSymbol( i );
351 : // Real part
352 0 : readInputLine( getShortcutLabel() + "_rmz-" + snum + ": SUM ARG=" + getShortcutLabel() + "_rmn-" + snum + " PERIODIC=NO");
353 : // Imaginary part
354 0 : readInputLine( getShortcutLabel() + "_imz-" + snum + ": SUM ARG=" + getShortcutLabel() + "_imn-" + snum + " PERIODIC=NO");
355 : }
356 : // Now calculate the total length of the vector
357 0 : createVectorNormInput( getShortcutLabel(), getShortcutLabel() + "_vsum", l, "_", "mz" );
358 : }
359 :
360 : // Now calculate the total length of the vector
361 84 : createVectorNormInput( getShortcutLabel() + "_sp", getShortcutLabel() + "_norm", l, ".", "m" );
362 : // And take average
363 84 : readInputLine( getShortcutLabel() + ": CUSTOM ARG=" + getShortcutLabel() + "_norm," + getShortcutLabel() + "_denom FUNC=x/y PERIODIC=NO");
364 84 : multicolvar::MultiColvarShortcuts::expandFunctions( getShortcutLabel(), getShortcutLabel(), "", this );
365 42 : }
366 :
367 43 : void Steinhardt::createVectorNormInput( const std::string& ilab, const std::string& olab, const int& l, const std::string& sep, const std::string& vlab ) {
368 43 : std::string arg_inp, norm_input = olab + "2: COMBINE PERIODIC=NO POWERS=2,2"; std::string snum = getSymbol( -l );
369 86 : arg_inp = " ARG=" + ilab + sep + "r" + vlab + "-" + snum +"," + ilab + sep + "i" + vlab + "-" + snum;
370 351 : for(int i=-l+1; i<=l; ++i) {
371 308 : snum = getSymbol( i );
372 616 : arg_inp += "," + ilab + sep + "r" + vlab + "-" + snum + "," + ilab + sep + "i" + vlab + "-" + snum;
373 : norm_input += ",2,2";
374 : }
375 43 : readInputLine( norm_input + arg_inp );
376 86 : readInputLine( olab + ": CUSTOM ARG=" + olab + "2 FUNC=sqrt(x) PERIODIC=NO");
377 43 : }
378 :
379 377 : std::string Steinhardt::getSymbol( const int& m ) const {
380 377 : if( m<0 ) {
381 166 : std::string num; Tools::convert( -1*m, num );
382 166 : return "n" + num;
383 211 : } else if( m>0 ) {
384 166 : std::string num; Tools::convert( m, num );
385 166 : return "p" + num;
386 : }
387 45 : return "0";
388 : }
389 :
390 : }
391 : }
392 :
|
