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