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 "core/ActionRegister.h"
25 : #include "vesselbase/LessThan.h"
26 : #include "vesselbase/Between.h"
27 :
28 : #include <string>
29 : #include <cmath>
30 :
31 : using namespace std;
32 :
33 : namespace PLMD {
34 : namespace multicolvar {
35 :
36 : //+PLUMEDOC MCOLVAR DISTANCES
37 : /*
38 : Calculate the distances between one or many pairs of atoms. You can then calculate functions of the distribution of
39 : distances such as the minimum, the number less than a certain quantity and so on.
40 :
41 : \par Examples
42 :
43 : The following input tells plumed to calculate the distances between atoms 3 and 5 and between atoms 1 and 2 and to
44 : print the minimum for these two distances.
45 : \plumedfile
46 : DISTANCES ATOMS1=3,5 ATOMS2=1,2 MIN={BETA=0.1} LABEL=d1
47 : PRINT ARG=d1.min
48 : \endplumedfile
49 : (See also \ref PRINT).
50 :
51 : The following input tells plumed to calculate the distances between atoms 3 and 5 and between atoms 1 and 2
52 : and then to calculate the number of these distances that are less than 0.1 nm. The number of distances
53 : less than 0.1nm is then printed to a file.
54 : \plumedfile
55 : DISTANCES ATOMS1=3,5 ATOMS2=1,2 LABEL=d1 LESS_THAN={RATIONAL R_0=0.1}
56 : PRINT ARG=d1.lt0.1
57 : \endplumedfile
58 : (See also \ref PRINT \ref switchingfunction).
59 :
60 : The following input tells plumed to calculate all the distances between atoms 1, 2 and 3 (i.e. the distances between atoms
61 : 1 and 2, atoms 1 and 3 and atoms 2 and 3). The average of these distances is then calculated.
62 : \plumedfile
63 : DISTANCES GROUP=1-3 MEAN LABEL=d1
64 : PRINT ARG=d1.mean
65 : \endplumedfile
66 : (See also \ref PRINT)
67 :
68 : The following input tells plumed to calculate all the distances between the atoms in GROUPA and the atoms in GROUPB.
69 : In other words the distances between atoms 1 and 2 and the distance between atoms 1 and 3. The number of distances
70 : more than 0.1 is then printed to a file.
71 : \plumedfile
72 : DISTANCES GROUPA=1 GROUPB=2,3 MORE_THAN={RATIONAL R_0=0.1}
73 : PRINT ARG=d1.gt0.1
74 : \endplumedfile
75 : (See also \ref PRINT \ref switchingfunction)
76 :
77 :
78 : \par Calculating minimum distances
79 :
80 : To calculate and print the minimum distance between two groups of atoms you use the following commands
81 :
82 : \plumedfile
83 : d1: DISTANCES GROUPA=1-10 GROUPB=11-20 MIN={BETA=500.}
84 : PRINT ARG=d1.min FILE=colvar STRIDE=10
85 : \endplumedfile
86 : (see \ref DISTANCES and \ref PRINT)
87 :
88 : In order to ensure differentiability the minimum is calculated using the following function:
89 :
90 : \f[
91 : s = \frac{\beta}{ \log \sum_i \exp\left( \frac{\beta}{s_i} \right) }
92 : \f]
93 :
94 : where \f$\beta\f$ is a user specified parameter.
95 :
96 : This input is used rather than a separate MINDIST colvar so that the same routine and the same input style can be
97 : used to calculate minimum coordinatetion numbers (see \ref COORDINATIONNUMBER), minimum
98 : angles (see \ref ANGLES) and many other variables.
99 :
100 : This new way of calculating mindist is part of plumed 2's multicolvar functionality. These special actions
101 : allow you to calculate multiple functions of a distribution of simple collective variables. As an example you
102 : can calculate the number of distances less than 1.0, the minimum distance, the number of distances more than
103 : 2.0 and the number of distances between 1.0 and 2.0 by using the following command:
104 :
105 : \plumedfile
106 : DISTANCES ...
107 : GROUPA=1-10 GROUPB=11-20
108 : LESS_THAN={RATIONAL R_0=1.0}
109 : MORE_THAN={RATIONAL R_0=2.0}
110 : BETWEEN={GAUSSIAN LOWER=1.0 UPPER=2.0}
111 : MIN={BETA=500.}
112 : ... DISTANCES
113 : PRINT ARG=d1.lessthan,d1.morethan,d1.between,d1.min FILE=colvar STRIDE=10
114 : \endplumedfile
115 : (see \ref DISTANCES and \ref PRINT)
116 :
117 : A calculation performed this way is fast because the expensive part of the calculation - the calculation of all the distances - is only
118 : done once per step. Furthermore, it can be made faster by using the TOL keyword to discard those distance that make only a small contributions
119 : to the final values together with the NL_STRIDE keyword, which ensures that the distances that make only a small contribution to the final values aren't
120 : calculated at every step.
121 :
122 : */
123 : //+ENDPLUMEDOC
124 :
125 :
126 102 : class Distances : public MultiColvarBase {
127 : private:
128 : public:
129 : static void registerKeywords( Keywords& keys );
130 : explicit Distances(const ActionOptions&);
131 : // active methods:
132 : virtual double compute( const unsigned& tindex, AtomValuePack& myatoms ) const ;
133 : /// Returns the number of coordinates of the field
134 111 : bool isPeriodic() { return false; }
135 : };
136 :
137 6503 : PLUMED_REGISTER_ACTION(Distances,"DISTANCES")
138 :
139 52 : void Distances::registerKeywords( Keywords& keys ) {
140 52 : MultiColvarBase::registerKeywords( keys );
141 208 : keys.use("ALT_MIN"); keys.use("LOWEST"); keys.use("HIGHEST");
142 260 : keys.use("MEAN"); keys.use("MIN"); keys.use("MAX"); keys.use("LESS_THAN"); // keys.use("DHENERGY");
143 260 : keys.use("MORE_THAN"); keys.use("BETWEEN"); keys.use("HISTOGRAM"); keys.use("MOMENTS");
144 208 : keys.add("numbered","ATOMS","the atoms involved in each of the distances you wish to calculate. "
145 : "Keywords like ATOMS1, ATOMS2, ATOMS3,... should be listed and one distance will be "
146 : "calculated for each ATOM keyword you specify (all ATOM keywords should "
147 : "specify the indices of two atoms). The eventual number of quantities calculated by this "
148 : "action will depend on what functions of the distribution you choose to calculate.");
149 156 : keys.reset_style("ATOMS","atoms");
150 208 : keys.add("atoms-1","GROUP","Calculate the distance between each distinct pair of atoms in the group");
151 208 : keys.add("atoms-2","GROUPA","Calculate the distances between all the atoms in GROUPA and all "
152 : "the atoms in GROUPB. This must be used in conjuction with GROUPB.");
153 208 : keys.add("atoms-2","GROUPB","Calculate the distances between all the atoms in GROUPA and all the atoms "
154 : "in GROUPB. This must be used in conjuction with GROUPA.");
155 52 : }
156 :
157 51 : Distances::Distances(const ActionOptions&ao):
158 : Action(ao),
159 51 : MultiColvarBase(ao)
160 : {
161 : // Read in the atoms
162 : std::vector<AtomNumber> all_atoms;
163 204 : readTwoGroups( "GROUP", "GROUPA", "GROUPB", all_atoms );
164 86 : if( atom_lab.size()==0 ) readAtomsLikeKeyword( "ATOMS", 2, all_atoms );
165 51 : setupMultiColvarBase( all_atoms );
166 : // And check everything has been read in correctly
167 51 : checkRead();
168 :
169 : // Now check if we can use link cells
170 51 : if( getNumberOfVessels()>0 ) {
171 : bool use_link=false; double rcut;
172 49 : vesselbase::LessThan* lt=dynamic_cast<vesselbase::LessThan*>( getPntrToVessel(0) );
173 49 : if( lt ) {
174 16 : use_link=true; rcut=lt->getCutoff();
175 : } else {
176 33 : vesselbase::Between* bt=dynamic_cast<vesselbase::Between*>( getPntrToVessel(0) );
177 33 : if( bt ) { use_link=true; rcut=bt->getCutoff(); }
178 : }
179 : if( use_link ) {
180 106 : for(unsigned i=1; i<getNumberOfVessels(); ++i) {
181 30 : vesselbase::LessThan* lt2=dynamic_cast<vesselbase::LessThan*>( getPntrToVessel(i) );
182 30 : vesselbase::Between* bt=dynamic_cast<vesselbase::Between*>( getPntrToVessel(i) );
183 30 : if( lt2 ) {
184 0 : double tcut=lt2->getCutoff();
185 0 : if( tcut>rcut ) rcut=tcut;
186 30 : } else if( bt ) {
187 30 : double tcut=bt->getCutoff();
188 30 : if( tcut>rcut ) rcut=tcut;
189 : } else {
190 : use_link=false;
191 : }
192 : }
193 : }
194 49 : if( use_link ) setLinkCellCutoff( rcut );
195 : }
196 51 : }
197 :
198 15275 : double Distances::compute( const unsigned& tindex, AtomValuePack& myatoms ) const {
199 15275 : Vector distance;
200 30550 : distance=getSeparation( myatoms.getPosition(0), myatoms.getPosition(1) );
201 15275 : const double value=distance.modulo();
202 15275 : const double invvalue=1.0/value;
203 :
204 : // And finish the calculation
205 15275 : addAtomDerivatives( 1, 0,-invvalue*distance, myatoms );
206 15275 : addAtomDerivatives( 1, 1, invvalue*distance, myatoms );
207 15275 : myatoms.addBoxDerivatives( 1, -invvalue*Tensor(distance,distance) );
208 15275 : return value;
209 : }
210 :
211 : }
212 4839 : }
213 :
|
