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 "MultiColvarBase.h"
23 : #include "AtomValuePack.h"
24 : #include "tools/NeighborList.h"
25 : #include "core/ActionRegister.h"
26 : #include "tools/SwitchingFunction.h"
27 :
28 : #include <string>
29 : #include <cmath>
30 :
31 : namespace PLMD {
32 : namespace multicolvar {
33 :
34 : //+PLUMEDOC MCOLVAR COORDINATIONNUMBER
35 : /*
36 : Calculate the coordination numbers of atoms so that you can then calculate functions of the distribution of
37 : coordination numbers such as the minimum, the number less than a certain quantity and so on.
38 :
39 : So that the calculated coordination numbers have continuous derivatives the following function is used:
40 :
41 : \f[
42 : s = \frac{ 1 - \left(\frac{r-d_0}{r_0}\right)^n } { 1 - \left(\frac{r-d_0}{r_0}\right)^m }
43 : \f]
44 :
45 : If R_POWER is set, this will use the product of pairwise distance
46 : raised to the R_POWER with the coordination number function defined
47 : above. This was used in White and Voth \cite white2014efficient as a
48 : way of indirectly biasing radial distribution functions. Note that in
49 : that reference this function is referred to as moments of coordination
50 : number, but here we call them powers to distinguish from the existing
51 : MOMENTS keyword of Multicolvars.
52 :
53 : \par Examples
54 :
55 : The following input tells plumed to calculate the coordination numbers of atoms 1-100 with themselves.
56 : The minimum coordination number is then calculated.
57 : \plumedfile
58 : COORDINATIONNUMBER SPECIES=1-100 R_0=1.0 MIN={BETA=0.1}
59 : \endplumedfile
60 :
61 : The following input tells plumed to calculate how many atoms from 1-100 are within 3.0 of each of the atoms
62 : from 101-110. In the first 101 is the central atom, in the second 102 is the central atom and so on. The
63 : number of coordination numbers more than 6 is then computed.
64 : \plumedfile
65 : 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}
66 : \endplumedfile
67 :
68 : The following input tells plumed to calculate the mean coordination number of all atoms with themselves
69 : and its powers. An explicit cutoff is set for each of 8.
70 : \plumedfile
71 : cn0: COORDINATIONNUMBER SPECIES=1-10 SWITCH={RATIONAL R_0=1.0 D_MAX=8} MEAN
72 : cn1: COORDINATIONNUMBER SPECIES=1-10 SWITCH={RATIONAL R_0=1.0 D_MAX=8} R_POWER=1 MEAN
73 : cn2: COORDINATIONNUMBER SPECIES=1-10 SWITCH={RATIONAL R_0=1.0 D_MAX=8} R_POWER=2 MEAN
74 : PRINT ARG=cn0.mean,cn1.mean,cn2.mean STRIDE=1 FILE=cn_out
75 : \endplumedfile
76 :
77 : */
78 : //+ENDPLUMEDOC
79 :
80 :
81 : class CoordinationNumbers : public MultiColvarBase {
82 : private:
83 : double rcut2;
84 : int r_power;
85 : SwitchingFunction switchingFunction;
86 : public:
87 : static void registerKeywords( Keywords& keys );
88 : explicit CoordinationNumbers(const ActionOptions&);
89 : // active methods:
90 : double compute( const unsigned& tindex, AtomValuePack& myatoms ) const override;
91 : /// Returns the number of coordinates of the field
92 856 : bool isPeriodic() override {
93 856 : return false;
94 : }
95 : };
96 :
97 13889 : PLUMED_REGISTER_ACTION(CoordinationNumbers,"COORDINATIONNUMBER")
98 :
99 56 : void CoordinationNumbers::registerKeywords( Keywords& keys ) {
100 56 : MultiColvarBase::registerKeywords( keys );
101 56 : keys.use("SPECIES");
102 56 : keys.use("SPECIESA");
103 56 : keys.use("SPECIESB");
104 112 : keys.add("compulsory","NN","6","The n parameter of the switching function ");
105 112 : keys.add("compulsory","MM","0","The m parameter of the switching function; 0 implies 2*NN");
106 112 : keys.add("compulsory","D_0","0.0","The d_0 parameter of the switching function");
107 112 : keys.add("compulsory","R_0","The r_0 parameter of the switching function");
108 112 : keys.add("optional","R_POWER","Multiply the coordination number function by a power of r, "
109 : "as done in White and Voth (see note above, default: no)");
110 112 : keys.add("optional","SWITCH","This keyword is used if you want to employ an alternative to the continuous switching function defined above. "
111 : "The following provides information on the \\ref switchingfunction that are available. "
112 : "When this keyword is present you no longer need the NN, MM, D_0 and R_0 keywords.");
113 : // Use actionWithDistributionKeywords
114 56 : keys.use("MEAN");
115 56 : keys.use("MORE_THAN");
116 56 : keys.use("LESS_THAN");
117 56 : keys.use("MAX");
118 56 : keys.use("MIN");
119 56 : keys.use("BETWEEN");
120 56 : keys.use("HISTOGRAM");
121 56 : keys.use("MOMENTS");
122 56 : keys.use("ALT_MIN");
123 56 : keys.use("LOWEST");
124 56 : keys.use("HIGHEST");
125 56 : }
126 :
127 52 : CoordinationNumbers::CoordinationNumbers(const ActionOptions&ao):
128 : Action(ao),
129 : MultiColvarBase(ao),
130 52 : r_power(0) {
131 :
132 : // Read in the switching function
133 : std::string sw, errors;
134 104 : parse("SWITCH",sw);
135 52 : if(sw.length()>0) {
136 50 : switchingFunction.set(sw,errors);
137 50 : if( errors.length()!=0 ) {
138 0 : error("problem reading SWITCH keyword : " + errors );
139 : }
140 : } else {
141 2 : double r_0=-1.0, d_0;
142 : int nn, mm;
143 2 : parse("NN",nn);
144 2 : parse("MM",mm);
145 2 : parse("R_0",r_0);
146 2 : parse("D_0",d_0);
147 2 : if( r_0<0.0 ) {
148 0 : error("you must set a value for R_0");
149 : }
150 2 : switchingFunction.set(nn,mm,r_0,d_0);
151 :
152 : }
153 52 : log.printf(" coordination of central atom and those within %s\n",( switchingFunction.description() ).c_str() );
154 :
155 : //get cutoff of switching function
156 52 : double rcut = switchingFunction.get_dmax();
157 :
158 : //parse power
159 52 : parse("R_POWER", r_power);
160 52 : if(r_power > 0) {
161 4 : log.printf(" Multiplying switching function by r^%d\n", r_power);
162 4 : double offset = switchingFunction.calculate(rcut*0.9999, rcut2) * std::pow(rcut*0.9999, r_power);
163 4 : log.printf(" You will have a discontinuous jump of %f to 0 near the cutoff of your switching function. "
164 : "Consider setting D_MAX or reducing R_POWER if this is large\n", offset);
165 : }
166 :
167 : // Set the link cell cutoff
168 52 : setLinkCellCutoff( rcut );
169 52 : rcut2 = rcut * rcut;
170 :
171 : // And setup the ActionWithVessel
172 : std::vector<AtomNumber> all_atoms;
173 52 : setupMultiColvarBase( all_atoms );
174 52 : checkRead();
175 52 : }
176 :
177 37088 : double CoordinationNumbers::compute( const unsigned& tindex, AtomValuePack& myatoms ) const {
178 : // Calculate the coordination number
179 : double dfunc, sw, d, raised;
180 3822803 : for(unsigned i=1; i<myatoms.getNumberOfAtoms(); ++i) {
181 : Vector& distance=myatoms.getPosition(i);
182 : double d2;
183 7107372 : if ( (d2=distance[0]*distance[0])<rcut2 &&
184 3321657 : (d2+=distance[1]*distance[1])<rcut2 &&
185 6676366 : (d2+=distance[2]*distance[2])<rcut2 &&
186 : d2>epsilon ) {
187 :
188 2633753 : sw = switchingFunction.calculateSqr( d2, dfunc );
189 2633753 : if(r_power > 0) {
190 19350 : d = std::sqrt(d2);
191 19350 : raised = std::pow( d, r_power - 1 );
192 19350 : accumulateSymmetryFunction( 1, i, sw * raised * d,
193 19350 : (dfunc * d * raised + sw * r_power * raised / d) * distance,
194 38700 : (-dfunc * d * raised - sw * r_power * raised / d) * Tensor(distance, distance),
195 : myatoms );
196 : } else {
197 2614403 : accumulateSymmetryFunction( 1, i, sw, (dfunc)*distance, (-dfunc)*Tensor(distance,distance), myatoms );
198 : }
199 : }
200 : }
201 :
202 37088 : return myatoms.getValue(1);
203 : }
204 :
205 : }
206 : }
|
