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 { return false; }
93 : };
94 :
95 10523 : PLUMED_REGISTER_ACTION(CoordinationNumbers,"COORDINATIONNUMBER")
96 :
97 53 : void CoordinationNumbers::registerKeywords( Keywords& keys ) {
98 53 : MultiColvarBase::registerKeywords( keys );
99 159 : keys.use("SPECIES"); keys.use("SPECIESA"); keys.use("SPECIESB");
100 106 : keys.add("compulsory","NN","6","The n parameter of the switching function ");
101 106 : keys.add("compulsory","MM","0","The m parameter of the switching function; 0 implies 2*NN");
102 106 : keys.add("compulsory","D_0","0.0","The d_0 parameter of the switching function");
103 106 : keys.add("compulsory","R_0","The r_0 parameter of the switching function");
104 106 : keys.add("optional","R_POWER","Multiply the coordination number function by a power of r, "
105 : "as done in White and Voth (see note above, default: no)");
106 106 : keys.add("optional","SWITCH","This keyword is used if you want to employ an alternative to the continuous switching function defined above. "
107 : "The following provides information on the \\ref switchingfunction that are available. "
108 : "When this keyword is present you no longer need the NN, MM, D_0 and R_0 keywords.");
109 : // Use actionWithDistributionKeywords
110 212 : keys.use("MEAN"); keys.use("MORE_THAN"); keys.use("LESS_THAN"); keys.use("MAX");
111 212 : keys.use("MIN"); keys.use("BETWEEN"); keys.use("HISTOGRAM"); keys.use("MOMENTS");
112 159 : keys.use("ALT_MIN"); keys.use("LOWEST"); keys.use("HIGHEST");
113 53 : }
114 :
115 52 : CoordinationNumbers::CoordinationNumbers(const ActionOptions&ao):
116 : Action(ao),
117 : MultiColvarBase(ao),
118 52 : r_power(0)
119 : {
120 :
121 : // Read in the switching function
122 104 : std::string sw, errors; parse("SWITCH",sw);
123 52 : if(sw.length()>0) {
124 50 : switchingFunction.set(sw,errors);
125 50 : if( errors.length()!=0 ) error("problem reading SWITCH keyword : " + errors );
126 : } else {
127 2 : double r_0=-1.0, d_0; int nn, mm;
128 4 : parse("NN",nn); parse("MM",mm);
129 4 : parse("R_0",r_0); parse("D_0",d_0);
130 2 : if( r_0<0.0 ) error("you must set a value for R_0");
131 2 : switchingFunction.set(nn,mm,r_0,d_0);
132 :
133 : }
134 52 : log.printf(" coordination of central atom and those within %s\n",( switchingFunction.description() ).c_str() );
135 :
136 : //get cutoff of switching function
137 52 : double rcut = switchingFunction.get_dmax();
138 :
139 : //parse power
140 52 : parse("R_POWER", r_power);
141 52 : if(r_power > 0) {
142 4 : log.printf(" Multiplying switching function by r^%d\n", r_power);
143 4 : double offset = switchingFunction.calculate(rcut*0.9999, rcut2) * std::pow(rcut*0.9999, r_power);
144 4 : log.printf(" You will have a discontinuous jump of %f to 0 near the cutoff of your switching function. "
145 : "Consider setting D_MAX or reducing R_POWER if this is large\n", offset);
146 : }
147 :
148 : // Set the link cell cutoff
149 52 : setLinkCellCutoff( rcut );
150 52 : rcut2 = rcut * rcut;
151 :
152 : // And setup the ActionWithVessel
153 52 : std::vector<AtomNumber> all_atoms; setupMultiColvarBase( all_atoms ); checkRead();
154 52 : }
155 :
156 37088 : double CoordinationNumbers::compute( const unsigned& tindex, AtomValuePack& myatoms ) const {
157 : // Calculate the coordination number
158 : double dfunc, sw, d, raised;
159 3822803 : for(unsigned i=1; i<myatoms.getNumberOfAtoms(); ++i) {
160 : Vector& distance=myatoms.getPosition(i);
161 : double d2;
162 7107372 : if ( (d2=distance[0]*distance[0])<rcut2 &&
163 3321657 : (d2+=distance[1]*distance[1])<rcut2 &&
164 6676366 : (d2+=distance[2]*distance[2])<rcut2 &&
165 : d2>epsilon ) {
166 :
167 2633753 : sw = switchingFunction.calculateSqr( d2, dfunc );
168 2633753 : if(r_power > 0) {
169 19350 : d = std::sqrt(d2); raised = std::pow( d, r_power - 1 );
170 19350 : accumulateSymmetryFunction( 1, i, sw * raised * d,
171 19350 : (dfunc * d * raised + sw * r_power * raised / d) * distance,
172 38700 : (-dfunc * d * raised - sw * r_power * raised / d) * Tensor(distance, distance),
173 : myatoms );
174 : } else {
175 2614403 : accumulateSymmetryFunction( 1, i, sw, (dfunc)*distance, (-dfunc)*Tensor(distance,distance), myatoms );
176 : }
177 : }
178 : }
179 :
180 37088 : return myatoms.getValue(1);
181 : }
182 :
183 : }
184 : }
|
