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 352 : keys.needsAction("CONTACT_MATRIX"); keys.needsAction("GROUP");
115 176 : }
116 :
117 108 : void CoordinationNumbers::expandMatrix( const bool& components, const std::string& lab, const std::string& sp_str,
118 : const std::string& spa_str, const std::string& spb_str, ActionShortcut* action ) {
119 108 : if( sp_str.length()==0 && spa_str.length()==0 ) return;
120 :
121 108 : std::string matinp = lab + "_mat: CONTACT_MATRIX";
122 108 : if( sp_str.length()>0 ) {
123 69 : matinp += " GROUP=" + sp_str;
124 138 : action->readInputLine( lab + "_grp: GROUP ATOMS=" + sp_str );
125 39 : } else if( spa_str.length()>0 ) {
126 78 : matinp += " GROUPA=" + spa_str + " GROUPB=" + spb_str;
127 78 : action->readInputLine( lab + "_grp: GROUP ATOMS=" + spa_str );
128 : }
129 :
130 216 : std::string sw_str; action->parse("SWITCH",sw_str);
131 108 : if( sw_str.length()>0 ) {
132 162 : matinp += " SWITCH={" + sw_str + "}";
133 : } else {
134 81 : std::string r0; action->parse("R_0",r0); std::string d0; action->parse("D_0",d0);
135 27 : if( r0.length()==0 ) action->error("missing switching function parameters use SWITCH/R_0");
136 54 : std::string nn; action->parse("NN",nn); std::string mm; action->parse("MM",mm);
137 54 : matinp += " R_0=" + r0 + " D_0=" + d0 + " NN=" + nn + " MM=" + mm;
138 : }
139 108 : if( components ) matinp += " COMPONENTS";
140 108 : action->readInputLine( matinp );
141 : }
142 :
143 67 : void CoordinationNumbers::registerKeywords( Keywords& keys ) {
144 67 : shortcutKeywords( keys );
145 134 : keys.add("compulsory","R_POWER","the power to which you want to raise the distance");
146 134 : keys.addFlag("LOWMEM",false,"this flag does nothing and is present only to ensure back-compatibility");
147 134 : keys.add("optional","MOMENTS","the list of moments that you would like to calculate");
148 134 : keys.addOutputComponent("moment","MOMENTS","the moments of the distribution");
149 201 : keys.needsAction("MATRIX_VECTOR_PRODUCT"); keys.needsAction("ONES"); keys.needsAction("MOMENTS");
150 67 : }
151 :
152 45 : CoordinationNumbers::CoordinationNumbers(const ActionOptions& ao):
153 : Action(ao),
154 45 : ActionShortcut(ao)
155 : {
156 45 : bool lowmem; parseFlag("LOWMEM",lowmem);
157 45 : if( lowmem ) warning("LOWMEM flag is deprecated and is no longer required for this action");
158 : // Setup the contract matrix if that is what is needed
159 : std::string matlab, sp_str, specA, specB;
160 180 : parse("SPECIES",sp_str); parse("SPECIESA",specA); parse("SPECIESB",specB);
161 45 : if( sp_str.length()>0 || specA.length()>0 ) {
162 90 : matlab = getShortcutLabel() + "_mat"; bool comp=false;
163 48 : if( getName()=="COORDINATION_MOMENTS" ) { comp=true; matlab = getShortcutLabel() + "_mat"; }
164 45 : expandMatrix( comp, getShortcutLabel(), sp_str, specA, specB, this );
165 0 : } else error("missing atoms input use SPECIES or SPECIESA/SPECIESB");
166 45 : ActionWithValue* mb=plumed.getActionSet().selectWithLabel<ActionWithValue*>( matlab );
167 45 : if( !mb ) error("could not find action with name " + matlab );
168 45 : Value* arg=mb->copyOutput(0);
169 45 : if( arg->getRank()!=2 || arg->hasDerivatives() ) error("the input to this action should be a matrix or scalar");
170 : // Create vector of ones to multiply input matrix by
171 45 : std::string nones; Tools::convert( arg->getShape()[1], nones );
172 90 : readInputLine( getShortcutLabel() + "_ones: ONES SIZE=" + nones );
173 45 : if( getName()=="COORDINATION_MOMENTS" ) {
174 : // Calculate the lengths of the vectors
175 3 : std::string r_power; parse("R_POWER",r_power);
176 9 : readInputLine( getShortcutLabel() + "_pow: CUSTOM ARG=" + matlab + ".x," + matlab + ".y," + matlab + ".z," + matlab + ".w VAR=x,y,z,w "
177 6 : + "PERIODIC=NO FUNC=w*(sqrt(x*x+y*y+z*z)^" + r_power +")");
178 6 : matlab = getShortcutLabel() + "_pow";
179 : }
180 : // Calcualte coordination numbers as matrix vector times vector of ones
181 90 : readInputLine( getShortcutLabel() + ": MATRIX_VECTOR_PRODUCT ARG=" + matlab + "," + getShortcutLabel() + "_ones");
182 90 : std::vector<std::string> moments; parseVector("MOMENTS",moments); Tools::interpretRanges( moments );
183 90 : readInputLine( getShortcutLabel() + "_caverage: MEAN ARG=" + getShortcutLabel() + " PERIODIC=NO");
184 47 : for(unsigned i=0; i<moments.size(); ++i) {
185 4 : readInputLine( getShortcutLabel() + "_diffpow-" + moments[i] + ": CUSTOM ARG=" + getShortcutLabel() + "," + getShortcutLabel() + "_caverage PERIODIC=NO FUNC=(x-y)^" + moments[i] );
186 4 : readInputLine( getShortcutLabel() + "_moment-" + moments[i] + ": MEAN ARG=" + getShortcutLabel() + "_diffpow-" + moments[i] + " PERIODIC=NO");
187 : }
188 : // Read in all the shortcut stuff
189 45 : std::map<std::string,std::string> keymap; multicolvar::MultiColvarShortcuts::readShortcutKeywords( keymap, this );
190 90 : multicolvar::MultiColvarShortcuts::expandFunctions( getShortcutLabel(), getShortcutLabel(), "", keymap, this );
191 90 : }
192 :
193 :
194 : }
195 : }
|
