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