Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2012-2017 The plumed team
3 : (see the PEOPLE file at the root of the distribution for a list of names)
4 :
5 : See http://www.plumed.org for more information.
6 :
7 : This file is part of plumed, version 2.
8 :
9 : plumed is free software: you can redistribute it and/or modify
10 : it under the terms of the GNU Lesser General Public License as published by
11 : the Free Software Foundation, either version 3 of the License, or
12 : (at your option) any later version.
13 :
14 : plumed is distributed in the hope that it will be useful,
15 : but WITHOUT ANY WARRANTY; without even the implied warranty of
16 : MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
17 : GNU Lesser General Public License for more details.
18 :
19 : You should have received a copy of the GNU Lesser General Public License
20 : along with plumed. If not, see <http://www.gnu.org/licenses/>.
21 : +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
22 : #include "CoordinationNumbers.h"
23 : #include "core/ActionShortcut.h"
24 : #include "multicolvar/MultiColvarShortcuts.h"
25 : #include "core/ActionWithValue.h"
26 : #include "core/PlumedMain.h"
27 : #include "core/ActionSet.h"
28 : #include "core/ActionRegister.h"
29 :
30 : #include <complex>
31 :
32 : namespace PLMD {
33 : namespace symfunc {
34 :
35 : //+PLUMEDOC MCOLVAR Q1
36 : /*
37 : Calculate 1st order Steinhardt parameters
38 :
39 : The 1st order Steinhardt parameters allow us to measure the degree to which the first coordination shell
40 : around an atom is ordered with the atoms aranged on a line. The Steinhardt parameter for atom, $i$ is complex vector whose components are
41 : calculated using the following formula:
42 :
43 : $$
44 : q_{1m}(i) = \frac{\sum_j \sigma( r_{ij} ) Y_{1m}(\mathbf{r}_{ij}) }{\sum_j \sigma( r_{ij} ) }
45 : $$
46 :
47 : where $Y_{1m}$ is one of the 1st order spherical harmonics so $m$ is a number that runs from $-1$ to
48 : $+1$. The function $\sigma( r_{ij} )$ is a switching function that acts on the distance between
49 : atoms $i$ and $j$. The parameters of this function should be set so that it the function is equal to one
50 : when atom $j$ is in the first coordination sphere of atom $i$ and is zero otherwise.
51 :
52 : As discussed on [this page](https://www.plumed-tutorials.org/lessons/23/001/data/Steinhardt.html), the Steinhardt parameters can
53 : be used to measure the degree of order in the system in a variety of different ways. The
54 : simplest way of measuring whether or not the coordination sphere is ordered is to simply take the norm of the above vector i.e.
55 :
56 : $$
57 : Q_1(i) = \sqrt{ \sum_{m=-1}^1 q_{1m}(i)^{*} q_{1m}(i) }
58 : $$
59 :
60 : This norm is small when the coordination shell is disordered and larger when the coordination shell is ordered. Furthermore, in inputs like
61 : the one shown below where the keywords LESS_THAN, MIN, MAX, HISTOGRAM, MEAN and so on are used with it is the distribution of these normed quantities
62 : that is investigated. You can investigate precisely what calculation is performed here by expaning the shortcuts in the input below.
63 :
64 : ```plumed
65 : q1: Q1 SPECIES=1-64 D_0=1.3 R_0=0.2 MEAN
66 : PRINT ARG=q1.mean FILE=colvar
67 : ```
68 :
69 : Another similar command is provided below. This time the histogram of Q1 parameters for the 64 atoms in a box of Lennard Jones is computed and printed to
70 : to a file called colvar:
71 :
72 : ```plumed
73 : q1: Q1 SPECIES=1-64 D_0=1.3 R_0=0.2 HISTOGRAM={GAUSSIAN LOWER=0.0 UPPER=1.0 NBINS=20 SMEAR=0.1}
74 : PRINT ARG=q1.* FILE=colvar
75 : ```
76 :
77 : The command below could be used to measure the Q1 parameters that describe the arrangement of chlorine ions around the
78 : sodium atoms in sodium chloride. The imagined system here is composed of 64 NaCl formula units and the atoms are arranged in the input
79 : with the 64 Na$^+$ ions followed by the 64 Cl$-$ ions. Once again the average Q1 parameter is calculated and output to a
80 : file called colvar
81 :
82 : ```plumed
83 : q1: Q1 SPECIESA=1-64 SPECIESB=65-128 D_0=1.3 R_0=0.2 MEAN
84 : PRINT ARG=q1.mean FILE=colvar
85 : ```
86 :
87 : If you simply want to examine the values of the Q1 parameters for each of the atoms in your system you can do so by exploiting the
88 : command [DUMPATOMS](DUMPATOMS.md) as shown in the example below. The following output file will output a file in an extended xyz format
89 : called q1.xyz for each frame of the analyzed MD trajectory. The first column in this file will contain a dummy name for each of the
90 : atoms, columns 2-4 will then contain the x, y and z positions of the atoms, column 5 will contain the value of the Q1 parameter, columns
91 : 6-12 will contain the real parts of the director of the $q_{1m}$ vector while columns 12-19 will contain the imaginary parts of this director.
92 :
93 : ```plumed
94 : q1: Q1 SPECIESA=1-64 SPECIESB=65-128 D_0=1.3 R_0=0.2 MEAN
95 : DUMPATOMS ATOMS=q1 ARG=q1 FILE=q1.xyz
96 : ```
97 :
98 : */
99 : //+ENDPLUMEDOC
100 :
101 : //+PLUMEDOC MCOLVAR Q3
102 : /*
103 : Calculate 3rd order Steinhardt parameters.
104 :
105 : The 3rd order Steinhardt parameters allow us to measure the degree to which the first coordination shell
106 : around an atom is ordered. The Steinhardt parameter for atom, $i$ is complex vector whose components are
107 : calculated using the following formula:
108 :
109 : $$
110 : q_{3m}(i) = \frac{\sum_j \sigma( r_{ij} ) Y_{3m}(\mathbf{r}_{ij}) }{\sum_j \sigma( r_{ij} ) }
111 : $$
112 :
113 : where $Y_{3m}$ is one of the 3rd order spherical harmonics so $m$ is a number that runs from $-3$ to
114 : $+3$. The function $\sigma( r_{ij} )$ is a switching function that acts on the distance between
115 : atoms $i$ and $j$. The parameters of this function should be set so that it the function is equal to one
116 : when atom $j$ is in the first coordination sphere of atom $i$ and is zero otherwise.
117 :
118 : As discussed on [this page](https://www.plumed-tutorials.org/lessons/23/001/data/Steinhardt.html), the Steinhardt parameters can
119 : be used to measure the degree of order in the system in a variety of different ways. The
120 : simplest way of measuring whether or not the coordination sphere is ordered is to simply take the norm of the above vector i.e.
121 :
122 : $$
123 : Q_3(i) = \sqrt{ \sum_{m=-3}^3 q_{3m}(i)^{*} q_{3m}(i) }
124 : $$
125 :
126 : This norm is small when the coordination shell is disordered and larger when the coordination shell is ordered. Furthermore, in inputs like
127 : the one shown below where the keywords LESS_THAN, MIN, MAX, HISTOGRAM, MEAN and so on are used with it is the distribution of these normed quantities
128 : that is investigated. You can investigate precisely what calculation is performed here by expaning the shortcuts in the input below.
129 :
130 : ```plumed
131 : q3: Q3 SPECIES=1-64 D_0=1.3 R_0=0.2 MEAN
132 : PRINT ARG=q3.mean FILE=colvar
133 : ```
134 :
135 : Another similar command is provided below. This time the histogram of Q3 parameters for the 64 atoms in a box of Lennard Jones is computed and printed to
136 : to a file called colvar:
137 :
138 : ```plumed
139 : q3: Q3 SPECIES=1-64 D_0=1.3 R_0=0.2 HISTOGRAM={GAUSSIAN LOWER=0.0 UPPER=1.0 NBINS=20 SMEAR=0.1}
140 : PRINT ARG=q3.* FILE=colvar
141 : ```
142 :
143 : The command below could be used to measure the Q3 parameters that describe the arrangement of chlorine ions around the
144 : sodium atoms in sodium chloride. The imagined system here is composed of 64 NaCl formula units and the atoms are arranged in the input
145 : with the 64 Na$^+$ ions followed by the 64 Cl$-$ ions. Once again the average Q3 parameter is calculated and output to a
146 : file called colvar
147 :
148 : ```plumed
149 : q3: Q3 SPECIESA=1-64 SPECIESB=65-128 D_0=1.3 R_0=0.2 MEAN
150 : PRINT ARG=q3.mean FILE=colvar
151 : ```
152 :
153 : 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
154 : command [DUMPATOMS](DUMPATOMS.md) as shown in the example below. The following output file will output a file in an extended xyz format
155 : 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
156 : 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
157 : 6-12 will contain the real parts of the director of the $q_{3m}$ vector while columns 12-19 will contain the imaginary parts of this director.
158 :
159 : ```plumed
160 : q3: Q3 SPECIESA=1-64 SPECIESB=65-128 D_0=1.3 R_0=0.2 MEAN
161 : DUMPATOMS ATOMS=q3 ARG=q3 FILE=q3.xyz
162 : ```
163 :
164 : */
165 : //+ENDPLUMEDOC
166 :
167 : //+PLUMEDOC MCOLVAR Q4
168 : /*
169 : Calculate fourth order Steinhardt parameters.
170 :
171 : The 4th order Steinhardt parameters allow us to measure the degree to which the first coordination shell
172 : around an atom is ordered. The Steinhardt parameter for atom, $i$ is complex vector whose components are
173 : calculated using the following formula:
174 :
175 : $$
176 : q_{4m}(i) = \frac{\sum_j \sigma( r_{ij} ) Y_{4m}(\mathbf{r}_{ij}) }{\sum_j \sigma( r_{ij} ) }
177 : $$
178 :
179 : where $Y_{4m}$ is one of the 4th order spherical harmonics so $m$ is a number that runs from $-4$ to
180 : $+4$. The function $\sigma( r_{ij} )$ is a switching function that acts on the distance between
181 : atoms $i$ and $j$. The parameters of this function should be set so that it the function is equal to one
182 : when atom $j$ is in the first coordination sphere of atom $i$ and is zero otherwise.
183 :
184 : As discussed on [this page](https://www.plumed-tutorials.org/lessons/23/001/data/Steinhardt.html), the Steinhardt parameters can
185 : be used to measure the degree of order in the system in a variety of different ways. The
186 : simplest way of measuring whether or not the coordination sphere is ordered is to simply take the norm of the above vector i.e.
187 :
188 : $$
189 : Q_4(i) = \sqrt{ \sum_{m=-4}^4 q_{4m}(i)^{*} q_{4m}(i) }
190 : $$
191 :
192 : This norm is small when the coordination shell is disordered and larger when the coordination shell is ordered. Furthermore, in inputs like
193 : the one shown below where the keywords LESS_THAN, MIN, MAX, HISTOGRAM, MEAN and so on are used with it is the distribution of these normed quantities
194 : that is investigated. You can investigate precisely what calculation is performed here by expaning the shortcuts in the input below.
195 :
196 : ```plumed
197 : q4: Q4 SPECIES=1-64 D_0=1.3 R_0=0.2 MEAN
198 : PRINT ARG=q4.mean FILE=colvar
199 : ```
200 :
201 : Another similar command is provided below. This time the histogram of Q4 parameters for the 64 atoms in a box of Lennard Jones is computed and printed to
202 : to a file called colvar:
203 :
204 : ```plumed
205 : q4: Q4 SPECIES=1-64 D_0=1.3 R_0=0.2 HISTOGRAM={GAUSSIAN LOWER=0.0 UPPER=1.0 NBINS=20 SMEAR=0.1}
206 : PRINT ARG=q4.* FILE=colvar
207 : ```
208 :
209 : The command below could be used to measure the Q4 parameters that describe the arrangement of chlorine ions around the
210 : sodium atoms in sodium chloride. The imagined system here is composed of 64 NaCl formula units and the atoms are arranged in the input
211 : with the 64 Na$^+$ ions followed by the 64 Cl$-$ ions. Once again the average Q4 parameter is calculated and output to a
212 : file called colvar
213 :
214 : ```plumed
215 : q4: Q4 SPECIESA=1-64 SPECIESB=65-128 D_0=1.3 R_0=0.2 MEAN
216 : PRINT ARG=q4.mean FILE=colvar
217 : ```
218 :
219 : 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
220 : command [DUMPATOMS](DUMPATOMS.md) as shown in the example below. The following output file will output a file in an extended xyz format
221 : called q4.xyz for each frame of the analyzed MD trajectory. The first column in this file will contain a dummy name for each of the
222 : 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
223 : 6-12 will contain the real parts of the director of the $q_{4m}$ vector while columns 12-19 will contain the imaginary parts of this director.
224 :
225 : ```plumed
226 : q4: Q4 SPECIESA=1-64 SPECIESB=65-128 D_0=1.3 R_0=0.2 MEAN
227 : DUMPATOMS ATOMS=q4 ARG=q4 FILE=q4.xyz
228 : ```
229 :
230 : */
231 : //+ENDPLUMEDOC
232 :
233 : //+PLUMEDOC MCOLVAR Q6
234 : /*
235 : Calculate sixth order Steinhardt parameters.
236 :
237 : The 6th order Steinhardt parameters allow us to measure the degree to which the first coordination shell
238 : around an atom is ordered. The Steinhardt parameter for atom, $i$ is complex vector whose components are
239 : calculated using the following formula:
240 :
241 : $$
242 : q_{6m}(i) = \frac{\sum_j \sigma( r_{ij} ) Y_{6m}(\mathbf{r}_{ij}) }{\sum_j \sigma( r_{ij} ) }
243 : $$
244 :
245 : where $Y_{6m}$ is one of the 6th order spherical harmonics so $m$ is a number that runs from $-6$ to
246 : $+6$. The function $\sigma( r_{ij} )$ is a switching function that acts on the distance between
247 : atoms $i$ and $j$. The parameters of this function should be set so that it the function is equal to one
248 : when atom $j$ is in the first coordination sphere of atom $i$ and is zero otherwise.
249 :
250 : As discussed on [this page](https://www.plumed-tutorials.org/lessons/23/001/data/Steinhardt.html), the Steinhardt parameters can
251 : be used to measure the degree of order in the system in a variety of different ways. The
252 : simplest way of measuring whether or not the coordination sphere is ordered is to simply take the norm of the above vector i.e.
253 :
254 : $$
255 : Q_6(i) = \sqrt{ \sum_{m=-6}^6 q_{6m}(i)^{*} q_{6m}(i) }
256 : $$
257 :
258 : This norm is small when the coordination shell is disordered and larger when the coordination shell is ordered. Furthermore, in inputs like
259 : the one shown below where the keywords LESS_THAN, MIN, MAX, HISTOGRAM, MEAN and so on are used with it is the distribution of these normed quantities
260 : that is investigated. You can investigate precisely what calculation is performed here by expaning the shortcuts in the input below.
261 :
262 : ```plumed
263 : q6: Q6 SPECIES=1-64 D_0=1.3 R_0=0.2 MEAN
264 : PRINT ARG=q6.mean FILE=colvar
265 : ```
266 :
267 : Another similar command is provided below. This time the histogram of Q6 parameters for the 64 atoms in a box of Lennard Jones is computed and printed to
268 : to a file called colvar:
269 :
270 : ```plumed
271 : q6: Q6 SPECIES=1-64 D_0=1.3 R_0=0.2 HISTOGRAM={GAUSSIAN LOWER=0.0 UPPER=1.0 NBINS=20 SMEAR=0.1}
272 : PRINT ARG=q6.* FILE=colvar
273 : ```
274 :
275 : The command below could be used to measure the Q6 parameters that describe the arrangement of chlorine ions around the
276 : sodium atoms in sodium chloride. The imagined system here is composed of 64 NaCl formula units and the atoms are arranged in the input
277 : with the 64 Na$^+$ ions followed by the 64 Cl$-$ ions. Once again the average Q6 parameter is calculated and output to a
278 : file called colvar
279 :
280 : ```plumed
281 : q6: Q6 SPECIESA=1-64 SPECIESB=65-128 D_0=1.3 R_0=0.2 MEAN
282 : PRINT ARG=q6.mean FILE=colvar
283 : ```
284 :
285 : 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
286 : command [DUMPATOMS](DUMPATOMS.md) as shown in the example below. The following output file will output a file in an extended xyz format
287 : 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
288 : 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
289 : 6-12 will contain the real parts of the director of the $q_{6m}$ vector while columns 12-19 will contain the imaginary parts of this director.
290 :
291 : ```plumed
292 : q6: Q6 SPECIESA=1-64 SPECIESB=65-128 D_0=1.3 R_0=0.2 MEAN
293 : DUMPATOMS ATOMS=q6 ARG=q6 FILE=q6.xyz
294 : ```
295 :
296 : */
297 : //+ENDPLUMEDOC
298 :
299 : class Steinhardt : public ActionShortcut {
300 : private:
301 : std::string getSymbol( const int& m ) const ;
302 : void createVectorNormInput( const std::string& ilab, const std::string& olab, const int& l, const std::string& sep, const std::string& vlab );
303 : public:
304 : static void registerKeywords( Keywords& keys );
305 : explicit Steinhardt(const ActionOptions&);
306 : };
307 :
308 : PLUMED_REGISTER_ACTION(Steinhardt,"Q1")
309 : PLUMED_REGISTER_ACTION(Steinhardt,"Q3")
310 : PLUMED_REGISTER_ACTION(Steinhardt,"Q4")
311 : PLUMED_REGISTER_ACTION(Steinhardt,"Q6")
312 :
313 60 : void Steinhardt::registerKeywords( Keywords& keys ) {
314 60 : CoordinationNumbers::shortcutKeywords( keys );
315 60 : keys.addFlag("LOWMEM",false,"this flag does nothing and is present only to ensure back-compatibility");
316 60 : keys.addFlag("VMEAN",false,"calculate the norm of the mean vector.");
317 120 : keys.addOutputComponent("_vmean","VMEAN","scalar","the norm of the mean vector");
318 60 : keys.addFlag("VSUM",false,"calculate the norm of the sum of all the vectors");
319 120 : keys.addOutputComponent("_vsum","VSUM","scalar","the norm of the mean vector");
320 60 : keys.needsAction("GROUP");
321 60 : keys.needsAction("CONTACT_MATRIX");
322 60 : keys.needsAction("SPHERICAL_HARMONIC");
323 60 : keys.needsAction("ONES");
324 60 : keys.needsAction("MATRIX_VECTOR_PRODUCT");
325 60 : keys.needsAction("COMBINE");
326 60 : keys.needsAction("CUSTOM");
327 60 : keys.needsAction("MEAN");
328 60 : keys.needsAction("SUM");
329 120 : keys.setValueDescription("vector","the norms of the vectors of spherical harmonic coefficients");
330 60 : }
331 :
332 42 : Steinhardt::Steinhardt( const ActionOptions& ao):
333 : Action(ao),
334 42 : ActionShortcut(ao) {
335 : bool lowmem;
336 42 : parseFlag("LOWMEM",lowmem);
337 42 : if( lowmem ) {
338 0 : warning("LOWMEM flag is deprecated and is no longer required for this action");
339 : }
340 : std::string sp_str, specA, specB;
341 42 : parse("SPECIES",sp_str);
342 42 : parse("SPECIESA",specA);
343 42 : parse("SPECIESB",specB);
344 42 : CoordinationNumbers::expandMatrix( true, getShortcutLabel(), sp_str, specA, specB, this );
345 : int l;
346 84 : std::string sph_input = getShortcutLabel() + "_sh: SPHERICAL_HARMONIC ARG=" + getShortcutLabel() + "_mat.x," + getShortcutLabel() + "_mat.y," + getShortcutLabel() + "_mat.z," + getShortcutLabel() + "_mat.w";
347 :
348 42 : if( getName()=="Q1" ) {
349 : sph_input +=" L=1";
350 19 : l=1;
351 23 : } else if( getName()=="Q3" ) {
352 : sph_input += " L=3";
353 1 : l=3;
354 22 : } else if( getName()=="Q4" ) {
355 : sph_input += " L=4";
356 3 : l=4;
357 19 : } else if( getName()=="Q6" ) {
358 : sph_input += " L=6";
359 19 : l=6;
360 : } else {
361 0 : plumed_merror("invalid input");
362 : }
363 42 : readInputLine( sph_input );
364 :
365 : // Input for denominator (coord)
366 42 : ActionWithValue* av = plumed.getActionSet().selectWithLabel<ActionWithValue*>( getShortcutLabel() + "_mat");
367 42 : plumed_assert( av && av->getNumberOfComponents()>0 && (av->copyOutput(0))->getRank()==2 );
368 : std::string size;
369 42 : Tools::convert( (av->copyOutput(0))->getShape()[1], size );
370 84 : readInputLine( getShortcutLabel() + "_denom_ones: ONES SIZE=" + size );
371 84 : readInputLine( getShortcutLabel() + "_denom: MATRIX_VECTOR_PRODUCT ARG=" + getShortcutLabel() + "_mat.w," + getShortcutLabel() + "_denom_ones" );
372 84 : readInputLine( getShortcutLabel() + "_sp: MATRIX_VECTOR_PRODUCT ARG=" + getShortcutLabel() + "_sh.*," + getShortcutLabel() + "_denom_ones");
373 :
374 : // If we are doing VMEAN determine sum of vector components
375 : std::string snum;
376 : bool do_vmean;
377 42 : parseFlag("VMEAN",do_vmean);
378 : bool do_vsum;
379 42 : parseFlag("VSUM",do_vsum);
380 42 : if( do_vmean || do_vsum ) {
381 : // Divide all components by coordination numbers
382 14 : for(int i=-l; i<=l; ++i) {
383 13 : snum = getSymbol( i );
384 : // Real part
385 26 : readInputLine( getShortcutLabel() + "_rmn-" + snum + ": CUSTOM ARG=" + getShortcutLabel() + "_sp.rm-" + snum + "," + getShortcutLabel() + "_denom FUNC=x/y PERIODIC=NO");
386 : // Imaginary part
387 26 : readInputLine( getShortcutLabel() + "_imn-" + snum + ": CUSTOM ARG=" + getShortcutLabel() + "_sp.im-" + snum + "," + getShortcutLabel() + "_denom FUNC=x/y PERIODIC=NO");
388 : }
389 : }
390 :
391 42 : if( do_vmean ) {
392 14 : for(int i=-l; i<=l; ++i) {
393 13 : snum = getSymbol( i );
394 : // Real part
395 26 : readInputLine( getShortcutLabel() + "_rms-" + snum + ": MEAN ARG=" + getShortcutLabel() + "_rmn-" + snum + " PERIODIC=NO");
396 : // Imaginary part
397 26 : readInputLine( getShortcutLabel() + "_ims-" + snum + ": MEAN ARG=" + getShortcutLabel() + "_imn-" + snum + " PERIODIC=NO");
398 : }
399 : // Now calculate the total length of the vector
400 2 : createVectorNormInput( getShortcutLabel(), getShortcutLabel() + "_vmean", l, "_", "ms" );
401 : }
402 42 : if( do_vsum ) {
403 0 : for(int i=-l; i<=l; ++i) {
404 0 : snum = getSymbol( i );
405 : // Real part
406 0 : readInputLine( getShortcutLabel() + "_rmz-" + snum + ": SUM ARG=" + getShortcutLabel() + "_rmn-" + snum + " PERIODIC=NO");
407 : // Imaginary part
408 0 : readInputLine( getShortcutLabel() + "_imz-" + snum + ": SUM ARG=" + getShortcutLabel() + "_imn-" + snum + " PERIODIC=NO");
409 : }
410 : // Now calculate the total length of the vector
411 0 : createVectorNormInput( getShortcutLabel(), getShortcutLabel() + "_vsum", l, "_", "mz" );
412 : }
413 :
414 : // Now calculate the total length of the vector
415 84 : createVectorNormInput( getShortcutLabel() + "_sp", getShortcutLabel() + "_norm", l, ".", "m" );
416 : // And take average
417 84 : readInputLine( getShortcutLabel() + ": CUSTOM ARG=" + getShortcutLabel() + "_norm," + getShortcutLabel() + "_denom FUNC=x/y PERIODIC=NO");
418 84 : multicolvar::MultiColvarShortcuts::expandFunctions( getShortcutLabel(), getShortcutLabel(), "", this );
419 42 : }
420 :
421 43 : void Steinhardt::createVectorNormInput( const std::string& ilab, const std::string& olab, const int& l, const std::string& sep, const std::string& vlab ) {
422 43 : std::string arg_inp, norm_input = olab + "2: COMBINE PERIODIC=NO POWERS=2,2";
423 43 : std::string snum = getSymbol( -l );
424 86 : arg_inp = " ARG=" + ilab + sep + "r" + vlab + "-" + snum +"," + ilab + sep + "i" + vlab + "-" + snum;
425 351 : for(int i=-l+1; i<=l; ++i) {
426 308 : snum = getSymbol( i );
427 616 : arg_inp += "," + ilab + sep + "r" + vlab + "-" + snum + "," + ilab + sep + "i" + vlab + "-" + snum;
428 : norm_input += ",2,2";
429 : }
430 43 : readInputLine( norm_input + arg_inp );
431 86 : readInputLine( olab + ": CUSTOM ARG=" + olab + "2 FUNC=sqrt(x) PERIODIC=NO");
432 43 : }
433 :
434 377 : std::string Steinhardt::getSymbol( const int& m ) const {
435 377 : if( m<0 ) {
436 : std::string num;
437 166 : Tools::convert( -1*m, num );
438 166 : return "n" + num;
439 211 : } else if( m>0 ) {
440 : std::string num;
441 166 : Tools::convert( m, num );
442 166 : return "p" + num;
443 : }
444 45 : return "0";
445 : }
446 :
447 : }
448 : }
449 :
|
