Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2011-2018 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 "CoordinationBase.h"
23 : #include "tools/SwitchingFunction.h"
24 : #include "ActionRegister.h"
25 :
26 : #include <string>
27 :
28 : using namespace std;
29 :
30 : namespace PLMD {
31 : namespace colvar {
32 :
33 : //+PLUMEDOC COLVAR COORDINATION
34 : /*
35 : Calculate coordination numbers.
36 :
37 : This keyword can be used to calculate the number of contacts between two groups of atoms
38 : and is defined as
39 : \f[
40 : \sum_{i\in A} \sum_{i\in B} s_{ij}
41 : \f]
42 : where \f$s_{ij}\f$ is 1 if the contact between atoms \f$i\f$ and \f$j\f$ is formed,
43 : zero otherwise.
44 : In practise, \f$s_{ij}\f$ is replaced with a switching function to make it differentiable.
45 : The default switching function is:
46 : \f[
47 : s_{ij} = \frac{ 1 - \left(\frac{{\bf r}_{ij}-d_0}{r_0}\right)^n } { 1 - \left(\frac{{\bf r}_{ij}-d_0}{r_0}\right)^m }
48 : \f]
49 : but it can be changed using the optional SWITCH option.
50 :
51 : To make your calculation faster you can use a neighbor list, which makes it that only a
52 : relevant subset of the pairwise distance are calculated at every step.
53 :
54 : If GROUPB is empty, it will sum the \f$\frac{N(N-1)}{2}\f$ pairs in GROUPA. This avoids computing
55 : twice permuted indexes (e.g. pair (i,j) and (j,i)) thus running at twice the speed.
56 :
57 : Notice that if there are common atoms between GROUPA and GROUPB the switching function should be
58 : equal to one. These "self contacts" are discarded by plumed (since version 2.1),
59 : so that they actually count as "zero".
60 :
61 :
62 : \par Examples
63 :
64 : The following example instructs plumed to calculate the total coordination number of the atoms in group 1-10 with the atoms in group 20-100. For atoms 1-10 coordination numbers are calculated that count the number of atoms from the second group that are within 0.3 nm of the central atom. A neighbour list is used to make this calculation faster, this neighbour list is updated every 100 steps.
65 : \verbatim
66 : COORDINATION GROUPA=1-10 GROUPB=20-100 R_0=0.3 NLIST NL_CUTOFF=0.5 NL_STRIDE=100
67 : \endverbatim
68 :
69 : The following is a dummy example which should compute the value 0 because the self interaction
70 : of atom 1 is skipped. Notice that in plumed 2.0 "self interactions" were not skipped, and the
71 : same calculation should return 1.
72 : \verbatim
73 : c: COORDINATION GROUPA=1 GROUPB=1 R_0=0.3
74 : PRINT ARG=c STRIDE=10
75 : \endverbatim
76 :
77 : \verbatim
78 : c1: COORDINATION GROUPA=1-10 GROUPB=1-10 R_0=0.3
79 : x: COORDINATION GROUPA=1-10 R_0=0.3
80 : c2: COMBINE ARG=x COEFFICIENTS=2
81 : # the two variables c1 and c2 should be identical, but the calculation of c2 is twice faster
82 : # since it runs on half of the pairs. Notice that to get the same result you
83 : # should double it
84 : PRINT ARG=c1,c2 STRIDE=10
85 : \endverbatim
86 : See also \ref PRINT and \ref COMBINE
87 :
88 :
89 :
90 : */
91 : //+ENDPLUMEDOC
92 :
93 232 : class Coordination : public CoordinationBase {
94 : SwitchingFunction switchingFunction;
95 :
96 : public:
97 : explicit Coordination(const ActionOptions&);
98 : // active methods:
99 : static void registerKeywords( Keywords& keys );
100 : virtual double pairing(double distance,double&dfunc,unsigned i,unsigned j)const;
101 : };
102 :
103 2639 : PLUMED_REGISTER_ACTION(Coordination,"COORDINATION")
104 :
105 117 : void Coordination::registerKeywords( Keywords& keys ) {
106 117 : CoordinationBase::registerKeywords(keys);
107 117 : keys.add("compulsory","NN","6","The n parameter of the switching function ");
108 117 : keys.add("compulsory","MM","0","The m parameter of the switching function; 0 implies 2*NN");
109 117 : keys.add("compulsory","D_0","0.0","The d_0 parameter of the switching function");
110 117 : keys.add("compulsory","R_0","The r_0 parameter of the switching function");
111 : keys.add("optional","SWITCH","This keyword is used if you want to employ an alternative to the continuous swiching function defined above. "
112 : "The following provides information on the \\ref switchingfunction that are available. "
113 117 : "When this keyword is present you no longer need the NN, MM, D_0 and R_0 keywords.");
114 117 : }
115 :
116 116 : Coordination::Coordination(const ActionOptions&ao):
117 : Action(ao),
118 116 : CoordinationBase(ao)
119 : {
120 :
121 232 : string sw,errors;
122 116 : parse("SWITCH",sw);
123 116 : if(sw.length()>0) {
124 90 : switchingFunction.set(sw,errors);
125 90 : if( errors.length()!=0 ) error("problem reading SWITCH keyword : " + errors );
126 : } else {
127 26 : int nn=6;
128 26 : int mm=0;
129 26 : double d0=0.0;
130 26 : double r0=0.0;
131 26 : parse("R_0",r0);
132 26 : if(r0<=0.0) error("R_0 should be explicitly specified and positive");
133 26 : parse("D_0",d0);
134 26 : parse("NN",nn);
135 26 : parse("MM",mm);
136 26 : switchingFunction.set(nn,mm,r0,d0);
137 : }
138 :
139 116 : checkRead();
140 :
141 232 : log<<" contacts are counted with cutoff "<<switchingFunction.description()<<"\n";
142 116 : }
143 :
144 5703806 : double Coordination::pairing(double distance,double&dfunc,unsigned i,unsigned j)const {
145 : (void) i; // avoid warnings
146 : (void) j; // avoid warnings
147 5703806 : return switchingFunction.calculateSqr(distance,dfunc);
148 : }
149 :
150 : }
151 :
152 2523 : }
|
