Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2011-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 "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 : \plumedfile
66 : COORDINATION GROUPA=1-10 GROUPB=20-100 R_0=0.3 NLIST NL_CUTOFF=0.5 NL_STRIDE=100
67 : \endplumedfile
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 : \plumedfile
73 : c: COORDINATION GROUPA=1 GROUPB=1 R_0=0.3
74 : PRINT ARG=c STRIDE=10
75 : \endplumedfile
76 :
77 : Here's an example that shows what happens when providing COORDINATION with
78 : a single group:
79 : \plumedfile
80 : # define some huge group:
81 : group: GROUP ATOMS=1-1000
82 : # Here's coordination of a group against itself:
83 : c1: COORDINATION GROUPA=group GROUPB=group R_0=0.3
84 : # Here's coordination within a single group:
85 : x: COORDINATION GROUPA=group R_0=0.3
86 : # This is just multiplying times 2 the variable x:
87 : c2: COMBINE ARG=x COEFFICIENTS=2
88 :
89 : # the two variables c1 and c2 should be identical, but the calculation of c2 is twice faster
90 : # since it runs on half of the pairs.
91 : PRINT ARG=c1,c2 STRIDE=10
92 : \endplumedfile
93 :
94 :
95 :
96 : */
97 : //+ENDPLUMEDOC
98 :
99 300 : class Coordination : public CoordinationBase {
100 : SwitchingFunction switchingFunction;
101 :
102 : public:
103 : explicit Coordination(const ActionOptions&);
104 : // active methods:
105 : static void registerKeywords( Keywords& keys );
106 : virtual double pairing(double distance,double&dfunc,unsigned i,unsigned j)const;
107 : };
108 :
109 6605 : PLUMED_REGISTER_ACTION(Coordination,"COORDINATION")
110 :
111 154 : void Coordination::registerKeywords( Keywords& keys ) {
112 154 : CoordinationBase::registerKeywords(keys);
113 770 : keys.add("compulsory","NN","6","The n parameter of the switching function ");
114 770 : keys.add("compulsory","MM","0","The m parameter of the switching function; 0 implies 2*NN");
115 770 : keys.add("compulsory","D_0","0.0","The d_0 parameter of the switching function");
116 616 : keys.add("compulsory","R_0","The r_0 parameter of the switching function");
117 616 : keys.add("optional","SWITCH","This keyword is used if you want to employ an alternative to the continuous swiching function defined above. "
118 : "The following provides information on the \\ref switchingfunction that are available. "
119 : "When this keyword is present you no longer need the NN, MM, D_0 and R_0 keywords.");
120 154 : }
121 :
122 153 : Coordination::Coordination(const ActionOptions&ao):
123 : Action(ao),
124 156 : CoordinationBase(ao)
125 : {
126 :
127 : string sw,errors;
128 306 : parse("SWITCH",sw);
129 153 : if(sw.length()>0) {
130 111 : switchingFunction.set(sw,errors);
131 114 : if( errors.length()!=0 ) error("problem reading SWITCH keyword : " + errors );
132 : } else {
133 42 : int nn=6;
134 42 : int mm=0;
135 42 : double d0=0.0;
136 42 : double r0=0.0;
137 84 : parse("R_0",r0);
138 42 : if(r0<=0.0) error("R_0 should be explicitly specified and positive");
139 84 : parse("D_0",d0);
140 84 : parse("NN",nn);
141 80 : parse("MM",mm);
142 40 : switchingFunction.set(nn,mm,r0,d0);
143 : }
144 :
145 150 : checkRead();
146 :
147 303 : log<<" contacts are counted with cutoff "<<switchingFunction.description()<<"\n";
148 150 : }
149 :
150 6074544 : double Coordination::pairing(double distance,double&dfunc,unsigned i,unsigned j)const {
151 : (void) i; // avoid warnings
152 : (void) j; // avoid warnings
153 6074544 : return switchingFunction.calculateSqr(distance,dfunc);
154 : }
155 :
156 : }
157 :
158 4839 : }
|
