Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2011-2023 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/ActionWithValue.h"
25 : #include "core/PlumedMain.h"
26 : #include "core/ActionSet.h"
27 : #include "multicolvar/MultiColvarShortcuts.h"
28 :
29 : //+PLUMEDOC MCOLVAR SMAC
30 : /*
31 : Calculate the SMAC order parameter for a set of molecules
32 :
33 : This shortcut action provides a quantity that can be thought of as a [symmetry function](https://www.plumed-tutorials.org/lessons/23/001/data/SymmetryFunction.html) for molecules.
34 : You can thus use it to detect whether the molecules in a sample are ordered in the way that they would expect to be ordered in the crystal. If you are using this CV you normally use
35 : the [COM](COM.md) command to define centers of mass for your molecules and the [DISTANCE](DISTANCE.md) or [PLANE](PLANE.md) command to define each molecules orientation. You can thus calculate
36 : the relative orientation molecules $i$ and $j$ by calculating the torsional angle, $\theta_{ij}$, between the two molecular orientations around the vector $r_{ij}$ connecting the centers of
37 : mass of the two molecules. The value of the molecular order parameter for molecule $i$ is thus calculated as:
38 :
39 : $$
40 : s_i = frac{1-\gamma(c_i)}{c_i} \sum_j \sigma(|r_{ij}|) \sum_k K\left( \frac{\theta_{ij} - \phi_k}{b} \right) \qquad \textrm{where} \qquad c_i = \sum_j \sigma(|r_{ij}|)
41 : $$
42 :
43 : In this expression $\sigma$ is a switching function that acts on the distance between the centers of mass of molecules $i$ and $j$. $c_i$ is thus the number of molecules that are within
44 : a certain cutoff of molecule $i$ and $\gamma$ is another switching function that acts upon this quantity. This switching function ensures that the symmetry function is zero for atoms that are
45 : regions where the density of molecules are low. $K$ is then a kernel function with bandwidth $b$ that is centered at $\phi_k$. The function above is thus only large if molecule $i$ is surrounded
46 : by molecules whose relative orientations are as the user has requested by specifying $\phi_k$ parameters.
47 :
48 : The following example illustrates how the SMAC order parameter in PLUMED is used:
49 :
50 : ```plumed
51 : m3: DISTANCES ...
52 : ATOMS1=9,10 LOCATION1=9
53 : ATOMS2=89,90 LOCATION2=89
54 : ATOMS3=473,474 LOCATION3=473
55 : ATOMS4=1161,1162 LOCATION4=1161
56 : ATOMS5=1521,1522 LOCATION5=1521
57 : ATOMS6=1593,1594 LOCATION6=1593
58 : ATOMS7=1601,1602 LOCATION7=1601
59 : ATOMS8=2201,2202 LOCATION8=2201
60 : COMPONENTS
61 : ...
62 :
63 : s2: SMAC SPECIES=m3 KERNEL1={GAUSSIAN CENTER=0 SIGMA=0.480} KERNEL2={GAUSSIAN CENTER=pi SIGMA=0.480} SWITCH={RATIONAL R_0=0.6} MORE_THAN={RATIONAL R_0=0.7} SWITCH_COORD={EXP R_0=4}
64 :
65 : PRINT ARG=s2_morethan FILE=colvar
66 : ```
67 :
68 : Here the orientations of the molecules are specified by calculating the vectors connecting pairs of atoms in the molecules. The LOCATION keywords in the distance command are used to specify the
69 : positions of the molecules from which the $r_{ij}$ vectors in the above expression are calculated. The sum over $k$ in the above expression has two terms corresponding to the two Gaussian kernels
70 : that have been specified using KERNEL keywords. The SWITCH keyword has been used to specify the parameters of the switching function $\sigma$, while the SWITCH_COORD keyword has been used to specify
71 : the parameters of the switching function $\gamma$.
72 :
73 : A vector of 8 values for $s_i$ is calculated. The elements of this vector are transformed by a switching function that is one if the $s_i$ value is larger than 0.7. The eight elements of the resulting vector
74 : of transformed $s_i$ are then added together to determine the final quantity that is output to the colvar file.
75 :
76 : Incidentally, the authors who designed the SMAC symmetry function have forgotten what the letters in this acronym stand for.
77 :
78 : ## Working with two types of molecule
79 :
80 : You can use two different sets of vectors in the input to this action as shown below:
81 :
82 : ```plumed
83 : m3: DISTANCES ...
84 : ATOMS1=9,10 LOCATION1=9
85 : ATOMS2=89,90 LOCATION2=89
86 : ATOMS3=473,474 LOCATION3=473
87 : ATOMS4=1161,1162 LOCATION4=1161
88 : COMPONENTS
89 : ...
90 :
91 : m4: DISTANCES ...
92 : ATOMS1=1521,1522 LOCATION1=1521
93 : ATOMS2=1593,1594 LOCATION2=1593
94 : ATOMS3=1601,1602 LOCATION3=1601
95 : ATOMS4=2201,2202 LOCATION4=2201
96 : COMPONENTS
97 : ...
98 :
99 : s2: SMAC SPECIESA=m3 SPECIESB=m4 KERNEL1={GAUSSIAN CENTER=0 SIGMA=0.480} KERNEL2={GAUSSIAN CENTER=pi SIGMA=0.480} SWITCH={RATIONAL R_0=0.6} MORE_THAN={RATIONAL R_0=0.7} SWITCH_COORD={EXP R_0=4}
100 :
101 : PRINT ARG=s2_morethan FILE=colvar
102 : ```
103 :
104 : In this input the $s_i$ values defined using the equation above is evaluated for the four vectors that are specified in the DISTANCE action with label m3. In evaluating the sum in that expression
105 : we run over the four vectors that are specified in the action with label `m4` above.
106 :
107 :
108 : */
109 : //+ENDPLUMEDOC
110 :
111 : namespace PLMD {
112 : namespace symfunc {
113 :
114 : class SMAC : public ActionShortcut {
115 : public:
116 : static void registerKeywords(Keywords& keys);
117 : explicit SMAC(const ActionOptions&);
118 : };
119 :
120 : PLUMED_REGISTER_ACTION(SMAC,"SMAC")
121 :
122 15 : void SMAC::registerKeywords(Keywords& keys) {
123 15 : ActionShortcut::registerKeywords( keys );
124 15 : keys.add("optional","SPECIES","the label of the DISTANCES or PLANES action that computes the orientations of the molecules for which you would like to compute SMAC. The coordination sphere contains the same list of molecules if you use this keyword.");
125 15 : keys.add("optional","SPECIESA","the label of the DISTANCES or PLANES action that computes the orientations of the molecules for which you would like to compute SMAC. This keyword must be used with SPECIESB.");
126 15 : keys.add("optional","SPECIESB","the label of the DISTANCES or PLANES action that computes the orientations of the molecules that should be considered as part of the coordination sphere. This keyword must be used with SPECIESA.");
127 15 : keys.add("optional","SWITCH","This keyword is used if you want to employ an alternative to the continuous switching function defined above. "
128 : "The following provides information on the switching functions that are available. "
129 : "When this keyword is present you no longer need the NN, MM, D_0 and R_0 keywords.");
130 30 : keys.linkActionInDocs("SWITCH","LESS_THAN");
131 15 : keys.add("numbered","KERNEL","The kernels used in the function of the angle");
132 15 : keys.add("optional","SWITCH_COORD","This keyword is used to define the coordination switching function.");
133 30 : keys.linkActionInDocs("SWITCH_COORD","LESS_THAN");
134 30 : keys.reset_style("KERNEL","optional");
135 30 : keys.setValueDescription("vector","the value of the smac parameter for each of the input molecules");
136 15 : multicolvar::MultiColvarShortcuts::shortcutKeywords( keys );
137 15 : keys.needsAction("VSTACK");
138 15 : keys.needsAction("TRANSPOSE");
139 15 : keys.needsAction("CONTACT_MATRIX");
140 15 : keys.needsAction("TORSIONS_MATRIX");
141 15 : keys.needsAction("COMBINE");
142 15 : keys.needsAction("CUSTOM");
143 15 : keys.needsAction("ONES");
144 15 : keys.needsAction("MATRIX_VECTOR_PRODUCT");
145 15 : keys.needsAction("MORE_THAN");
146 15 : keys.addDOI("10.1016/j.ces.2014.08.032");
147 15 : keys.addDOI("10.1021/acs.jctc.6b01073");
148 15 : }
149 :
150 5 : SMAC::SMAC(const ActionOptions& ao):
151 : Action(ao),
152 5 : ActionShortcut(ao) {
153 : // Create the matrices
154 : std::string sw_input;
155 10 : parse("SWITCH",sw_input);
156 : std::string sp_lab, sp_laba;
157 5 : parse("SPECIES",sp_lab);
158 10 : parse("SPECIESA",sp_laba);
159 5 : if( sp_lab.length()>0 ) {
160 16 : readInputLine( getShortcutLabel() + "_vecs: VSTACK ARG=" + sp_lab + ".x," + sp_lab + ".y," + sp_lab + ".z" );
161 8 : readInputLine( getShortcutLabel() + "_vecsT: TRANSPOSE ARG=" + getShortcutLabel() + "_vecs" );
162 12 : readInputLine( getShortcutLabel() + "_cmap: CONTACT_MATRIX GROUP=" + sp_lab + " SWITCH={" + sw_input + "}");
163 28 : readInputLine( getShortcutLabel() + "_tpmat: TORSIONS_MATRIX ARG=" + getShortcutLabel() + "_vecs," + getShortcutLabel() + "_vecsT POSITIONS1=" + sp_lab + " POSITIONS2=" + sp_lab + " MASK=" + getShortcutLabel() + "_cmap");
164 1 : } else if( sp_laba.length()>0 ) {
165 : std::string sp_labb;
166 1 : parse("SPECIESB",sp_labb);
167 4 : readInputLine( getShortcutLabel() + "_vecsa: VSTACK ARG=" + sp_laba + ".x," + sp_laba + ".y," + sp_laba + ".z" );
168 4 : readInputLine( getShortcutLabel() + "_vecsb: VSTACK ARG=" + sp_labb + ".x," + sp_labb + ".y," + sp_labb + ".z" );
169 2 : readInputLine( getShortcutLabel() + "_vecsbT: TRANSPOSE ARG=" + getShortcutLabel() + "_vecsb" );
170 4 : readInputLine( getShortcutLabel() + "_cmap: CONTACT_MATRIX GROUPA=" + sp_laba + " GROUPB=" + sp_labb + " SWITCH={" + sw_input + "}");
171 7 : readInputLine( getShortcutLabel() + "_tpmat: TORSIONS_MATRIX ARG=" + getShortcutLabel() + "_vecsa," + getShortcutLabel() + "_vecsbT POSITIONS1=" + sp_laba + " POSITIONS2=" + sp_labb + " MASK=" + getShortcutLabel() + "_cmap");
172 : }
173 : // Now need the Gaussians
174 5 : std::string kmap_input= getShortcutLabel() + "_ksum: COMBINE MASK=" + getShortcutLabel() + "_cmap PERIODIC=NO";
175 5 : for(unsigned i=1;; ++i) {
176 : std::string kstr_inpt, istr;
177 15 : Tools::convert( i, istr );
178 30 : if( !parseNumbered("KERNEL",i,kstr_inpt ) ) {
179 : break;
180 : }
181 10 : std::vector<std::string> words = Tools::getWords(kstr_inpt);
182 : std::string center, var;
183 10 : Tools::parse(words,"CENTER",center);
184 10 : Tools::parse(words,"SIGMA",var);
185 : double nsig;
186 10 : Tools::convert( var, nsig );
187 : std::string coeff;
188 10 : Tools::convert( 1/(nsig*nsig), coeff );
189 70 : readInputLine( getShortcutLabel() + "_kf" + istr + "_r2: COMBINE MASK=" + getShortcutLabel() + "_cmap PERIODIC=NO ARG=" + getShortcutLabel() + "_tpmat COEFFICIENTS=" + coeff + " PARAMETERS=" + center + " POWERS=2");
190 10 : if( words[0]=="GAUSSIAN" ) {
191 48 : readInputLine( getShortcutLabel() + "_kf" + istr + ": CUSTOM MASK=" + getShortcutLabel() + "_cmap PERIODIC=NO FUNC=exp(-x/2) ARG=" + getShortcutLabel() + "_kf" + istr + "_r2" );
192 2 : } else if( words[0]=="TRIANGULAR" ) {
193 12 : readInputLine( getShortcutLabel() + "_kf" + istr + ": CUSTOM MASK=" + getShortcutLabel() + "_cmap PERIODIC=NO FUNC=step(1-sqrt(x))*(1-sqrt(x)) ARG=" + getShortcutLabel() + "_kf" + istr + "_r2" );
194 : } else {
195 0 : readInputLine( getShortcutLabel() + "_kf" + istr + ": CUSTOM PERIODIC=NO MASK=" + getShortcutLabel() + "_cmap FUNC=" + words[0] + " ARG=" + getShortcutLabel() + "_kf" + istr + "_r2" );
196 : }
197 10 : if( i==1 ) {
198 15 : kmap_input += " ARG=" + getShortcutLabel() + "_kf" + istr;
199 : } else {
200 15 : kmap_input += "," + getShortcutLabel() + "_kf" + istr;
201 : }
202 20 : }
203 5 : readInputLine( kmap_input );
204 : // Now create the product matrix
205 15 : readInputLine( getShortcutLabel() + "_prod: CUSTOM ARG=" + getShortcutLabel() + "_cmap," + getShortcutLabel() + "_ksum FUNC=x*y PERIODIC=NO");
206 : // Now the sum of coordination numbers times the switching functions
207 5 : ActionWithValue* av = plumed.getActionSet().selectWithLabel<ActionWithValue*>( getShortcutLabel() + "_cmap");
208 5 : plumed_assert( av && av->getNumberOfComponents()>0 && (av->copyOutput(0))->getRank()==2 );
209 : std::string size;
210 5 : Tools::convert( (av->copyOutput(0))->getShape()[1], size );
211 5 : readInputLine( getShortcutLabel() + "_ones: ONES SIZE=" + size );
212 15 : readInputLine( getShortcutLabel() + ": MATRIX_VECTOR_PRODUCT ARG=" + getShortcutLabel() + "_prod," + getShortcutLabel() + "_ones");
213 : // And just the sum of the coordination numbers
214 20 : readInputLine( getShortcutLabel() + "_denom: MATRIX_VECTOR_PRODUCT ARG=" + getShortcutLabel() + "_cmap," + getShortcutLabel() + "_ones");
215 : // And the transformed switching functions
216 : std::string swcoord_str;
217 5 : parse("SWITCH_COORD",swcoord_str);
218 15 : readInputLine( getShortcutLabel() + "_mtdenom: MORE_THAN ARG=" + getShortcutLabel() + "_denom SWITCH={" + swcoord_str +"}");
219 : // And matheval to get the final quantity
220 20 : readInputLine( getShortcutLabel() + "_smac: CUSTOM ARG=" + getShortcutLabel() + "," + getShortcutLabel() + "_mtdenom," + getShortcutLabel() + "_denom FUNC=(x*y)/z PERIODIC=NO");
221 : // And this expands everything
222 10 : multicolvar::MultiColvarShortcuts::expandFunctions( getShortcutLabel(), getShortcutLabel() + "_smac", "", this );
223 5 : }
224 :
225 : }
226 : }
|
