Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2012-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 "CoordinationNumbers.h"
23 : #include "multicolvar/MultiColvarShortcuts.h"
24 : #include "core/ActionWithValue.h"
25 : #include "core/ActionRegister.h"
26 : #include "core/PlumedMain.h"
27 : #include "core/ActionSet.h"
28 :
29 : #include <string>
30 : #include <cmath>
31 :
32 : namespace PLMD {
33 : namespace symfunc {
34 :
35 : //+PLUMEDOC MCOLVAR COORDINATIONNUMBER
36 : /*
37 : Calculate the coordination numbers of atoms so that you can then calculate functions of the distribution of
38 : coordination numbers such as the minimum, the number less than a certain quantity and so on.
39 :
40 : The coordination number of a atom $i$ is the number of atoms that are within a certain cutoff distance of it.
41 : This quantity can be calculated as follows:
42 :
43 : $$
44 : s_i = \sum_j \sigma(r_{ij})
45 : $$
46 :
47 : where $r_{ij}$ is the distance between atoms $i$ and $j$ and $\sigma$ is a switching function. The typical switching
48 : function that is used in metadynamics is this one:
49 :
50 : $$
51 : s(r) = \frac{ 1 - \left(\frac{r-d_0}{r_0}\right)^n } { 1 - \left(\frac{r-d_0}{r_0}\right)^m }
52 : $$
53 :
54 : The following example shows how you can use this shortcut action to calculate and print the coordination numbers of
55 : one hundred atoms with themselves:
56 :
57 : ```plumed
58 : c: COORDINATIONNUMBER SPECIES=1-100 R_0=1.0
59 : DUMPATOMS ATOMS=c ARG=c FILE=coords.xyz
60 : ```
61 :
62 : This input will produce an output file called coords that contains the coordination numbers of the 100 input atoms. The cutoff
63 : that is used to calculate the coordination number in this case is 1.0.
64 :
65 : The vectors that are output by the COORDINATIONNUMBER shortcut can be used in the input for many other functions that are within
66 : PLUMED. In addition, in order to ensure compatibility with older versions of PLUMED you can add additional keywords on the input
67 : line for COORDINATIONNUMBER in order to calculate various functions of the input. For example, the following input tells plumed
68 : ato calculate the coordination numbers of atoms 1-100 with themselves. The minimum coordination number is then calculated.
69 :
70 : ```plumed
71 : c: COORDINATIONNUMBER SPECIES=1-100 R_0=1.0 MIN={BETA=0.1}
72 : ```
73 :
74 : By constrast, this input tells plumed to calculate how many atoms from 1-100 are within 3.0 of each of the atoms
75 : from 101-110. In the first 101 is the central atom, in the second 102 is the central atom and so on. The
76 : number of coordination numbers that are more than 6 is then computed.
77 :
78 : ```plumed
79 : c: COORDINATIONNUMBER SPECIESA=101-110 SPECIESB=1-100 R_0=3.0 MORE_THAN={RATIONAL R_0=6.0 NN=6 MM=12 D_0=0}
80 : ```
81 :
82 : Notice that these inputs both use shortcuts. If you expand the inputs above you can determine the set of actions
83 : that are being used to calculate each of the quantities of interest.
84 :
85 : */
86 : //+ENDPLUMEDOC
87 :
88 : //+PLUMEDOC MCOLVAR COORDINATION_MOMENTS
89 : /*
90 : Calculate moments of the distribution of distances in the first coordination sphere
91 :
92 : This is the CV that was developed by White and Voth and is described in the paper in the bibliograhy below. This action provides a way of indirectly biasing radial distribution functions and computes the following function
93 :
94 : $$
95 : s_i = \sum_j \sigma(r_{ij})r_{ij}^k
96 : $$
97 :
98 : where $k$ is the value that is input using the R_POWER keyword, $r_{ij}$ is the distance between atoms $i$ and $j$ and $\sigma$ is a switching function.
99 :
100 : The following example shows how this action can be used.
101 :
102 : ```plumed
103 : cn0: COORDINATIONNUMBER SPECIES=1-10 SWITCH={RATIONAL R_0=1.0 D_MAX=8} MEAN
104 : cn1: COORDINATION_MOMENTS SPECIES=1-10 SWITCH={RATIONAL R_0=1.0 D_MAX=8} R_POWER=1 MEAN
105 : cn2: COORDINATION_MOMENTS SPECIES=1-10 SWITCH={RATIONAL R_0=1.0 D_MAX=8} R_POWER=2 MEAN
106 : PRINT ARG=cn0.mean,cn1.mean,cn2.mean STRIDE=1 FILE=cn_out
107 : ```
108 :
109 : As you can see, it works similarlly to [COORDINATIONNUMBER](COORDINATIONNUMBER.md).
110 :
111 :
112 : */
113 : //+ENDPLUMEDOC
114 :
115 :
116 : PLUMED_REGISTER_ACTION(CoordinationNumbers,"COORDINATIONNUMBER")
117 : PLUMED_REGISTER_ACTION(CoordinationNumbers,"COORDINATION_MOMENTS")
118 :
119 176 : void CoordinationNumbers::shortcutKeywords( Keywords& keys ) {
120 176 : ActionShortcut::registerKeywords( keys );
121 176 : keys.add("atoms-3","SPECIES","the list of atoms for which the symmetry function is being calculated and the atoms that can be in the environments");
122 176 : keys.add("atoms-4","SPECIESA","the list of atoms for which the symmetry function is being calculated. This keyword must be used in conjunction with SPECIESB, which specifies the atoms that are in the environment.");
123 176 : keys.add("atoms-4","SPECIESB","the list of atoms that can be in the environments of each of the atoms for which the symmetry function is being calculated. This keyword must be used in conjunction with SPECIESA, which specifies the atoms for which the symmetry function is being calculated.");
124 176 : keys.add("compulsory","NN","6","The n parameter of the switching function ");
125 176 : keys.add("compulsory","MM","0","The m parameter of the switching function; 0 implies 2*NN");
126 176 : keys.add("compulsory","D_0","0.0","The d_0 parameter of the switching function");
127 176 : keys.add("compulsory","R_0","The r_0 parameter of the switching function");
128 176 : keys.add("optional","SWITCH","the switching function that it used in the construction of the contact matrix");
129 352 : keys.linkActionInDocs("SWITCH","LESS_THAN");
130 176 : multicolvar::MultiColvarShortcuts::shortcutKeywords( keys );
131 176 : keys.needsAction("CONTACT_MATRIX");
132 176 : keys.needsAction("GROUP");
133 176 : keys.addDOI("10.1021/ct500320c");
134 176 : }
135 :
136 108 : void CoordinationNumbers::expandMatrix( const bool& components, const std::string& lab, const std::string& sp_str,
137 : const std::string& spa_str, const std::string& spb_str, ActionShortcut* action ) {
138 108 : if( sp_str.length()==0 && spa_str.length()==0 ) {
139 0 : return;
140 : }
141 :
142 108 : std::string matinp = lab + "_mat: CONTACT_MATRIX";
143 108 : if( sp_str.length()>0 ) {
144 69 : matinp += " GROUP=" + sp_str;
145 138 : action->readInputLine( lab + "_grp: GROUP ATOMS=" + sp_str );
146 39 : } else if( spa_str.length()>0 ) {
147 78 : matinp += " GROUPA=" + spa_str + " GROUPB=" + spb_str;
148 78 : action->readInputLine( lab + "_grp: GROUP ATOMS=" + spa_str );
149 : }
150 :
151 : std::string sw_str;
152 216 : action->parse("SWITCH",sw_str);
153 108 : if( sw_str.length()>0 ) {
154 162 : matinp += " SWITCH={" + sw_str + "}";
155 : } else {
156 : std::string r0;
157 54 : action->parse("R_0",r0);
158 : std::string d0;
159 54 : action->parse("D_0",d0);
160 27 : if( r0.length()==0 ) {
161 0 : action->error("missing switching function parameters use SWITCH/R_0");
162 : }
163 : std::string nn;
164 54 : action->parse("NN",nn);
165 : std::string mm;
166 27 : action->parse("MM",mm);
167 54 : matinp += " R_0=" + r0 + " D_0=" + d0 + " NN=" + nn + " MM=" + mm;
168 : }
169 108 : if( components ) {
170 : matinp += " COMPONENTS";
171 : }
172 108 : action->readInputLine( matinp );
173 : }
174 :
175 67 : void CoordinationNumbers::registerKeywords( Keywords& keys ) {
176 67 : shortcutKeywords( keys );
177 67 : keys.add("compulsory","R_POWER","the power to which you want to raise the distance");
178 67 : keys.addFlag("LOWMEM",false,"this flag does nothing and is present only to ensure back-compatibility");
179 67 : keys.add("optional","MOMENTS","the list of moments that you would like to calculate");
180 134 : keys.addOutputComponent("moment","MOMENTS","scalar","the moments of the distribution");
181 67 : keys.needsAction("MATRIX_VECTOR_PRODUCT");
182 67 : keys.needsAction("ONES");
183 67 : keys.needsAction("MOMENTS");
184 134 : keys.setValueDescription("vector","the coordination numbers of the specified atoms");
185 67 : }
186 :
187 45 : CoordinationNumbers::CoordinationNumbers(const ActionOptions& ao):
188 : Action(ao),
189 45 : ActionShortcut(ao) {
190 : bool lowmem;
191 45 : parseFlag("LOWMEM",lowmem);
192 45 : if( lowmem ) {
193 0 : warning("LOWMEM flag is deprecated and is no longer required for this action");
194 : }
195 : // Setup the contract matrix if that is what is needed
196 : std::string matlab, sp_str, specA, specB;
197 45 : parse("SPECIES",sp_str);
198 45 : parse("SPECIESA",specA);
199 90 : parse("SPECIESB",specB);
200 45 : if( sp_str.length()>0 || specA.length()>0 ) {
201 45 : matlab = getShortcutLabel() + "_mat";
202 45 : bool comp=false;
203 45 : if( getName()=="COORDINATION_MOMENTS" ) {
204 3 : comp=true;
205 6 : matlab = getShortcutLabel() + "_mat";
206 : }
207 45 : expandMatrix( comp, getShortcutLabel(), sp_str, specA, specB, this );
208 : } else {
209 0 : error("missing atoms input use SPECIES or SPECIESA/SPECIESB");
210 : }
211 45 : ActionWithValue* mb=plumed.getActionSet().selectWithLabel<ActionWithValue*>( matlab );
212 45 : if( !mb ) {
213 0 : error("could not find action with name " + matlab );
214 : }
215 45 : Value* arg=mb->copyOutput(0);
216 45 : if( arg->getRank()!=2 || arg->hasDerivatives() ) {
217 0 : error("the input to this action should be a matrix or scalar");
218 : }
219 : // Create vector of ones to multiply input matrix by
220 : std::string nones;
221 45 : Tools::convert( arg->getShape()[1], nones );
222 90 : readInputLine( getShortcutLabel() + "_ones: ONES SIZE=" + nones );
223 45 : if( getName()=="COORDINATION_MOMENTS" ) {
224 : // Calculate the lengths of the vectors
225 : std::string r_power;
226 3 : parse("R_POWER",r_power);
227 9 : readInputLine( getShortcutLabel() + "_pow: CUSTOM ARG=" + matlab + ".x," + matlab + ".y," + matlab + ".z," + matlab + ".w VAR=x,y,z,w "
228 6 : + "PERIODIC=NO FUNC=w*(sqrt(x*x+y*y+z*z)^" + r_power +")");
229 6 : matlab = getShortcutLabel() + "_pow";
230 : }
231 : // Calcualte coordination numbers as matrix vector times vector of ones
232 90 : readInputLine( getShortcutLabel() + ": MATRIX_VECTOR_PRODUCT ARG=" + matlab + "," + getShortcutLabel() + "_ones");
233 : std::vector<std::string> moments;
234 45 : parseVector("MOMENTS",moments);
235 45 : Tools::interpretRanges( moments );
236 90 : readInputLine( getShortcutLabel() + "_caverage: MEAN ARG=" + getShortcutLabel() + " PERIODIC=NO");
237 47 : for(unsigned i=0; i<moments.size(); ++i) {
238 4 : readInputLine( getShortcutLabel() + "_diffpow-" + moments[i] + ": CUSTOM ARG=" + getShortcutLabel() + "," + getShortcutLabel() + "_caverage PERIODIC=NO FUNC=(x-y)^" + moments[i] );
239 4 : readInputLine( getShortcutLabel() + "_moment-" + moments[i] + ": MEAN ARG=" + getShortcutLabel() + "_diffpow-" + moments[i] + " PERIODIC=NO");
240 : }
241 : // Read in all the shortcut stuff
242 : std::map<std::string,std::string> keymap;
243 45 : multicolvar::MultiColvarShortcuts::readShortcutKeywords( keymap, this );
244 90 : multicolvar::MultiColvarShortcuts::expandFunctions( getShortcutLabel(), getShortcutLabel(), "", keymap, this );
245 90 : }
246 :
247 :
248 : }
249 : }
|
