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 "core/ActionRegister.h"
24 : #include "core/PlumedMain.h"
25 : #include "core/ActionSet.h"
26 : #include "core/ActionWithValue.h"
27 : #include "CoordinationNumbers.h"
28 : #include "multicolvar/MultiColvarShortcuts.h"
29 :
30 : #include <complex>
31 :
32 : namespace PLMD {
33 : namespace symfunc {
34 :
35 : //+PLUMEDOC MCOLVAR LOCAL_CRYSTALINITY
36 : /*
37 : Calculate the local crystalinity symmetry function
38 :
39 : This shortcut provides an implementation of the local crystalinity order parameter that is described in the paper in the bibliography.
40 : To use this CV you define a series of unit vectors, $g_k$, using multiple instances of the GVECTOR keyword. This allows you to define a
41 : [symmetry function](https://www.plumed-tutorials.org/lessons/23/001/data/SymmetryFunction.html) for the $i$th atom as:
42 :
43 : $$
44 : s_i = \sum_k \left| \frac{\sum_j \sigma(|r_{ij}|) e^{ig_k r_{ij}}}{\sum_j \sigma(|r_{ij}|)} \right|^2
45 : $$
46 :
47 : In this expression $r_{ij}$ is the vector connecting atom $i$ to atom $j$ and $\sigma$ is a switching function that acts upon the magnidue of this vector, $|r_{ij}|$.
48 : The following input is an example that shows how this symmetry function can be used in practice.
49 :
50 : ```plumed
51 : lc: LOCAL_CRYSTALINITY ...
52 : SPECIES=1-64 D_0=3.0 R_0=1.5
53 : GVECTOR1=1,1,1 GVECTOR2=1,0.5,0.5 GVECTOR3=0.5,1.0,1.0 SUM
54 : ...
55 : PRINT ARG=lc_sum FILE=colvar
56 : ```
57 :
58 : As you can see if you expand the shortcut in this input, the sum over $k$ in the above expression has three terms in this input as 3 GVECTORS are specified.
59 : Sixty four values for the expression above are computed. These sixty four numbers are then added together in order to give a global mesuare of the crystallinity
60 : for the simulated system.
61 :
62 : In the input above $\sigma(|r_{ij}|)$ is a rational [switching function](LESS_THAN.md)
63 : with the parameters above. We would recommend using SWITCH syntax
64 : rather than the syntax above when giving the parameters for the switching function as you can then use any of the switching functions described
65 : 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
66 : switching function as demonstrated below:
67 :
68 : ```plumed
69 : lc: LOCAL_CRYSTALINITY ...
70 : SPECIES=1-64 SWITCH={RATIONAL D_0=3.0 R_0=1.5 D_MAX=3.0}
71 : GVECTOR1=1,1,1 GVECTOR2=1,0.5,0.5 GVECTOR3=0.5,1.0,1.0 SUM
72 : ...
73 : PRINT ARG=lc_sum FILE=colvar
74 : ```
75 :
76 : ## Working with two types of atom
77 :
78 : If you would like to calculate the function above for the bonds connecting the atoms in GROUPA to the atoms in GROUPB you can use an input like the one
79 : shown below:
80 :
81 : ```plumed
82 : lc: LOCAL_CRYSTALINITY ...
83 : SPECIESA=1-5 SPECIESB=1-64 SWITCH={RATIONAL D_0=3.0 R_0=1.5 D_MAX=3.0}
84 : GVECTOR1=1,1,1 GVECTOR2=1,0.5,0.5 GVECTOR3=0.5,1.0,1.0 SUM
85 : ...
86 : PRINT ARG=lc_sum FILE=colvar
87 : ```
88 :
89 : ## The MASK keyword
90 :
91 : 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
92 : 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
93 : from the [volumes module](module_volumes.md) to calculate the average value of the TETRAHEDRAL parameter for only those atoms that lie in a certain part of the simulation box.
94 :
95 : ```plumed
96 : # Fixed virtual atom which serves as the probe volume's center (pos. in nm)
97 : center: FIXEDATOM AT=2.5,2.5,2.5
98 : # Vector in which element i is one if atom i is in sphere of interest and zero otherwise
99 : sphere: INSPHERE ATOMS=1-64 CENTER=center RADIUS={GAUSSIAN D_0=0.5 R_0=0.01 D_MAX=0.52}
100 : # Calculate the local crystallinity parameters
101 : lc: LOCAL_CRYSTALINITY ...
102 : SPECIES=1-64 SWITCH={RATIONAL D_0=3.0 R_0=1.5 D_MAX=3.0}
103 : GVECTOR1=1,1,1 GVECTOR2=1,0.5,0.5 GVECTOR3=0.5,1.0,1.0 MASK=sphere
104 : ...
105 : # Multiply local crystallinity parameters numbers by sphere vector
106 : prod: CUSTOM ARG=lc,sphere FUNC=x*y PERIODIC=NO
107 : # Sum of coordination numbers for atoms that are in the sphere of interest
108 : numer: SUM ARG=prod PERIODIC=NO
109 : # Number of atoms that are in sphere of interest
110 : denom: SUM ARG=sphere PERIODIC=NO
111 : # Average coordination number for atoms in sphere of interest
112 : av: CUSTOM ARG=numer,denom FUNC=x/y PERIODIC=NO
113 : # And print out final CV to a file
114 : PRINT ARG=av FILE=colvar STRIDE=1
115 : ```
116 :
117 : This input calculate the average local crystallinity parameter for only those atoms that are within a spherical region that is centered on the point
118 : $(2.5,2.5,2.5)$.
119 :
120 : */
121 : //+ENDPLUMEDOC
122 :
123 :
124 : class LocalCrystallinity : public ActionShortcut {
125 : public:
126 : static void registerKeywords( Keywords& keys );
127 : explicit LocalCrystallinity(const ActionOptions&);
128 : };
129 :
130 : PLUMED_REGISTER_ACTION(LocalCrystallinity,"LOCAL_CRYSTALINITY")
131 :
132 6 : void LocalCrystallinity::registerKeywords( Keywords& keys ) {
133 6 : CoordinationNumbers::shortcutKeywords( keys );
134 6 : keys.add("numbered","GVECTOR","the coefficients of the linear combinations to compute for the CV");
135 12 : keys.setValueDescription("vector","the value of the local crystalinity for each of the input atoms");
136 6 : keys.needsAction("ONES");
137 6 : keys.needsAction("MATRIX_VECTOR_PRODUCT");
138 6 : keys.needsAction("COMBINE");
139 6 : keys.needsAction("CUSTOM");
140 6 : keys.addDOI("10.1063/1.4822877");
141 6 : }
142 :
143 2 : LocalCrystallinity::LocalCrystallinity( const ActionOptions& ao):
144 : Action(ao),
145 2 : ActionShortcut(ao) {
146 : // This builds an adjacency matrix
147 : std::string sp_str, specA, specB;
148 2 : parse("SPECIES",sp_str);
149 2 : parse("SPECIESA",specA);
150 2 : parse("SPECIESB",specB);
151 2 : CoordinationNumbers::expandMatrix( true, getShortcutLabel(), sp_str, specA, specB, this );
152 : // Input for denominator (coord)
153 2 : ActionWithValue* av = plumed.getActionSet().selectWithLabel<ActionWithValue*>( getShortcutLabel() + "_mat");
154 3 : plumed_assert( av && av->getNumberOfComponents()>0 && (av->copyOutput(0))->getRank()==2 );
155 : std::string size;
156 2 : Tools::convert( (av->copyOutput(0))->getShape()[1], size );
157 2 : readInputLine( getShortcutLabel() + "_ones: ONES SIZE=" + size );
158 6 : readInputLine( getShortcutLabel() + "_denom: MATRIX_VECTOR_PRODUCT ARG=" + getShortcutLabel() + "_mat.w," + getShortcutLabel() + "_ones");
159 : // Input for numerator
160 3 : std::string finput = "";
161 2 : for(unsigned i=1;; ++i) {
162 5 : std::vector<std::string> gvec;
163 : std::string istr;
164 5 : Tools::convert( i, istr );
165 10 : if( !parseNumberedVector("GVECTOR",i,gvec) ) {
166 : break;
167 : }
168 3 : if( gvec.size()!=3 ) {
169 0 : error("gvectors should have size 3");
170 : }
171 : // This is the dot product between the input gvector and the bond vector
172 15 : readInputLine( getShortcutLabel() + "_dot" + istr + ": COMBINE ARG=" + getShortcutLabel() + "_mat.x," + getShortcutLabel() + "_mat.y," + getShortcutLabel() + "_mat.z "
173 9 : "PERIODIC=NO COEFFICIENTS=" + gvec[0] + "," + gvec[1] + "," + gvec[2] );
174 : // Now calculate the sine and cosine of the dot product
175 15 : readInputLine( getShortcutLabel() + "_cos" + istr + ": CUSTOM ARG=" + getShortcutLabel() +"_mat.w," + getShortcutLabel() + "_dot" + istr + " FUNC=x*cos(y) PERIODIC=NO");
176 15 : readInputLine( getShortcutLabel() + "_sin" + istr + ": CUSTOM ARG=" + getShortcutLabel() +"_mat.w," + getShortcutLabel() + "_dot" + istr + " FUNC=x*sin(y) PERIODIC=NO");
177 : // And sum up the sine and cosine over the coordination spheres
178 15 : readInputLine( getShortcutLabel() + "_cossum" + istr + ": MATRIX_VECTOR_PRODUCT ARG=" + getShortcutLabel() + "_cos" + istr + "," + getShortcutLabel() + "_ones");
179 15 : readInputLine( getShortcutLabel() + "_sinsum" + istr + ": MATRIX_VECTOR_PRODUCT ARG=" + getShortcutLabel() + "_sin" + istr + "," + getShortcutLabel() + "_ones");
180 : // And average the sine and cosine over the number of bonds
181 15 : readInputLine( getShortcutLabel() + "_cosmean" + istr + ": CUSTOM FUNC=x/y PERIODIC=NO ARG=" + getShortcutLabel() + "_cossum" + istr + "," + getShortcutLabel() + "_denom");
182 15 : readInputLine( getShortcutLabel() + "_sinmean" + istr + ": CUSTOM FUNC=x/y PERIODIC=NO ARG=" + getShortcutLabel() + "_sinsum" + istr + "," + getShortcutLabel() + "_denom");
183 : // And work out the square modulus of this complex number
184 15 : readInputLine( getShortcutLabel() + "_" + istr + ": CUSTOM FUNC=x*x+y*y PERIODIC=NO ARG=" + getShortcutLabel() + "_cosmean" + istr + "," + getShortcutLabel() + "_sinmean" + istr);
185 : // These are all the kvectors that we are adding together in the final combine for the final CV
186 3 : if( i>1 ) {
187 : finput += ",";
188 : }
189 6 : finput += getShortcutLabel() + "_" + istr;
190 5 : }
191 2 : if( finput.empty() ) {
192 1 : error("no GVECTOR was specified: the numbered GVECTOR keywords start from GVECTOR1");
193 : }
194 : // This computes the final CV
195 1 : readInputLine( getShortcutLabel() + ": COMBINE NORMALIZE PERIODIC=NO ARG=" + finput );
196 : // Now calculate the total length of the vector
197 3 : multicolvar::MultiColvarShortcuts::expandFunctions( getShortcutLabel(), getShortcutLabel(), "", this );
198 3 : }
199 :
200 : }
201 : }
202 :
|
