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 : So that the calculated coordination numbers have continuous derivatives the following function is used:
41 :
42 : \f[
43 : s = \frac{ 1 - \left(\frac{r-d_0}{r_0}\right)^n } { 1 - \left(\frac{r-d_0}{r_0}\right)^m }
44 : \f]
45 :
46 : If R_POWER is set, this will use the product of pairwise distance
47 : raised to the R_POWER with the coordination number function defined
48 : above. This was used in White and Voth \cite white2014efficient as a
49 : way of indirectly biasing radial distribution functions. Note that in
50 : that reference this function is referred to as moments of coordination
51 : number, but here we call them powers to distinguish from the existing
52 : MOMENTS keyword of Multicolvars.
53 :
54 : \par Examples
55 :
56 : The following input tells plumed to calculate the coordination numbers of atoms 1-100 with themselves.
57 : The minimum coordination number is then calculated.
58 : \plumedfile
59 : COORDINATIONNUMBER SPECIES=1-100 R_0=1.0 MIN={BETA=0.1}
60 : \endplumedfile
61 :
62 : The following input tells plumed to calculate how many atoms from 1-100 are within 3.0 of each of the atoms
63 : from 101-110. In the first 101 is the central atom, in the second 102 is the central atom and so on. The
64 : number of coordination numbers more than 6 is then computed.
65 : \plumedfile
66 : 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}
67 : \endplumedfile
68 :
69 : The following input tells plumed to calculate the mean coordination number of all atoms with themselves
70 : and its powers. An explicit cutoff is set for each of 8.
71 : \plumedfile
72 : cn0: COORDINATIONNUMBER SPECIES=1-10 SWITCH={RATIONAL R_0=1.0 D_MAX=8} MEAN
73 : cn1: COORDINATIONNUMBER SPECIES=1-10 SWITCH={RATIONAL R_0=1.0 D_MAX=8} R_POWER=1 MEAN
74 : cn2: COORDINATIONNUMBER SPECIES=1-10 SWITCH={RATIONAL R_0=1.0 D_MAX=8} R_POWER=2 MEAN
75 : PRINT ARG=cn0.mean,cn1.mean,cn2.mean STRIDE=1 FILE=cn_out
76 : \endplumedfile
77 :
78 : */
79 : //+ENDPLUMEDOC
80 :
81 : //+PLUMEDOC MCOLVAR COORDINATION_MOMENTS
82 : /*
83 : Calculate moments of the distribution of distances in the first coordination sphere
84 :
85 : \par Examples
86 :
87 : */
88 : //+ENDPLUMEDOC
89 :
90 :
91 : PLUMED_REGISTER_ACTION(CoordinationNumbers,"COORDINATIONNUMBER")
92 : PLUMED_REGISTER_ACTION(CoordinationNumbers,"COORDINATION_MOMENTS")
93 :
94 176 : void CoordinationNumbers::shortcutKeywords( Keywords& keys ) {
95 176 : ActionShortcut::registerKeywords( keys );
96 352 : keys.add("atoms-3","SPECIES","this keyword is used for colvars such as coordination number. In that context it specifies that plumed should calculate "
97 : "one coordination number for each of the atoms specified. Each of these coordination numbers specifies how many of the "
98 : "other specified atoms are within a certain cutoff of the central atom. You can specify the atoms here as another multicolvar "
99 : "action or using a MultiColvarFilter or ActionVolume action. When you do so the quantity is calculated for those atoms specified "
100 : "in the previous multicolvar. This is useful if you would like to calculate the Steinhardt parameter for those atoms that have a "
101 : "coordination number more than four for example");
102 352 : keys.add("atoms-4","SPECIESA","this keyword is used for colvars such as the coordination number. In that context it species that plumed should calculate "
103 : "one coordination number for each of the atoms specified in SPECIESA. Each of these cooordination numbers specifies how many "
104 : "of the atoms specifies using SPECIESB is within the specified cutoff. As with the species keyword the input can also be specified "
105 : "using the label of another multicolvar");
106 352 : keys.add("atoms-4","SPECIESB","this keyword is used for colvars such as the coordination number. It must appear with SPECIESA. For a full explanation see "
107 : "the documentation for that keyword");
108 352 : keys.add("compulsory","NN","6","The n parameter of the switching function ");
109 352 : keys.add("compulsory","MM","0","The m parameter of the switching function; 0 implies 2*NN");
110 352 : keys.add("compulsory","D_0","0.0","The d_0 parameter of the switching function");
111 352 : keys.add("compulsory","R_0","The r_0 parameter of the switching function");
112 352 : keys.add("optional","SWITCH","the switching function that it used in the construction of the contact matrix");
113 176 : multicolvar::MultiColvarShortcuts::shortcutKeywords( keys );
114 176 : keys.needsAction("CONTACT_MATRIX");
115 176 : keys.needsAction("GROUP");
116 176 : }
117 :
118 108 : void CoordinationNumbers::expandMatrix( const bool& components, const std::string& lab, const std::string& sp_str,
119 : const std::string& spa_str, const std::string& spb_str, ActionShortcut* action ) {
120 108 : if( sp_str.length()==0 && spa_str.length()==0 ) {
121 0 : return;
122 : }
123 :
124 108 : std::string matinp = lab + "_mat: CONTACT_MATRIX";
125 108 : if( sp_str.length()>0 ) {
126 69 : matinp += " GROUP=" + sp_str;
127 138 : action->readInputLine( lab + "_grp: GROUP ATOMS=" + sp_str );
128 39 : } else if( spa_str.length()>0 ) {
129 78 : matinp += " GROUPA=" + spa_str + " GROUPB=" + spb_str;
130 78 : action->readInputLine( lab + "_grp: GROUP ATOMS=" + spa_str );
131 : }
132 :
133 : std::string sw_str;
134 216 : action->parse("SWITCH",sw_str);
135 108 : if( sw_str.length()>0 ) {
136 162 : matinp += " SWITCH={" + sw_str + "}";
137 : } else {
138 : std::string r0;
139 54 : action->parse("R_0",r0);
140 : std::string d0;
141 54 : action->parse("D_0",d0);
142 27 : if( r0.length()==0 ) {
143 0 : action->error("missing switching function parameters use SWITCH/R_0");
144 : }
145 : std::string nn;
146 54 : action->parse("NN",nn);
147 : std::string mm;
148 27 : action->parse("MM",mm);
149 54 : matinp += " R_0=" + r0 + " D_0=" + d0 + " NN=" + nn + " MM=" + mm;
150 : }
151 108 : if( components ) {
152 : matinp += " COMPONENTS";
153 : }
154 108 : action->readInputLine( matinp );
155 : }
156 :
157 67 : void CoordinationNumbers::registerKeywords( Keywords& keys ) {
158 67 : shortcutKeywords( keys );
159 134 : keys.add("compulsory","R_POWER","the power to which you want to raise the distance");
160 134 : keys.addFlag("LOWMEM",false,"this flag does nothing and is present only to ensure back-compatibility");
161 134 : keys.add("optional","MOMENTS","the list of moments that you would like to calculate");
162 134 : keys.addOutputComponent("moment","MOMENTS","the moments of the distribution");
163 67 : keys.needsAction("MATRIX_VECTOR_PRODUCT");
164 67 : keys.needsAction("ONES");
165 67 : keys.needsAction("MOMENTS");
166 67 : }
167 :
168 45 : CoordinationNumbers::CoordinationNumbers(const ActionOptions& ao):
169 : Action(ao),
170 45 : ActionShortcut(ao) {
171 : bool lowmem;
172 45 : parseFlag("LOWMEM",lowmem);
173 45 : if( lowmem ) {
174 0 : warning("LOWMEM flag is deprecated and is no longer required for this action");
175 : }
176 : // Setup the contract matrix if that is what is needed
177 : std::string matlab, sp_str, specA, specB;
178 45 : parse("SPECIES",sp_str);
179 45 : parse("SPECIESA",specA);
180 90 : parse("SPECIESB",specB);
181 45 : if( sp_str.length()>0 || specA.length()>0 ) {
182 45 : matlab = getShortcutLabel() + "_mat";
183 45 : bool comp=false;
184 45 : if( getName()=="COORDINATION_MOMENTS" ) {
185 3 : comp=true;
186 6 : matlab = getShortcutLabel() + "_mat";
187 : }
188 45 : expandMatrix( comp, getShortcutLabel(), sp_str, specA, specB, this );
189 : } else {
190 0 : error("missing atoms input use SPECIES or SPECIESA/SPECIESB");
191 : }
192 45 : ActionWithValue* mb=plumed.getActionSet().selectWithLabel<ActionWithValue*>( matlab );
193 45 : if( !mb ) {
194 0 : error("could not find action with name " + matlab );
195 : }
196 45 : Value* arg=mb->copyOutput(0);
197 45 : if( arg->getRank()!=2 || arg->hasDerivatives() ) {
198 0 : error("the input to this action should be a matrix or scalar");
199 : }
200 : // Create vector of ones to multiply input matrix by
201 : std::string nones;
202 45 : Tools::convert( arg->getShape()[1], nones );
203 90 : readInputLine( getShortcutLabel() + "_ones: ONES SIZE=" + nones );
204 45 : if( getName()=="COORDINATION_MOMENTS" ) {
205 : // Calculate the lengths of the vectors
206 : std::string r_power;
207 3 : parse("R_POWER",r_power);
208 9 : readInputLine( getShortcutLabel() + "_pow: CUSTOM ARG=" + matlab + ".x," + matlab + ".y," + matlab + ".z," + matlab + ".w VAR=x,y,z,w "
209 6 : + "PERIODIC=NO FUNC=w*(sqrt(x*x+y*y+z*z)^" + r_power +")");
210 6 : matlab = getShortcutLabel() + "_pow";
211 : }
212 : // Calcualte coordination numbers as matrix vector times vector of ones
213 90 : readInputLine( getShortcutLabel() + ": MATRIX_VECTOR_PRODUCT ARG=" + matlab + "," + getShortcutLabel() + "_ones");
214 : std::vector<std::string> moments;
215 45 : parseVector("MOMENTS",moments);
216 45 : Tools::interpretRanges( moments );
217 90 : readInputLine( getShortcutLabel() + "_caverage: MEAN ARG=" + getShortcutLabel() + " PERIODIC=NO");
218 47 : for(unsigned i=0; i<moments.size(); ++i) {
219 4 : readInputLine( getShortcutLabel() + "_diffpow-" + moments[i] + ": CUSTOM ARG=" + getShortcutLabel() + "," + getShortcutLabel() + "_caverage PERIODIC=NO FUNC=(x-y)^" + moments[i] );
220 4 : readInputLine( getShortcutLabel() + "_moment-" + moments[i] + ": MEAN ARG=" + getShortcutLabel() + "_diffpow-" + moments[i] + " PERIODIC=NO");
221 : }
222 : // Read in all the shortcut stuff
223 : std::map<std::string,std::string> keymap;
224 45 : multicolvar::MultiColvarShortcuts::readShortcutKeywords( keymap, this );
225 90 : multicolvar::MultiColvarShortcuts::expandFunctions( getShortcutLabel(), getShortcutLabel(), "", keymap, this );
226 90 : }
227 :
228 :
229 : }
230 : }
|
