Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2011-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 "Colvar.h"
23 : #include "ColvarShortcut.h"
24 : #include "MultiColvarTemplate.h"
25 : #include "core/ActionRegister.h"
26 : #include "tools/Pbc.h"
27 :
28 : namespace PLMD {
29 : namespace colvar {
30 :
31 : //+PLUMEDOC COLVAR DISTANCE
32 : /*
33 : Calculate the distance between a pair of atoms.
34 :
35 : By default the distance is computed taking into account periodic
36 : boundary conditions. This behavior can be changed with the NOPBC flag.
37 : Moreover, single components in Cartesian space (x,y, and z, with COMPONENTS)
38 : or single components projected to the three lattice vectors (a,b, and c, with SCALED_COMPONENTS)
39 : can be also computed.
40 :
41 : Notice that Cartesian components will not have the proper periodicity!
42 : If you have to study e.g. the permeation of a molecule across a membrane,
43 : better to use SCALED_COMPONENTS.
44 :
45 : \par Examples
46 :
47 : The following input tells plumed to print the distance between atoms 3 and 5,
48 : the distance between atoms 2 and 4 and the x component of the distance between atoms 2 and 4.
49 : \plumedfile
50 : d1: DISTANCE ATOMS=3,5
51 : d2: DISTANCE ATOMS=2,4
52 : d2c: DISTANCE ATOMS=2,4 COMPONENTS
53 : PRINT ARG=d1,d2,d2c.x
54 : \endplumedfile
55 :
56 : The following input computes the end-to-end distance for a polymer
57 : of 100 atoms and keeps it at a value around 5.
58 : \plumedfile
59 : WHOLEMOLECULES ENTITY0=1-100
60 : e2e: DISTANCE ATOMS=1,100 NOPBC
61 : RESTRAINT ARG=e2e KAPPA=1 AT=5
62 : \endplumedfile
63 :
64 : Notice that NOPBC is used
65 : to be sure that if the end-to-end distance is larger than half the simulation
66 : box the distance is compute properly. Also notice that, since many MD
67 : codes break molecules across cell boundary, it might be necessary to
68 : use the \ref WHOLEMOLECULES keyword (also notice that it should be
69 : _before_ distance). The list of atoms provided to \ref WHOLEMOLECULES
70 : here contains all the atoms between 1 and 100. Strictly speaking, this
71 : is not necessary. If you know for sure that atoms with difference in
72 : the index say equal to 10 are _not_ going to be farther than half cell
73 : you can e.g. use
74 : \plumedfile
75 : WHOLEMOLECULES ENTITY0=1,10,20,30,40,50,60,70,80,90,100
76 : e2e: DISTANCE ATOMS=1,100 NOPBC
77 : RESTRAINT ARG=e2e KAPPA=1 AT=5
78 : \endplumedfile
79 : Just be sure that the ordered list provide to \ref WHOLEMOLECULES has the following
80 : properties:
81 : - Consecutive atoms should be closer than half-cell throughout the entire simulation.
82 : - Atoms required later for the distance (e.g. 1 and 100) should be included in the list
83 :
84 : The following example shows how to take into account periodicity e.g.
85 : in z-component of a distance
86 : \plumedfile
87 : # this is a center of mass of a large group
88 : c: COM ATOMS=1-100
89 : # this is the distance between atom 101 and the group
90 : d: DISTANCE ATOMS=c,101 COMPONENTS
91 : # this makes a new variable, dd, equal to d and periodic, with domain -10,10
92 : # this is the right choise if e.g. the cell is orthorombic and its size in
93 : # z direction is 20.
94 : dz: COMBINE ARG=d.z PERIODIC=-10,10
95 : # metadynamics on dd
96 : METAD ARG=dz SIGMA=0.1 HEIGHT=0.1 PACE=200
97 : \endplumedfile
98 :
99 : Using SCALED_COMPONENTS this problem should not arise because they are always periodic
100 : with domain (-0.5,+0.5).
101 :
102 :
103 :
104 :
105 : */
106 : //+ENDPLUMEDOC
107 :
108 : //+PLUMEDOC MCOLVAR DISTANCE_SCALAR
109 : /*
110 : Calculate the distance between a pair of atoms
111 :
112 : \par Examples
113 :
114 : */
115 : //+ENDPLUMEDOC
116 :
117 : //+PLUMEDOC MCOLVAR DISTANCE_VECTOR
118 : /*
119 : Calculate a vector containing the distances between various pairs of atoms
120 :
121 : \par Examples
122 :
123 : */
124 : //+ENDPLUMEDOC
125 :
126 : class Distance : public Colvar {
127 : bool components;
128 : bool scaled_components;
129 : bool pbc;
130 :
131 : std::vector<double> value, masses, charges;
132 : std::vector<std::vector<Vector> > derivs;
133 : std::vector<Tensor> virial;
134 : public:
135 : static void registerKeywords( Keywords& keys );
136 : explicit Distance(const ActionOptions&);
137 : static void parseAtomList( const int& num, std::vector<AtomNumber>& t, ActionAtomistic* aa );
138 : static unsigned getModeAndSetupValues( ActionWithValue* av );
139 : // active methods:
140 : void calculate() override;
141 : static void calculateCV( const unsigned& mode, const std::vector<double>& masses, const std::vector<double>& charges,
142 : const std::vector<Vector>& pos, std::vector<double>& vals, std::vector<std::vector<Vector> >& derivs,
143 : std::vector<Tensor>& virial, const ActionAtomistic* aa );
144 : };
145 :
146 : typedef ColvarShortcut<Distance> DistanceShortcut;
147 : PLUMED_REGISTER_ACTION(DistanceShortcut,"DISTANCE")
148 : PLUMED_REGISTER_ACTION(Distance,"DISTANCE_SCALAR")
149 : typedef MultiColvarTemplate<Distance> DistanceMulti;
150 : PLUMED_REGISTER_ACTION(DistanceMulti,"DISTANCE_VECTOR")
151 :
152 2380 : void Distance::registerKeywords( Keywords& keys ) {
153 2380 : Colvar::registerKeywords( keys ); keys.setDisplayName("DISTANCE");
154 4760 : keys.add("atoms","ATOMS","the pair of atom that we are calculating the distance between");
155 4760 : keys.addFlag("COMPONENTS",false,"calculate the x, y and z components of the distance separately and store them as label.x, label.y and label.z");
156 4760 : keys.addFlag("SCALED_COMPONENTS",false,"calculate the a, b and c scaled components of the distance separately and store them as label.a, label.b and label.c");
157 4760 : keys.addOutputComponent("x","COMPONENTS","scalar/vector","the x-component of the vector connecting the two atoms");
158 4760 : keys.addOutputComponent("y","COMPONENTS","scalar/vector","the y-component of the vector connecting the two atoms");
159 4760 : keys.addOutputComponent("z","COMPONENTS","scalar/vector","the z-component of the vector connecting the two atoms");
160 4760 : keys.addOutputComponent("a","SCALED_COMPONENTS","scalar/vector","the normalized projection on the first lattice vector of the vector connecting the two atoms");
161 4760 : keys.addOutputComponent("b","SCALED_COMPONENTS","scalar/vector","the normalized projection on the second lattice vector of the vector connecting the two atoms");
162 4760 : keys.addOutputComponent("c","SCALED_COMPONENTS","scalar/vector","the normalized projection on the third lattice vector of the vector connecting the two atoms");
163 4760 : keys.add("hidden","NO_ACTION_LOG","suppresses printing from action on the log");
164 4760 : keys.setValueDescription("scalar/vector","the DISTANCE between this pair of atoms");
165 2380 : }
166 :
167 623 : Distance::Distance(const ActionOptions&ao):
168 : PLUMED_COLVAR_INIT(ao),
169 623 : components(false),
170 623 : scaled_components(false),
171 623 : pbc(true),
172 623 : value(1),
173 625 : derivs(1),
174 1246 : virial(1)
175 : {
176 623 : derivs[0].resize(2);
177 : std::vector<AtomNumber> atoms;
178 623 : parseAtomList(-1,atoms,this);
179 623 : if(atoms.size()!=2)
180 1 : error("Number of specified atoms should be 2");
181 :
182 622 : bool nopbc=!pbc;
183 624 : parseFlag("NOPBC",nopbc);
184 622 : pbc=!nopbc;
185 :
186 622 : if(pbc) log.printf(" using periodic boundary conditions\n");
187 4 : else log.printf(" without periodic boundary conditions\n");
188 :
189 622 : unsigned mode = getModeAndSetupValues( this );
190 621 : if(mode==1) components=true; else if(mode==2) scaled_components=true;
191 621 : if( components || scaled_components ) {
192 40 : value.resize(3); derivs.resize(3); virial.resize(3);
193 160 : for(unsigned i=0; i<3; ++i) derivs[i].resize(2);
194 : }
195 621 : requestAtoms(atoms);
196 627 : }
197 :
198 30941 : void Distance::parseAtomList( const int& num, std::vector<AtomNumber>& t, ActionAtomistic* aa ) {
199 61882 : aa->parseAtomList("ATOMS",num,t);
200 30941 : if( t.size()==2 ) aa->log.printf(" between atoms %d %d\n",t[0].serial(),t[1].serial());
201 30941 : }
202 :
203 730 : unsigned Distance::getModeAndSetupValues( ActionWithValue* av ) {
204 730 : bool c; av->parseFlag("COMPONENTS",c);
205 730 : bool sc; av->parseFlag("SCALED_COMPONENTS",sc);
206 730 : if( c && sc ) av->error("COMPONENTS and SCALED_COMPONENTS are not compatible");
207 :
208 729 : if(c) {
209 160 : av->addComponentWithDerivatives("x"); av->componentIsNotPeriodic("x");
210 160 : av->addComponentWithDerivatives("y"); av->componentIsNotPeriodic("y");
211 160 : av->addComponentWithDerivatives("z"); av->componentIsNotPeriodic("z");
212 80 : av->log<<" WARNING: components will not have the proper periodicity - see manual\n";
213 80 : return 1;
214 649 : } else if(sc) {
215 18 : av->addComponentWithDerivatives("a"); av->componentIsPeriodic("a","-0.5","+0.5");
216 18 : av->addComponentWithDerivatives("b"); av->componentIsPeriodic("b","-0.5","+0.5");
217 18 : av->addComponentWithDerivatives("c"); av->componentIsPeriodic("c","-0.5","+0.5");
218 6 : return 2;
219 : }
220 1286 : av->addValueWithDerivatives(); av->setNotPeriodic();
221 : return 0;
222 : }
223 :
224 : // calculator
225 74406 : void Distance::calculate() {
226 :
227 74406 : if(pbc) makeWhole();
228 :
229 74406 : if( components ) {
230 409 : calculateCV( 1, masses, charges, getPositions(), value, derivs, virial, this );
231 409 : Value* valuex=getPntrToComponent("x");
232 409 : Value* valuey=getPntrToComponent("y");
233 409 : Value* valuez=getPntrToComponent("z");
234 :
235 1227 : for(unsigned i=0; i<2; ++i) setAtomsDerivatives(valuex,i,derivs[0][i] );
236 409 : setBoxDerivatives(valuex,virial[0]);
237 409 : valuex->set(value[0]);
238 :
239 1227 : for(unsigned i=0; i<2; ++i) setAtomsDerivatives(valuey,i,derivs[1][i] );
240 409 : setBoxDerivatives(valuey,virial[1]);
241 409 : valuey->set(value[1]);
242 :
243 1227 : for(unsigned i=0; i<2; ++i) setAtomsDerivatives(valuez,i,derivs[2][i] );
244 409 : setBoxDerivatives(valuez,virial[2]);
245 409 : valuez->set(value[2]);
246 73997 : } else if( scaled_components ) {
247 26 : calculateCV( 2, masses, charges, getPositions(), value, derivs, virial, this );
248 :
249 26 : Value* valuea=getPntrToComponent("a");
250 26 : Value* valueb=getPntrToComponent("b");
251 26 : Value* valuec=getPntrToComponent("c");
252 78 : for(unsigned i=0; i<2; ++i) setAtomsDerivatives(valuea,i,derivs[0][i] );
253 26 : valuea->set(value[0]);
254 78 : for(unsigned i=0; i<2; ++i) setAtomsDerivatives(valueb,i,derivs[1][i] );
255 26 : valueb->set(value[1]);
256 78 : for(unsigned i=0; i<2; ++i) setAtomsDerivatives(valuec,i,derivs[2][i] );
257 26 : valuec->set(value[2]);
258 : } else {
259 73971 : calculateCV( 0, masses, charges, getPositions(), value, derivs, virial, this );
260 221913 : for(unsigned i=0; i<2; ++i) setAtomsDerivatives(i,derivs[0][i] );
261 73971 : setBoxDerivatives(virial[0]);
262 73971 : setValue (value[0]);
263 : }
264 74406 : }
265 :
266 406099 : void Distance::calculateCV( const unsigned& mode, const std::vector<double>& masses, const std::vector<double>& charges,
267 : const std::vector<Vector>& pos, std::vector<double>& vals, std::vector<std::vector<Vector> >& derivs,
268 : std::vector<Tensor>& virial, const ActionAtomistic* aa ) {
269 406099 : Vector distance=delta(pos[0],pos[1]);
270 406099 : const double value=distance.modulo();
271 406099 : const double invvalue=1.0/value;
272 :
273 406099 : if(mode==1) {
274 94377 : derivs[0][0] = Vector(-1,0,0);
275 94377 : derivs[0][1] = Vector(+1,0,0);
276 94377 : vals[0] = distance[0];
277 :
278 94377 : derivs[1][0] = Vector(0,-1,0);
279 94377 : derivs[1][1] = Vector(0,+1,0);
280 94377 : vals[1] = distance[1];
281 :
282 94377 : derivs[2][0] = Vector(0,0,-1);
283 94377 : derivs[2][1] = Vector(0,0,+1);
284 94377 : vals[2] = distance[2];
285 94377 : setBoxDerivativesNoPbc( pos, derivs, virial );
286 311722 : } else if(mode==2) {
287 41 : Vector d=aa->getPbc().realToScaled(distance);
288 41 : derivs[0][0] = matmul(aa->getPbc().getInvBox(),Vector(-1,0,0));
289 41 : derivs[0][1] = matmul(aa->getPbc().getInvBox(),Vector(+1,0,0));
290 41 : vals[0] = Tools::pbc(d[0]);
291 41 : derivs[1][0] = matmul(aa->getPbc().getInvBox(),Vector(0,-1,0));
292 41 : derivs[1][1] = matmul(aa->getPbc().getInvBox(),Vector(0,+1,0));
293 41 : vals[1] = Tools::pbc(d[1]);
294 41 : derivs[2][0] = matmul(aa->getPbc().getInvBox(),Vector(0,0,-1));
295 41 : derivs[2][1] = matmul(aa->getPbc().getInvBox(),Vector(0,0,+1));
296 41 : vals[2] = Tools::pbc(d[2]);
297 : } else {
298 311681 : derivs[0][0] = -invvalue*distance;
299 311681 : derivs[0][1] = invvalue*distance;
300 311681 : setBoxDerivativesNoPbc( pos, derivs, virial );
301 311681 : vals[0] = value;
302 : }
303 406099 : }
304 :
305 : }
306 : }
307 :
308 :
309 :
|
