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 "core/ActionShortcut.h"
23 : #include "multicolvar/MultiColvarShortcuts.h"
24 : #include "core/PlumedMain.h"
25 : #include "core/ActionSet.h"
26 : #include "core/ActionRegister.h"
27 : #include "core/ActionWithValue.h"
28 : #include "CoordinationNumbers.h"
29 :
30 : #include <complex>
31 :
32 : namespace PLMD {
33 : namespace symfunc {
34 :
35 : //+PLUMEDOC MCOLVAR HEXACTIC_PARAMETER
36 : /*
37 : Calculate the hexatic order parameter
38 :
39 : This [symmetry function](https://www.plumed-tutorials.org/lessons/23/001/data/SymmetryFunction.html) can be used to understand phase transitions in two dimensional systems.
40 : The symmetry function for atom $k$ is calculated using:
41 :
42 : $$
43 : s_k = \left| \frac{\sum_j \sigma(r_{kj}) e^{6i\theta_j} }{\sum_j \sigma(r_{kj})} \right|
44 : $$
45 :
46 : In this expression, the sum run over all the atoms of interest and $r_{kj}$ is the distance between atom $k$ and atom $j$ and $\sigma$ is
47 : a switching function that acts upon this distance. $\theta_j$ is the angle between either the $x$, $y$ or $z$ axis and the bond connecting
48 : atom $k$ and atom $j$. This angle is multiplied by the imaginary number $i$ - the square root of minus one. In the code, we thus calculate
49 : $e^{i\theta_j}$ as follows:
50 :
51 : $$
52 : e^{i\theta_j} = \frac{x_{kj}}{r_{kj}} + i \frac{y_{kj}}{r_{kj}}
53 : $$
54 :
55 : We then take the 6th power of this complex number directly before compupting the magnitude by multiplying the result by its complex conjugate.
56 : Notice, furthermore, that we can replace $x_{kj}$ or $y_{kj}$ with $z_{kj}$ by using PLANE=xz or PLANE=yz in place of PLANE=xy.
57 :
58 : An example that shows how you can use this shortcut is shown below:
59 :
60 : ```plumed
61 : hex: HEXACTIC_PARAMETER SPECIES=1-400 PLANE=xy R_0=0.2 D_0=1.4 NN=6 MM=12
62 : hex_mean: MEAN ARG=hex_norm PERIODIC=NO
63 : PRINT ARG=hex_mean FILE=colvar
64 : ```
65 :
66 : As you can see if you expand the shortcut above, this input calculates the quantity defined in the equation above for the 400 atoms in the simulated system and stores them in a vector.
67 : The elements of this vector are then added together so the mean value can be computed.
68 :
69 : In the input above we use a rational [switching function](LESS_THAN.md) with the parameters above. We would recommend using SWITCH syntax
70 : rather than the syntax above when giving the parameters for the switching function as you can then use any of the switching functions described
71 : in the documentation for [LESS_THAN](LESS_THAN.md). More importantly, however, using this syntax allows you to set the D_MAX parameter for the
72 : switching function as demonstrated below:
73 :
74 :
75 : ```plumed
76 : hex: HEXACTIC_PARAMETER SPECIES=1-400 PLANE=xy SWITCH={RATIONAL D_0=1.4 R_0=0.2 D_MAX=3.0}
77 : hex_mean: MEAN ARG=hex_norm PERIODIC=NO
78 : PRINT ARG=hex_mean FILE=colvar
79 : ```
80 :
81 : Setting the `D_MAX` can substantially improve PLUMED performance as it turns on the linked list algorithm that is discussed in the optimisation details part
82 : of the documentation for [CONTACT_MATRIX](CONTACT_MATRIX.md).
83 :
84 : ## The VMEAN and VSUM options
85 :
86 : If you expand the inputs in the previous section you will see that the value `hex_norm` that we are calculating the average of is a vector that contains the magnitude of the
87 : complex number defined in the equation above. If you would like to add up the vector of complex numbers directly (or take an average of the vector) __before__ computing the
88 : magnitude of the complex number you can use the VSUM and VMEAN flags as sillstrated below:
89 :
90 : ```plumed
91 : hex: HEXACTIC_PARAMETER SPECIES=1-400 PLANE=xy SWITCH={RATIONAL D_0=1.4 R_0=0.2} VMEAN VSUM
92 : PRINT ARG=hex_vmean,hex_vsum FILE=colvar
93 : ```
94 :
95 : If you expand the shortcut in the input above you will see that this input adds the complex numbers together directly before calculating the modulus. This contrasts with the approach
96 : in the inputs above, which computes a vector that contains the square moduli for each of the individual hexatic order parameters and then computes the average from this vector.
97 :
98 : ## Using two types of atom
99 :
100 : If you would like to calculate the hexatic parameters using the directors of the bonds connecting the atoms in GROUPA to the atoms in GROUPB you can use an input like the one
101 : shown below:
102 :
103 : ```plumed
104 : hex: HEXACTIC_PARAMETER SPECIESA=1-200 SPECIESB=201-400 PLANE=xy SWITCH={RATIONAL D_0=1.4 R_0=0.2 D_MAX=3.0}
105 : hex_mean: MEAN ARG=hex_norm PERIODIC=NO
106 : PRINT ARG=hex_mean FILE=colvar
107 : ```
108 :
109 : ## The MASK keyword
110 :
111 : You can use the MASK keyword with this action in the same way that it is used with [COORDINATIONNUMBER](COORDINATIONNUMBER.md). This keyword thus expects a vector in
112 : input, which tells PLUMED the atoms for which you do not need to calculate the function. As illustrated below, this is useful if you are using functionality
113 : from the [volumes module](module_volumes.md) to calculate the average value of the hexatic parameter for only those atoms that lie in a certain part of the simulation box.
114 :
115 : ```plumed
116 : # Fixed virtual atom which serves as the probe volume's center (pos. in nm)
117 : center: FIXEDATOM AT=2.5,2.5,2.5
118 : # Vector in which element i is one if atom i is in sphere of interest and zero otherwise
119 : sphere: INSPHERE ATOMS=1-400 CENTER=center RADIUS={GAUSSIAN D_0=0.5 R_0=0.01 D_MAX=0.52}
120 : # Calculate the tetrahedral parameter of the atoms
121 : hex: HEXACTIC_PARAMETER ...
122 : SPECIES=1-400 MASK=sphere PLANE=xy
123 : SWITCH={RATIONAL D_0=3.0 R_0=1.5 D_MAX=6.0}
124 : ...
125 : # Multiply hexatic parameters by sphere vector
126 : prod_rm: CUSTOM ARG=hex_rm,sphere FUNC=x*y PERIODIC=NO
127 : prod_im: CUSTOM ARG=hex_im,sphere FUNC=x*y PERIODIC=NO
128 : # Sum of coordination numbers for atoms that are in the sphere of interest
129 : numer_rm: SUM ARG=prod_rm PERIODIC=NO
130 : numer_im: SUM ARG=prod_im PERIODIC=NO
131 : # Number of atoms that are in sphere of interest
132 : denom: SUM ARG=sphere PERIODIC=NO
133 : # Average coordination number for atoms in sphere of interest
134 : av_rm: CUSTOM ARG=numer_rm,denom FUNC=x/y PERIODIC=NO
135 : av_im: CUSTOM ARG=numer_im,denom FUNC=x/y PERIODIC=NO
136 : # Take the square modulus
137 : av: CUSTOM ARG=av_rm,av_im FUNC=x*x+y*y PERIODIC=NO
138 : # And print out final CV to a file
139 : PRINT ARG=av FILE=colvar STRIDE=1
140 : ```
141 :
142 : This input calculate the average value of the (complex) hexatic parameter for only those atoms that are within a spherical region that is centered on the point
143 : $(2.5,2.5,2.5)$. The square modulus of this average is then output.
144 :
145 : ## Using nearest neighbours
146 :
147 : In papers where symmetry functions similar to this one have been used a switching function is not employed. The sums over $j$ in the expression above are replaced by sums over the
148 : six nearest neighbours to each atom. If you would like to calculate this quantity using PLUMED you can use an input like this:
149 :
150 : ```plumed
151 : dmat: DISTANCE_MATRIX GROUP=1-400 CUTOFF=3.0 COMPONENTS
152 : neigh: NEIGHBORS ARG=dmat.w NLOWEST=6
153 : harm: CYLINDRICAL_HARMONIC DEGREE=6 ARG=dmat.x,dmat.y
154 : rprod: CUSTOM ARG=neigh,harm.rm FUNC=x*y PERIODIC=NO
155 : iprod: CUSTOM ARG=neigh,harm.im FUNC=x*y PERIODIC=NO
156 : hex2_ones: ONES SIZE=400
157 : hex2_denom: MATRIX_VECTOR_PRODUCT ARG=neigh,hex2_ones
158 : harm_rm: MATRIX_VECTOR_PRODUCT ARG=rprod,hex2_ones
159 : harm_im: MATRIX_VECTOR_PRODUCT ARG=iprod,hex2_ones
160 : hex2_rmn: CUSTOM ARG=harm_rm,hex2_denom FUNC=x/y PERIODIC=NO
161 : hex2_imn: CUSTOM ARG=harm_im,hex2_denom FUNC=x/y PERIODIC=NO
162 : DUMPATOMS ATOMS=1-400 ARG=hex2_rmn,hex2_imn,hex2_denom FILE=hexparam.xyz
163 : ```
164 :
165 : This input outputs the values of the order parameters for all the atoms to an extended xyz file .
166 :
167 : !!! warning "Broken virial"
168 :
169 : Virial is not working currently
170 :
171 :
172 : */
173 : //+ENDPLUMEDOC
174 :
175 :
176 : class HexacticParameter : public ActionShortcut {
177 : private:
178 : void createVectorNormInput( const std::string& ilab, const std::string& olab, const std::string& vlab );
179 : public:
180 : static void registerKeywords( Keywords& keys );
181 : explicit HexacticParameter(const ActionOptions&);
182 : };
183 :
184 : PLUMED_REGISTER_ACTION(HexacticParameter,"HEXACTIC_PARAMETER")
185 :
186 5 : void HexacticParameter::registerKeywords( Keywords& keys ) {
187 5 : CoordinationNumbers::shortcutKeywords( keys );
188 5 : keys.add("compulsory","PLANE","the plane to use when calculating the value of the order parameter should be xy, xz or yz");
189 10 : keys.setValueDescription("matrix","the value of the cylindrical harmonic for each bond vector specified");
190 5 : keys.addFlag("VMEAN",false,"calculate the norm of the mean vector.");
191 10 : keys.addOutputComponent("_vmean","VMEAN","scalar","the norm of the mean vector");
192 5 : keys.addFlag("VSUM",false,"calculate the norm of the sum of all the vectors");
193 10 : keys.addOutputComponent("_vsum","VSUM","scalar","the norm of the mean vector");
194 5 : keys.needsAction("CYLINDRICAL_HARMONIC_MATRIX");
195 5 : keys.needsAction("ONES");
196 5 : keys.needsAction("MATRIX_VECTOR_PRODUCT");
197 5 : keys.needsAction("CUSTOM");
198 5 : keys.needsAction("MEAN");
199 5 : keys.needsAction("SUM");
200 5 : keys.needsAction("COMBINE");
201 5 : keys.addDOI("10.1103/PhysRevLett.99.215701");
202 5 : }
203 :
204 1 : HexacticParameter::HexacticParameter( const ActionOptions& ao):
205 : Action(ao),
206 1 : ActionShortcut(ao) {
207 : std::string sp_str, specA, specB;
208 1 : parse("SPECIES",sp_str);
209 1 : parse("SPECIESA",specA);
210 1 : parse("SPECIESB",specB);
211 1 : CoordinationNumbers::expandMatrix( true, getShortcutLabel(), sp_str, specA, specB, this );
212 : std::string myplane;
213 2 : parse("PLANE",myplane);
214 1 : if( myplane=="xy" ) {
215 2 : readInputLine( getShortcutLabel() + ": CYLINDRICAL_HARMONIC_MATRIX DEGREE=6 ARG=" + getShortcutLabel() + "_mat.x," + getShortcutLabel() + "_mat.y," + getShortcutLabel() + "_mat.w" );
216 0 : } else if( myplane=="xz" ) {
217 0 : readInputLine( getShortcutLabel() + ": CYLINDRICAL_HARMONIC_MATRIX DEGREE=6 ARG=" + getShortcutLabel() + "_mat.x," + getShortcutLabel() + "_mat.z," + getShortcutLabel() + "_mat.w" );
218 0 : } else if( myplane=="yz" ) {
219 0 : readInputLine( getShortcutLabel() + ": CYLINDRICAL_HARMONIC_MATRIX DEGREE=6 ARG=" + getShortcutLabel() + "_mat.y," + getShortcutLabel() + "_mat.z," + getShortcutLabel() + "_mat.w" );
220 : } else {
221 0 : error("invalid input for plane -- should be xy, xz or yz");
222 : }
223 : // And coordination number
224 1 : ActionWithValue* av = plumed.getActionSet().selectWithLabel<ActionWithValue*>( getShortcutLabel() + "_mat");
225 1 : plumed_assert( av && av->getNumberOfComponents()>0 && (av->copyOutput(0))->getRank()==2 );
226 : std::string size;
227 1 : Tools::convert( (av->copyOutput(0))->getShape()[1], size );
228 2 : readInputLine( getShortcutLabel() + "_ones: ONES SIZE=" + size );
229 2 : readInputLine( getShortcutLabel() + "_rm: MATRIX_VECTOR_PRODUCT ARG=" + getShortcutLabel() + ".rm," + getShortcutLabel() + "_ones");
230 2 : readInputLine( getShortcutLabel() + "_im: MATRIX_VECTOR_PRODUCT ARG=" + getShortcutLabel() + ".im," + getShortcutLabel() + "_ones");
231 : // Input for denominator (coord)
232 2 : readInputLine( getShortcutLabel() + "_denom: MATRIX_VECTOR_PRODUCT ARG=" + getShortcutLabel() + "_mat.w," + getShortcutLabel() + "_ones");
233 : // Divide real part by coordination numbers
234 2 : readInputLine( getShortcutLabel() + "_rmn: CUSTOM ARG=" + getShortcutLabel() + "_rm," + getShortcutLabel() + "_denom FUNC=x/y PERIODIC=NO");
235 : // Devide imaginary part by coordination number
236 2 : readInputLine( getShortcutLabel() + "_imn: CUSTOM ARG=" + getShortcutLabel() + "_im," + getShortcutLabel() + "_denom FUNC=x/y PERIODIC=NO");
237 :
238 : // If we are doing VMEAN determine sum of vector components
239 : bool do_vmean;
240 1 : parseFlag("VMEAN",do_vmean);
241 1 : if( do_vmean ) {
242 : // Real part
243 0 : readInputLine( getShortcutLabel() + "_rms: MEAN ARG=" + getShortcutLabel() + "_rmn PERIODIC=NO");
244 : // Imaginary part
245 0 : readInputLine( getShortcutLabel() + "_ims: MEAN ARG=" + getShortcutLabel() + "_imn PERIODIC=NO");
246 : // Now calculate the total length of the vector
247 0 : createVectorNormInput( getShortcutLabel(), getShortcutLabel() + "_vmean", "ms" );
248 : }
249 : bool do_vsum;
250 1 : parseFlag("VSUM",do_vsum);
251 1 : if( do_vsum ) {
252 : // Real part
253 0 : readInputLine( getShortcutLabel() + "_rmz: SUM ARG=" + getShortcutLabel() + "_rmn PERIODIC=NO");
254 : // Imaginary part
255 0 : readInputLine( getShortcutLabel() + "_imz: SUM ARG=" + getShortcutLabel() + "_imn PERIODIC=NO");
256 : // Now calculate the total length of the vector
257 0 : createVectorNormInput( getShortcutLabel(), getShortcutLabel() + "_vsum", "mz" );
258 : }
259 :
260 : // Now calculate the total length of the vector
261 2 : createVectorNormInput( getShortcutLabel(), getShortcutLabel() + "_norm", "mn" );
262 2 : multicolvar::MultiColvarShortcuts::expandFunctions( getShortcutLabel(), getShortcutLabel() + "_norm", "", this );
263 1 : }
264 :
265 1 : void HexacticParameter::createVectorNormInput( const std::string& ilab, const std::string& olab, const std::string& vlab ) {
266 2 : readInputLine( olab + "2: COMBINE PERIODIC=NO ARG=" + ilab + "_r" + vlab + "," + ilab + "_i" + vlab + " POWERS=2,2" );
267 2 : readInputLine( olab + ": CUSTOM ARG=" + olab + "2 FUNC=sqrt(x) PERIODIC=NO");
268 1 : }
269 :
270 : }
271 : }
272 :
|
