Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2015-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 "core/ActionShortcut.h"
23 : #include "core/ActionRegister.h"
24 : #include "multicolvar/MultiColvarShortcuts.h"
25 : #include "core/PlumedMain.h"
26 : #include "core/ActionSet.h"
27 : #include "core/ActionWithValue.h"
28 :
29 : namespace PLMD {
30 : namespace symfunc {
31 :
32 : //+PLUMEDOC MCOLVARF LOCAL_Q1
33 : /*
34 : Calculate the local degree of order around an atoms by taking the average dot product between the q_1 vector on the central atom and the q_3 vector on the atoms in the first coordination sphere.
35 :
36 : \par Examples
37 :
38 : */
39 : //+ENDPLUMEDOC
40 :
41 : //+PLUMEDOC MCOLVARF LOCAL_Q3
42 : /*
43 : Calculate the local degree of order around an atoms by taking the average dot product between the q_3 vector on the central atom and the q_3 vector on the atoms in the first coordination sphere.
44 :
45 : The \ref Q3 command allows one to calculate one complex vectors for each of the atoms in your system that describe the degree of order in the coordination sphere
46 : around a particular atom. The difficulty with these vectors comes when combining the order parameters from all of the individual atoms/molecules so as to get a
47 : measure of the global degree of order for the system. The simplest way of doing this - calculating the average Steinhardt parameter - can be problematic. If one is
48 : examining nucleation say only the order parameters for those atoms in the nucleus will change significantly when the nucleus forms. The order parameters for the
49 : atoms in the surrounding liquid will remain pretty much the same. As such if one models a small nucleus embedded in a very large amount of solution/melt any
50 : change in the average order parameter will be negligible. Substantial changes in the value of this average can be observed in simulations of nucleation but only
51 : because the number of atoms is relatively small.
52 :
53 : When the average \ref Q3 parameter is used to bias the dynamics a problems
54 : can occur. These averaged coordinates cannot distinguish between the correct,
55 : single-nucleus pathway and a concerted pathway in which all the atoms rearrange
56 : themselves into their solid-like configuration simultaneously. This second type
57 : of pathway would be impossible in reality because there is a large entropic
58 : barrier that prevents concerted processes like this from happening. However,
59 : in the finite sized systems that are commonly simulated this barrier is reduced
60 : substantially. As a result in simulations where average Steinhardt parameters
61 : are biased there are often quite dramatic system size effects
62 :
63 : If one wants to simulate nucleation using some form on
64 : biased dynamics what is really required is an order parameter that measures:
65 :
66 : - Whether or not the coordination spheres around atoms are ordered
67 : - Whether or not the atoms that are ordered are clustered together in a crystalline nucleus
68 :
69 : \ref LOCAL_AVERAGE and \ref NLINKS are variables that can be combined with the Steinhardt parameters allow to calculate variables that satisfy these requirements.
70 : LOCAL_Q3 is another variable that can be used in these sorts of calculations. The LOCAL_Q3 parameter for a particular atom is a number that measures the extent to
71 : which the orientation of the atoms in the first coordination sphere of an atom match the orientation of the central atom. It does this by calculating the following
72 : quantity for each of the atoms in the system:
73 :
74 : \f[
75 : s_i = \frac{ \sum_j \sigma( r_{ij} ) \sum_{m=-3}^3 q_{3m}^{*}(i)q_{3m}(j) }{ \sum_j \sigma( r_{ij} ) }
76 : \f]
77 :
78 : where \f$q_{3m}(i)\f$ and \f$q_{3m}(j)\f$ are the 3rd order Steinhardt vectors calculated for atom \f$i\f$ and atom \f$j\f$ respectively and the asterisk denotes complex
79 : conjugation. The function
80 : \f$\sigma( r_{ij} )\f$ is a \ref switchingfunction that acts on the distance between atoms \f$i\f$ and \f$j\f$. The parameters of this function should be set
81 : so that it the function is equal to one when atom \f$j\f$ is in the first coordination sphere of atom \f$i\f$ and is zero otherwise. The sum in the numerator
82 : of this expression is the dot product of the Steinhardt parameters for atoms \f$i\f$ and \f$j\f$ and thus measures the degree to which the orientations of these
83 : adjacent atoms is correlated.
84 :
85 : \par Examples
86 :
87 : The following command calculates the average value of the LOCAL_Q3 parameter for the 64 Lennard Jones atoms in the system under study and prints this
88 : quantity to a file called colvar.
89 :
90 : \plumedfile
91 : Q3 SPECIES=1-64 D_0=1.3 R_0=0.2 LABEL=q3
92 : LOCAL_Q3 SPECIES=q3 SWITCH={RATIONAL D_0=1.3 R_0=0.2} MEAN LABEL=lq3
93 : PRINT ARG=lq3.mean FILE=colvar
94 : \endplumedfile
95 :
96 : The following input calculates the distribution of LOCAL_Q3 parameters at any given time and outputs this information to a file.
97 :
98 : \plumedfile
99 : Q3 SPECIES=1-64 D_0=1.3 R_0=0.2 LABEL=q3
100 : LOCAL_Q3 SPECIES=q3 SWITCH={RATIONAL D_0=1.3 R_0=0.2} HISTOGRAM={GAUSSIAN LOWER=0.0 UPPER=1.0 NBINS=20 SMEAR=0.1} LABEL=lq3
101 : PRINT ARG=lq3.* FILE=colvar
102 : \endplumedfile
103 :
104 : The following calculates the LOCAL_Q3 parameters for atoms 1-5 only. For each of these atoms comparisons of the geometry of the coordination sphere
105 : are done with those of all the other atoms in the system. The final quantity is the average and is outputted to a file
106 :
107 : \plumedfile
108 : Q3 SPECIESA=1-5 SPECIESB=1-64 D_0=1.3 R_0=0.2 LABEL=q3a
109 : Q3 SPECIESA=6-64 SPECIESB=1-64 D_0=1.3 R_0=0.2 LABEL=q3b
110 :
111 : LOCAL_Q3 SPECIES=q3a,q3b SWITCH={RATIONAL D_0=1.3 R_0=0.2} MEAN LOWMEM LABEL=w3
112 : PRINT ARG=w3.* FILE=colvar
113 : \endplumedfile
114 :
115 : */
116 : //+ENDPLUMEDOC
117 :
118 : //+PLUMEDOC MCOLVARF LOCAL_Q4
119 : /*
120 : Calculate the local degree of order around an atoms by taking the average dot product between the q_4 vector on the central atom and the q_4 vector on the atoms in the first coordination sphere.
121 :
122 : The \ref Q4 command allows one to calculate one complex vectors for each of the atoms in your system that describe the degree of order in the coordination sphere
123 : around a particular atom. The difficulty with these vectors comes when combining the order parameters from all of the individual atoms/molecules so as to get a
124 : measure of the global degree of order for the system. The simplest way of doing this - calculating the average Steinhardt parameter - can be problematic. If one is
125 : examining nucleation say only the order parameters for those atoms in the nucleus will change significantly when the nucleus forms. The order parameters for the
126 : atoms in the surrounding liquid will remain pretty much the same. As such if one models a small nucleus embedded in a very large amount of solution/melt any
127 : change in the average order parameter will be negligible. Substantial changes in the value of this average can be observed in simulations of nucleation but only
128 : because the number of atoms is relatively small.
129 :
130 : When the average \ref Q4 parameter is used to bias the dynamics a problems
131 : can occur. These averaged coordinates cannot distinguish between the correct,
132 : single-nucleus pathway and a concerted pathway in which all the atoms rearrange
133 : themselves into their solid-like configuration simultaneously. This second type
134 : of pathway would be impossible in reality because there is a large entropic
135 : barrier that prevents concerted processes like this from happening. However,
136 : in the finite sized systems that are commonly simulated this barrier is reduced
137 : substantially. As a result in simulations where average Steinhardt parameters
138 : are biased there are often quite dramatic system size effects
139 :
140 : If one wants to simulate nucleation using some form on
141 : biased dynamics what is really required is an order parameter that measures:
142 :
143 : - Whether or not the coordination spheres around atoms are ordered
144 : - Whether or not the atoms that are ordered are clustered together in a crystalline nucleus
145 :
146 : \ref LOCAL_AVERAGE and \ref NLINKS are variables that can be combined with the Steinhardt parameters allow to calculate variables that satisfy these requirements.
147 : LOCAL_Q4 is another variable that can be used in these sorts of calculations. The LOCAL_Q4 parameter for a particular atom is a number that measures the extent to
148 : which the orientation of the atoms in the first coordination sphere of an atom match the orientation of the central atom. It does this by calculating the following
149 : quantity for each of the atoms in the system:
150 :
151 : \f[
152 : s_i = \frac{ \sum_j \sigma( r_{ij} ) \sum_{m=-4}^4 q_{4m}^{*}(i)q_{4m}(j) }{ \sum_j \sigma( r_{ij} ) }
153 : \f]
154 :
155 : where \f$q_{4m}(i)\f$ and \f$q_{4m}(j)\f$ are the fourth order Steinhardt vectors calculated for atom \f$i\f$ and atom \f$j\f$ respectively and the asterisk denotes
156 : complex conjugation. The function
157 : \f$\sigma( r_{ij} )\f$ is a \ref switchingfunction that acts on the distance between atoms \f$i\f$ and \f$j\f$. The parameters of this function should be set
158 : so that it the function is equal to one when atom \f$j\f$ is in the first coordination sphere of atom \f$i\f$ and is zero otherwise. The sum in the numerator
159 : of this expression is the dot product of the Steinhardt parameters for atoms \f$i\f$ and \f$j\f$ and thus measures the degree to which the orientations of these
160 : adjacent atoms is correlated.
161 :
162 : \par Examples
163 :
164 : The following command calculates the average value of the LOCAL_Q4 parameter for the 64 Lennard Jones atoms in the system under study and prints this
165 : quantity to a file called colvar.
166 :
167 : \plumedfile
168 : Q4 SPECIES=1-64 D_0=1.3 R_0=0.2 LABEL=q4
169 : LOCAL_Q4 SPECIES=q4 SWITCH={RATIONAL D_0=1.3 R_0=0.2} MEAN LABEL=lq4
170 : PRINT ARG=lq4.mean FILE=colvar
171 : \endplumedfile
172 :
173 : The following input calculates the distribution of LOCAL_Q4 parameters at any given time and outputs this information to a file.
174 :
175 : \plumedfile
176 : Q4 SPECIES=1-64 D_0=1.3 R_0=0.2 LABEL=q4
177 : LOCAL_Q4 SPECIES=q4 SWITCH={RATIONAL D_0=1.3 R_0=0.2} HISTOGRAM={GAUSSIAN LOWER=0.0 UPPER=1.0 NBINS=20 SMEAR=0.1} LABEL=lq4
178 : PRINT ARG=lq4.* FILE=colvar
179 : \endplumedfile
180 :
181 : The following calculates the LOCAL_Q4 parameters for atoms 1-5 only. For each of these atoms comparisons of the geometry of the coordination sphere
182 : are done with those of all the other atoms in the system. The final quantity is the average and is outputted to a file
183 :
184 : \plumedfile
185 : Q4 SPECIESA=1-5 SPECIESB=1-64 D_0=1.3 R_0=0.2 LABEL=q4a
186 : Q4 SPECIESA=6-64 SPECIESB=1-64 D_0=1.3 R_0=0.2 LABEL=q4b
187 :
188 : LOCAL_Q4 SPECIES=q4a,q4b SWITCH={RATIONAL D_0=1.3 R_0=0.2} MEAN LOWMEM LABEL=w4
189 : PRINT ARG=w4.* FILE=colvar
190 : \endplumedfile
191 :
192 : */
193 : //+ENDPLUMEDOC
194 :
195 : //+PLUMEDOC MCOLVARF LOCAL_Q6
196 : /*
197 : Calculate the local degree of order around an atoms by taking the average dot product between the q_6 vector on the central atom and the q_6 vector on the atoms in the first coordination sphere.
198 :
199 : The \ref Q6 command allows one to calculate one complex vectors for each of the atoms in your system that describe the degree of order in the coordination sphere
200 : around a particular atom. The difficulty with these vectors comes when combining the order parameters from all of the individual atoms/molecules so as to get a
201 : measure of the global degree of order for the system. The simplest way of doing this - calculating the average Steinhardt parameter - can be problematic. If one is
202 : examining nucleation say only the order parameters for those atoms in the nucleus will change significantly when the nucleus forms. The order parameters for the
203 : atoms in the surrounding liquid will remain pretty much the same. As such if one models a small nucleus embedded in a very large amount of solution/melt any
204 : change in the average order parameter will be negligible. Substantial changes in the value of this average can be observed in simulations of nucleation but only
205 : because the number of atoms is relatively small.
206 :
207 : When the average \ref Q6 parameter is used to bias the dynamics a problems
208 : can occur. These averaged coordinates cannot distinguish between the correct,
209 : single-nucleus pathway and a concerted pathway in which all the atoms rearrange
210 : themselves into their solid-like configuration simultaneously. This second type
211 : of pathway would be impossible in reality because there is a large entropic
212 : barrier that prevents concerted processes like this from happening. However,
213 : in the finite sized systems that are commonly simulated this barrier is reduced
214 : substantially. As a result in simulations where average Steinhardt parameters
215 : are biased there are often quite dramatic system size effects
216 :
217 : If one wants to simulate nucleation using some form on
218 : biased dynamics what is really required is an order parameter that measures:
219 :
220 : - Whether or not the coordination spheres around atoms are ordered
221 : - Whether or not the atoms that are ordered are clustered together in a crystalline nucleus
222 :
223 : \ref LOCAL_AVERAGE and \ref NLINKS are variables that can be combined with the Steinhardt parameters allow to calculate variables that satisfy these requirements.
224 : LOCAL_Q6 is another variable that can be used in these sorts of calculations. The LOCAL_Q6 parameter for a particular atom is a number that measures the extent to
225 : which the orientation of the atoms in the first coordination sphere of an atom match the orientation of the central atom. It does this by calculating the following
226 : quantity for each of the atoms in the system:
227 :
228 : \f[
229 : s_i = \frac{ \sum_j \sigma( r_{ij} ) \sum_{m=-6}^6 q_{6m}^{*}(i)q_{6m}(j) }{ \sum_j \sigma( r_{ij} ) }
230 : \f]
231 :
232 : where \f$q_{6m}(i)\f$ and \f$q_{6m}(j)\f$ are the sixth order Steinhardt vectors calculated for atom \f$i\f$ and atom \f$j\f$ respectively and the asterisk denotes
233 : complex conjugation. The function
234 : \f$\sigma( r_{ij} )\f$ is a \ref switchingfunction that acts on the distance between atoms \f$i\f$ and \f$j\f$. The parameters of this function should be set
235 : so that it the function is equal to one when atom \f$j\f$ is in the first coordination sphere of atom \f$i\f$ and is zero otherwise. The sum in the numerator
236 : of this expression is the dot product of the Steinhardt parameters for atoms \f$i\f$ and \f$j\f$ and thus measures the degree to which the orientations of these
237 : adjacent atoms is correlated.
238 :
239 : \par Examples
240 :
241 : The following command calculates the average value of the LOCAL_Q6 parameter for the 64 Lennard Jones atoms in the system under study and prints this
242 : quantity to a file called colvar.
243 :
244 : \plumedfile
245 : Q6 SPECIES=1-64 D_0=1.3 R_0=0.2 LABEL=q6
246 : LOCAL_Q6 SPECIES=q6 SWITCH={RATIONAL D_0=1.3 R_0=0.2} MEAN LABEL=lq6
247 : PRINT ARG=lq6.mean FILE=colvar
248 : \endplumedfile
249 :
250 : The following input calculates the distribution of LOCAL_Q6 parameters at any given time and outputs this information to a file.
251 :
252 : \plumedfile
253 : Q6 SPECIES=1-64 D_0=1.3 R_0=0.2 LABEL=q6
254 : LOCAL_Q6 SPECIES=q6 SWITCH={RATIONAL D_0=1.3 R_0=0.2} HISTOGRAM={GAUSSIAN LOWER=0.0 UPPER=1.0 NBINS=20 SMEAR=0.1} LABEL=lq6
255 : PRINT ARG=lq6.* FILE=colvar
256 : \endplumedfile
257 :
258 : The following calculates the LOCAL_Q6 parameters for atoms 1-5 only. For each of these atoms comparisons of the geometry of the coordination sphere
259 : are done with those of all the other atoms in the system. The final quantity is the average and is outputted to a file
260 :
261 : \plumedfile
262 : Q6 SPECIESA=1-5 SPECIESB=1-64 D_0=1.3 R_0=0.2 LABEL=q6a
263 : Q6 SPECIESA=6-64 SPECIESB=1-64 D_0=1.3 R_0=0.2 LABEL=q6b
264 :
265 : LOCAL_Q6 SPECIES=q6a,q6b SWITCH={RATIONAL D_0=1.3 R_0=0.2} MEAN LOWMEM LABEL=w6
266 : PRINT ARG=w6.* FILE=colvar
267 : \endplumedfile
268 :
269 : */
270 : //+ENDPLUMEDOC
271 :
272 : class LocalSteinhardt : public ActionShortcut {
273 : private:
274 : std::string getSymbol( const int& m ) const ;
275 : std::string getArgsForStack( const int& l, const std::string& lab ) const;
276 : public:
277 : static void registerKeywords( Keywords& keys );
278 : explicit LocalSteinhardt(const ActionOptions&);
279 : };
280 :
281 : PLUMED_REGISTER_ACTION(LocalSteinhardt,"LOCAL_Q1")
282 : PLUMED_REGISTER_ACTION(LocalSteinhardt,"LOCAL_Q3")
283 : PLUMED_REGISTER_ACTION(LocalSteinhardt,"LOCAL_Q4")
284 : PLUMED_REGISTER_ACTION(LocalSteinhardt,"LOCAL_Q6")
285 :
286 20 : void LocalSteinhardt::registerKeywords( Keywords& keys ) {
287 20 : ActionShortcut::registerKeywords( keys );
288 40 : keys.add("optional","SPECIES","");
289 40 : keys.add("optional","SPECIESA","");
290 40 : keys.add("optional","SPECIESB","");
291 40 : keys.add("optional","SWITCH","This keyword is used if you want to employ an alternative to the continuous swiching function defined above. "
292 : "The following provides information on the \\ref switchingfunction that are available. "
293 : "When this keyword is present you no longer need the NN, MM, D_0 and R_0 keywords.");
294 40 : keys.addFlag("LOWMEM",false,"this flag does nothing and is present only to ensure back-compatibility");
295 20 : multicolvar::MultiColvarShortcuts::shortcutKeywords( keys );
296 20 : keys.needsAction("CONTACT_MATRIX");
297 20 : keys.needsAction("MATRIX_PRODUCT");
298 20 : keys.needsAction("GROUP");
299 20 : keys.needsAction("ONES");
300 20 : keys.needsAction("OUTER_PRODUCT");
301 20 : keys.needsAction("VSTACK");
302 20 : keys.needsAction("CONCATENATE");
303 20 : keys.needsAction("CUSTOM");
304 20 : keys.needsAction("TRANSPOSE");
305 20 : keys.needsAction("MATRIX_VECTOR_PRODUCT");
306 20 : }
307 :
308 98 : std::string LocalSteinhardt::getSymbol( const int& m ) const {
309 98 : if( m<0 ) {
310 : std::string num;
311 40 : Tools::convert( -1*m, num );
312 40 : return "n" + num;
313 58 : } else if( m>0 ) {
314 : std::string num;
315 49 : Tools::convert( m, num );
316 49 : return "p" + num;
317 : }
318 9 : return "0";
319 : }
320 :
321 9 : std::string LocalSteinhardt::getArgsForStack( const int& l, const std::string& sp_lab ) const {
322 : std::string numstr;
323 9 : Tools::convert( l, numstr );
324 18 : std::string data_mat = " ARG=" + sp_lab + "_sp.rm-n" + numstr + "," + sp_lab + "_sp.im-n" + numstr;
325 107 : for(int i=-l+1; i<=l; ++i) {
326 98 : numstr = getSymbol( i );
327 196 : data_mat += "," + sp_lab + "_sp.rm-" + numstr + "," + sp_lab + "_sp.im-" + numstr;
328 : }
329 9 : return data_mat;
330 : }
331 :
332 7 : LocalSteinhardt::LocalSteinhardt(const ActionOptions& ao):
333 : Action(ao),
334 7 : ActionShortcut(ao) {
335 : bool lowmem;
336 7 : parseFlag("LOWMEM",lowmem);
337 7 : if( lowmem ) {
338 0 : warning("LOWMEM flag is deprecated and is no longer required for this action");
339 : }
340 : // Get the Q value
341 : int l;
342 14 : Tools::convert( getName().substr(7), l);
343 : // Create a vector filled with ones
344 : std::string twolplusone;
345 7 : Tools::convert( 2*(2*l+1), twolplusone );
346 14 : readInputLine( getShortcutLabel() + "_uvec: ONES SIZE=" + twolplusone );
347 : // Read in species keyword
348 : std::string sp_str;
349 14 : parse("SPECIES",sp_str);
350 : std::string spa_str;
351 14 : parse("SPECIESA",spa_str);
352 7 : if( sp_str.length()>0 ) {
353 : // Create a group with these atoms
354 12 : readInputLine( getShortcutLabel() + "_grp: GROUP ATOMS=" + sp_str );
355 6 : std::vector<std::string> sp_lab = Tools::getWords(sp_str, "\t\n ,");
356 : // This creates the stash to hold all the vectors
357 6 : if( sp_lab.size()==1 ) {
358 : // The lengths of all the vectors in a vector
359 12 : readInputLine( getShortcutLabel() + "_nmat: OUTER_PRODUCT ARG=" + sp_lab[0] + "_norm," + getShortcutLabel() + "_uvec");
360 : // The unormalised vectors
361 12 : readInputLine( getShortcutLabel() + "_uvecs: VSTACK" + getArgsForStack( l, sp_lab[0] ) );
362 : } else {
363 0 : std::string len_vec = getShortcutLabel() + "_mags: CONCATENATE ARG=" + sp_lab[0] + "_norm";
364 0 : for(unsigned i=1; i<sp_lab.size(); ++i) {
365 0 : len_vec += "," + sp_lab[i] + "_norm";
366 : }
367 : // This is the vector that contains all the magnitudes
368 0 : readInputLine( len_vec );
369 0 : std::string concat_str = getShortcutLabel() + "_uvecs: CONCATENATE";
370 0 : for(unsigned i=0; i<sp_lab.size(); ++i) {
371 : std::string snum;
372 0 : Tools::convert( i+1, snum );
373 0 : concat_str += " MATRIX" + snum + "1=" + getShortcutLabel() + "_uvecs" + snum;
374 0 : readInputLine( getShortcutLabel() + "_uvecs" + snum + ": VSTACK" + getArgsForStack( l, sp_lab[i] ) );
375 : }
376 : // And the normalising matrix by taking the column vector of magnitudes and multiplying by the row vector of ones
377 0 : readInputLine( getShortcutLabel() + "_nmat: OUTER_PRODUCT ARG=" + getShortcutLabel() + "_mags," + getShortcutLabel() + "_uvec");
378 : // The unormalised vectors
379 0 : readInputLine( concat_str );
380 : }
381 : // Now normalise all the vectors by doing Hadammard "product" with normalising matrix
382 12 : readInputLine( getShortcutLabel() + "_vecs: CUSTOM ARG=" + getShortcutLabel() + "_uvecs," + getShortcutLabel() + "_nmat FUNC=x/y PERIODIC=NO");
383 : // And transpose the matrix
384 12 : readInputLine( getShortcutLabel() + "_vecsT: TRANSPOSE ARG=" + getShortcutLabel() + "_vecs" );
385 : std::string sw_str;
386 6 : parse("SWITCH",sw_str);
387 12 : readInputLine( getShortcutLabel() + "_cmap: CONTACT_MATRIX GROUP=" + sp_str + " SWITCH={" + sw_str + "}");
388 : // And the matrix of dot products
389 12 : readInputLine( getShortcutLabel() + "_dpmat: MATRIX_PRODUCT ARG=" + getShortcutLabel() + "_vecs," + getShortcutLabel() + "_vecsT" );
390 7 : } else if( spa_str.length()>0 ) {
391 : // Create a group with these atoms
392 2 : readInputLine( getShortcutLabel() + "_grp: GROUP ATOMS=" + spa_str );
393 : std::string spb_str;
394 2 : parse("SPECIESB",spb_str);
395 1 : if( spb_str.length()==0 ) {
396 0 : plumed_merror("need both SPECIESA and SPECIESB in input");
397 : }
398 1 : std::vector<std::string> sp_laba = Tools::getWords(spa_str, "\t\n ,");
399 1 : std::vector<std::string> sp_labb = Tools::getWords(spb_str, "\t\n ,");
400 1 : if( sp_laba.size()==1 ) {
401 : // The matrix that is used for normalising
402 2 : readInputLine( getShortcutLabel() + "_nmatA: OUTER_PRODUCT ARG=" + sp_laba[0] + "_norm," + getShortcutLabel() + "_uvec");
403 : // The unormalised vectors
404 2 : readInputLine( getShortcutLabel() + "_uvecsA: VSTACK" + getArgsForStack( l, sp_laba[0] ) );
405 : } else {
406 0 : std::string len_vec = getShortcutLabel() + "_magsA: CONCATENATE ARG=" + sp_laba[0] + "_norm";
407 0 : for(unsigned i=1; i<sp_laba.size(); ++i) {
408 0 : len_vec += "," + sp_laba[i] + "_norm";
409 : }
410 : // This is the vector that contains all the magnitudes
411 0 : readInputLine( len_vec );
412 0 : std::string concat_str = getShortcutLabel() + "_uvecsA: CONCATENATE";
413 0 : for(unsigned i=0; i<sp_laba.size(); ++i) {
414 : std::string snum;
415 0 : Tools::convert( i+1, snum );
416 0 : concat_str += " MATRIX" + snum + "1=" + getShortcutLabel() + "_uvecsA" + snum;
417 0 : readInputLine( getShortcutLabel() + "_uvecsA" + snum + ": VSTACK" + getArgsForStack( l, sp_laba[i] ) );
418 : }
419 : // And the normalising matrix by taking the column vector of magnitudes and multiplying by the row vector of ones
420 0 : readInputLine( getShortcutLabel() + "_nmatA: OUTER_PRODUCT ARG=" + getShortcutLabel() + "_magsA," + getShortcutLabel() + "_uvec");
421 : // The unormalised vector
422 0 : readInputLine( concat_str );
423 : }
424 : // Now normalise all the vectors by doing Hadammard "product" with normalising matrix
425 2 : readInputLine( getShortcutLabel() + "_vecsA: CUSTOM ARG=" + getShortcutLabel() + "_uvecsA," + getShortcutLabel() + "_nmatA FUNC=x/y PERIODIC=NO");
426 : // Now do second matrix
427 1 : if( sp_labb.size()==1 ) {
428 0 : readInputLine( getShortcutLabel() + "_nmatB: OUTER_PRODUCT ARG=" + sp_laba[0] + "_norm," + getShortcutLabel() + "_uvec");
429 0 : readInputLine( getShortcutLabel() + "_uvecsBT: VSTACK" + getArgsForStack( l, sp_labb[0] ) );
430 0 : readInputLine( getShortcutLabel() + "_uvecsB: TRANSPOSE ARG=" + getShortcutLabel() + "_uvecsBT");
431 : } else {
432 2 : std::string len_vec = getShortcutLabel() + "_magsB: CONCATENATE ARG=" + sp_labb[0] + "_norm";
433 2 : for(unsigned i=1; i<sp_labb.size(); ++i) {
434 2 : len_vec += "," + sp_labb[i] + "_norm";
435 : }
436 : // This is the vector that contains all the magnitudes
437 1 : readInputLine( len_vec );
438 1 : std::string concat_str = getShortcutLabel() + "_uvecsB: CONCATENATE";
439 3 : for(unsigned i=0; i<sp_labb.size(); ++i) {
440 : std::string snum;
441 2 : Tools::convert( i+1, snum );
442 4 : concat_str += " MATRIX1" + snum + "=" + getShortcutLabel() + "_uvecsB" + snum;
443 4 : readInputLine( getShortcutLabel() + "_uvecsBT" + snum + ": VSTACK" + getArgsForStack( l, sp_labb[i] ) );
444 4 : readInputLine( getShortcutLabel() + "_uvecsB" + snum + ": TRANSPOSE ARG=" + getShortcutLabel() + "_uvecsBT" + snum );
445 : }
446 : // And the normalising matrix
447 2 : readInputLine( getShortcutLabel() + "_nmatB: OUTER_PRODUCT ARG=" + getShortcutLabel() + "_uvec," + getShortcutLabel() + "_magsB");
448 : // The unormalised vectors
449 1 : readInputLine( concat_str );
450 : }
451 : // Now normalise all the vectors by doing Hadammard "product" with normalising matrix
452 2 : readInputLine( getShortcutLabel() + "_vecsB: CUSTOM ARG=" + getShortcutLabel() + "_uvecsB," + getShortcutLabel() + "_nmatB FUNC=x/y PERIODIC=NO");
453 : std::string sw_str;
454 1 : parse("SWITCH",sw_str);
455 2 : readInputLine( getShortcutLabel() + "_cmap: CONTACT_MATRIX GROUPA=" + spa_str + " GROUPB=" + spb_str + " SWITCH={" + sw_str + "}");
456 2 : readInputLine( getShortcutLabel() + "_dpmat: MATRIX_PRODUCT ARG=" + getShortcutLabel() + "_vecsA," + getShortcutLabel() + "_vecsB" );
457 1 : }
458 :
459 : // Now create the product matrix
460 14 : readInputLine( getShortcutLabel() + "_prod: CUSTOM ARG=" + getShortcutLabel() + "_cmap," + getShortcutLabel() + "_dpmat FUNC=x*y PERIODIC=NO");
461 : // Now the sum of coordination numbers times the switching functions
462 7 : ActionWithValue* av = plumed.getActionSet().selectWithLabel<ActionWithValue*>( getShortcutLabel() + "_cmap");
463 7 : plumed_assert( av && av->getNumberOfComponents()>0 && (av->copyOutput(0))->getRank()==2 );
464 : std::string size;
465 7 : Tools::convert( (av->copyOutput(0))->getShape()[1], size );
466 14 : readInputLine( getShortcutLabel() + "_ones: ONES SIZE=" + size );
467 14 : readInputLine( getShortcutLabel() + ": MATRIX_VECTOR_PRODUCT ARG=" + getShortcutLabel() +"_prod," + getShortcutLabel() +"_ones");
468 : // And just the sum of the coordination numbers
469 14 : readInputLine( getShortcutLabel() + "_denom: MATRIX_VECTOR_PRODUCT ARG=" + getShortcutLabel() + "_cmap," + getShortcutLabel() +"_ones");
470 : // And matheval to get the final quantity
471 14 : readInputLine( getShortcutLabel() + "_av: CUSTOM ARG=" + getShortcutLabel() + "," + getShortcutLabel() + "_denom FUNC=x/y PERIODIC=NO");
472 : // And this expands everything
473 14 : multicolvar::MultiColvarShortcuts::expandFunctions( getShortcutLabel(), getShortcutLabel() + "_av", "", this );
474 7 : }
475 :
476 : }
477 : }
|
