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 SWITCH={RATIONAL D_0=1.4 R_0=0.2} MEAN
62 : PRINT ARG=hex.mean FILE=colvar
63 : ```
64 :
65 : 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.
66 : The elements of this vector are then added together so the mean value can be computed.
67 :
68 : 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
69 : six nearest neighbours to each atom. If you would like to calculate this quantity using PLUMED you can use an input like this:
70 :
71 : ```plumed
72 : dmat: DISTANCE_MATRIX GROUP=1-400 CUTOFF=3.0 COMPONENTS
73 : neigh: NEIGHBORS ARG=dmat.w NLOWEST=6
74 : harm: CYLINDRICAL_HARMONIC DEGREE=6 ARG=dmat.x,dmat.y
75 : rprod: CUSTOM ARG=neigh,harm.rm FUNC=x*y PERIODIC=NO
76 : iprod: CUSTOM ARG=neigh,harm.im FUNC=x*y PERIODIC=NO
77 : hex2_ones: ONES SIZE=400
78 : hex2_denom: MATRIX_VECTOR_PRODUCT ARG=neigh,hex2_ones
79 : harm_rm: MATRIX_VECTOR_PRODUCT ARG=rprod,hex2_ones
80 : harm_im: MATRIX_VECTOR_PRODUCT ARG=iprod,hex2_ones
81 : hex2_rmn: CUSTOM ARG=harm_rm,hex2_denom FUNC=x/y PERIODIC=NO
82 : hex2_imn: CUSTOM ARG=harm_im,hex2_denom FUNC=x/y PERIODIC=NO
83 : DUMPATOMS ATOMS=1-400 ARG=hex2_rmn,hex2_imn,hex2_denom FILE=hexparam.xyz
84 : ```
85 :
86 : This input outputs the values of the order parameters for all the atoms to an extended xyz file .
87 :
88 : > ![CAUTION]
89 : > Virial is not working currently
90 :
91 :
92 : */
93 : //+ENDPLUMEDOC
94 :
95 :
96 : class HexacticParameter : public ActionShortcut {
97 : private:
98 : void createVectorNormInput( const std::string& ilab, const std::string& olab, const std::string& vlab );
99 : public:
100 : static void registerKeywords( Keywords& keys );
101 : explicit HexacticParameter(const ActionOptions&);
102 : };
103 :
104 : PLUMED_REGISTER_ACTION(HexacticParameter,"HEXACTIC_PARAMETER")
105 :
106 5 : void HexacticParameter::registerKeywords( Keywords& keys ) {
107 5 : CoordinationNumbers::shortcutKeywords( keys );
108 5 : keys.add("compulsory","PLANE","the plane to use when calculating the value of the order parameter should be xy, xz or yz");
109 10 : keys.setValueDescription("matrix","the value of the cylindrical harmonic for each bond vector specified");
110 5 : keys.addFlag("VMEAN",false,"calculate the norm of the mean vector.");
111 10 : keys.addOutputComponent("_vmean","VMEAN","scalar","the norm of the mean vector");
112 5 : keys.addFlag("VSUM",false,"calculate the norm of the sum of all the vectors");
113 10 : keys.addOutputComponent("_vsum","VSUM","scalar","the norm of the mean vector");
114 5 : keys.needsAction("CYLINDRICAL_HARMONIC_MATRIX");
115 5 : keys.needsAction("ONES");
116 5 : keys.needsAction("MATRIX_VECTOR_PRODUCT");
117 5 : keys.needsAction("CUSTOM");
118 5 : keys.needsAction("MEAN");
119 5 : keys.needsAction("SUM");
120 5 : keys.needsAction("COMBINE");
121 5 : keys.addDOI("https://journals.aps.org/prl/abstract/10.1103/PhysRevLett.99.215701");
122 5 : }
123 :
124 1 : HexacticParameter::HexacticParameter( const ActionOptions& ao):
125 : Action(ao),
126 1 : ActionShortcut(ao) {
127 : std::string sp_str, specA, specB;
128 1 : parse("SPECIES",sp_str);
129 1 : parse("SPECIESA",specA);
130 1 : parse("SPECIESB",specB);
131 1 : CoordinationNumbers::expandMatrix( true, getShortcutLabel(), sp_str, specA, specB, this );
132 : std::string myplane;
133 2 : parse("PLANE",myplane);
134 1 : if( myplane=="xy" ) {
135 2 : readInputLine( getShortcutLabel() + ": CYLINDRICAL_HARMONIC_MATRIX DEGREE=6 ARG=" + getShortcutLabel() + "_mat.x," + getShortcutLabel() + "_mat.y," + getShortcutLabel() + "_mat.w" );
136 0 : } else if( myplane=="xz" ) {
137 0 : readInputLine( getShortcutLabel() + ": CYLINDRICAL_HARMONIC_MATRIX DEGREE=6 ARG=" + getShortcutLabel() + "_mat.x," + getShortcutLabel() + "_mat.z," + getShortcutLabel() + "_mat.w" );
138 0 : } else if( myplane=="yz" ) {
139 0 : readInputLine( getShortcutLabel() + ": CYLINDRICAL_HARMONIC_MATRIX DEGREE=6 ARG=" + getShortcutLabel() + "_mat.y," + getShortcutLabel() + "_mat.z," + getShortcutLabel() + "_mat.w" );
140 : } else {
141 0 : error("invalid input for plane -- should be xy, xz or yz");
142 : }
143 : // And coordination number
144 1 : ActionWithValue* av = plumed.getActionSet().selectWithLabel<ActionWithValue*>( getShortcutLabel() + "_mat");
145 1 : plumed_assert( av && av->getNumberOfComponents()>0 && (av->copyOutput(0))->getRank()==2 );
146 : std::string size;
147 1 : Tools::convert( (av->copyOutput(0))->getShape()[1], size );
148 2 : readInputLine( getShortcutLabel() + "_ones: ONES SIZE=" + size );
149 2 : readInputLine( getShortcutLabel() + "_rm: MATRIX_VECTOR_PRODUCT ARG=" + getShortcutLabel() + ".rm," + getShortcutLabel() + "_ones");
150 2 : readInputLine( getShortcutLabel() + "_im: MATRIX_VECTOR_PRODUCT ARG=" + getShortcutLabel() + ".im," + getShortcutLabel() + "_ones");
151 : // Input for denominator (coord)
152 2 : readInputLine( getShortcutLabel() + "_denom: MATRIX_VECTOR_PRODUCT ARG=" + getShortcutLabel() + "_mat.w," + getShortcutLabel() + "_ones");
153 : // Divide real part by coordination numbers
154 2 : readInputLine( getShortcutLabel() + "_rmn: CUSTOM ARG=" + getShortcutLabel() + "_rm," + getShortcutLabel() + "_denom FUNC=x/y PERIODIC=NO");
155 : // Devide imaginary part by coordination number
156 2 : readInputLine( getShortcutLabel() + "_imn: CUSTOM ARG=" + getShortcutLabel() + "_im," + getShortcutLabel() + "_denom FUNC=x/y PERIODIC=NO");
157 :
158 : // If we are doing VMEAN determine sum of vector components
159 : bool do_vmean;
160 1 : parseFlag("VMEAN",do_vmean);
161 1 : if( do_vmean ) {
162 : // Real part
163 0 : readInputLine( getShortcutLabel() + "_rms: MEAN ARG=" + getShortcutLabel() + "_rmn PERIODIC=NO");
164 : // Imaginary part
165 0 : readInputLine( getShortcutLabel() + "_ims: MEAN ARG=" + getShortcutLabel() + "_imn PERIODIC=NO");
166 : // Now calculate the total length of the vector
167 0 : createVectorNormInput( getShortcutLabel(), getShortcutLabel() + "_vmean", "ms" );
168 : }
169 : bool do_vsum;
170 1 : parseFlag("VSUM",do_vsum);
171 1 : if( do_vsum ) {
172 : // Real part
173 0 : readInputLine( getShortcutLabel() + "_rmz: SUM ARG=" + getShortcutLabel() + "_rmn PERIODIC=NO");
174 : // Imaginary part
175 0 : readInputLine( getShortcutLabel() + "_imz: SUM ARG=" + getShortcutLabel() + "_imn PERIODIC=NO");
176 : // Now calculate the total length of the vector
177 0 : createVectorNormInput( getShortcutLabel(), getShortcutLabel() + "_vsum", "mz" );
178 : }
179 :
180 : // Now calculate the total length of the vector
181 2 : createVectorNormInput( getShortcutLabel(), getShortcutLabel() + "_norm", "mn" );
182 2 : multicolvar::MultiColvarShortcuts::expandFunctions( getShortcutLabel(), getShortcutLabel() + "_norm", "", this );
183 1 : }
184 :
185 1 : void HexacticParameter::createVectorNormInput( const std::string& ilab, const std::string& olab, const std::string& vlab ) {
186 2 : readInputLine( olab + "2: COMBINE PERIODIC=NO ARG=" + ilab + "_r" + vlab + "," + ilab + "_i" + vlab + " POWERS=2,2" );
187 2 : readInputLine( olab + ": CUSTOM ARG=" + olab + "2 FUNC=sqrt(x) PERIODIC=NO");
188 1 : }
189 :
190 : }
191 : }
192 :
|
