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/s between pairs of atoms.
34 :
35 : The following example illustrates how this action can be used to calculate and print the distance between atom 1
36 : and atom 2.
37 :
38 : ```plumed
39 : d: DISTANCE ATOMS=1,2
40 : PRINT ARG=d FILE=colvar
41 : ```
42 :
43 : By default the distance is computed in a way that takes periodic
44 : boundary conditions in account. This behavior can be changed by using the NOPBC flag.
45 : Furthermore, if you wish to calculate the vector connecting a pair of atoms you can use the
46 : `COMPONENTS` flag as shown below:
47 :
48 : ```plumed
49 : d: DISTANCE ATOMS=1,2 COMPONENTS
50 : PRINT ARG=d.x,d.y,d.z FILE=colvar
51 : ```
52 :
53 : Alternatively, you can calculate the components projected on the lattice vector by using the `SCALED_COMPONENTS`
54 : flag as shown below;
55 :
56 : ```plumed
57 : d: DISTANCE ATOMS=1,2 SCALED_COMPONENTS
58 : PRINT ARG=d.a,d.b,d.c FILE=colvar
59 : ```
60 :
61 : The advantage of using `SCALED_COMPONENTS` over `COMPONENTS` is that the a, b and c variables
62 : that are calculated when `SCALED_COMPONENTS` is employed have the proper periodicity. This feature is useful
63 : if you wish to study the motion of a molecule across a membrane.
64 :
65 : You can also use this command to calculate multiple indistinguishable distances or vectors with a single
66 : line of PLUMED input. For example, the following input calculates and outputs the distances between four
67 : pairs of atoms:
68 :
69 : ```plumed
70 : d: DISTANCE ATOMS1=1,2 ATOMS2=3,4 ATOMS3=5,6 ATOMS4=7,8
71 : PRINT ARG=d FILE=colvar
72 : ```
73 :
74 : By a similar token, the following input outputs three four dimensional vectors that contain the x, y and z
75 : components of the vectors connecting the four atoms:
76 :
77 : ```plumed
78 : d: DISTANCE COMPONENTS ATOMS1=1,2 ATOMS2=3,4 ATOMS3=5,6 ATOMS4=7,8
79 : PRINT ARG=d.x,d.y,d.z FILE=colvar
80 : ```
81 :
82 : You can also replace COMPONENTS with SCALED_COMPONENTS in the above input and obtain the projects of these vectors
83 : on the lattice vectors.
84 :
85 : ## Managing periodic boundary conditions
86 :
87 : When using the DISTANCE command to calculate the end-to-end distance for a large polymer you need to ensure that you
88 : are managing PBCs correctly. This problems that can occur with these calculations are explained at length in the
89 : early parts of the document that is referenced in the bibliography. Notice, however, that the input
90 : provides an example of an input that could be used to compute the end-to-end distance for a polymer
91 : of 100 atoms and keeps it at a value around 5.
92 :
93 : ```plumed
94 : WHOLEMOLECULES ENTITY0=1-100
95 : e2e: DISTANCE ATOMS=1,100 NOPBC
96 : RESTRAINT ARG=e2e KAPPA=1 AT=5
97 : ```
98 :
99 : Notice that NOPBC is used here so as to ensure that the distance is calculated correctely even if the end-to-end distance is larger than half the simulation
100 : Also notice that, since many MD codes break molecules across cell boundary, it might be necessary to
101 : use the [WHOLEMOLECULES](WHOLEMOLECULES.md) keyword (also notice that it should be _before_ distance). The list of atoms provided to [WHOLEMOLECULES](WHOLEMOLECULES.md)
102 : here contains all the atoms between 1 and 100. Strictly speaking, this
103 : is not necessary. If you know for sure that atoms with difference in
104 : the index say equal to 10 are _not_ going to be farther than half cell
105 : you can e.g. use
106 :
107 : ```plumed
108 : WHOLEMOLECULES ENTITY0=1,10,20,30,40,50,60,70,80,90,100
109 : e2e: DISTANCE ATOMS=1,100 NOPBC
110 : RESTRAINT ARG=e2e KAPPA=1 AT=5
111 : ```
112 :
113 : Just be sure that the ordered list provide to [WHOLEMOLECULES](WHOLEMOLECULES.md) has the following
114 : properties:
115 :
116 : - Consecutive atoms should be closer than half-cell throughout the entire simulation.
117 : - Atoms required later for the distance (e.g. 1 and 100) should be included in the list
118 :
119 : The following example shows how to take periodicity into account when computing the z-component of a distance
120 :
121 : ```plumed
122 : # this is a center of mass of a large group
123 : c: COM ATOMS=1-100
124 : # this is the distance between atom 101 and the group
125 : d: DISTANCE ATOMS=c,101 COMPONENTS
126 : # this makes a new variable, dd, equal to d and periodic, with domain -10,10
127 : # this is the right choise if e.g. the cell is orthorombic and its size in
128 : # z direction is 20.
129 : dz: COMBINE ARG=d.z PERIODIC=-10,10
130 : # metadynamics on dd
131 : METAD ARG=dz SIGMA=0.1 HEIGHT=0.1 PACE=200
132 : ```
133 :
134 : You can use the same input even if the DISTANCE command is calculating the vectors connecting multiple pairs of atoms.
135 : However, using SCALED_COMPONENTS ensures this problem does not arise because these variables are always periodic
136 : with domain (-0.5,+0.5).
137 :
138 : */
139 : //+ENDPLUMEDOC
140 :
141 : //+PLUMEDOC MCOLVAR DISTANCE_SCALAR
142 : /*
143 : Calculate the distance between a pair of atoms
144 :
145 : \par Examples
146 :
147 : */
148 : //+ENDPLUMEDOC
149 :
150 : //+PLUMEDOC MCOLVAR DISTANCE_VECTOR
151 : /*
152 : Calculate a vector containing the distances between various pairs of atoms
153 :
154 : \par Examples
155 :
156 : */
157 : //+ENDPLUMEDOC
158 :
159 : class Distance : public Colvar {
160 : bool components;
161 : bool scaled_components;
162 : bool pbc;
163 :
164 : std::vector<double> value, masses, charges;
165 : std::vector<std::vector<Vector> > derivs;
166 : std::vector<Tensor> virial;
167 : public:
168 : static void registerKeywords( Keywords& keys );
169 : explicit Distance(const ActionOptions&);
170 : static void parseAtomList( const int& num, std::vector<AtomNumber>& t, ActionAtomistic* aa );
171 : static unsigned getModeAndSetupValues( ActionWithValue* av );
172 : // active methods:
173 : void calculate() override;
174 : static void calculateCV( const unsigned& mode, const std::vector<double>& masses, const std::vector<double>& charges,
175 : const std::vector<Vector>& pos, std::vector<double>& vals, std::vector<std::vector<Vector> >& derivs,
176 : std::vector<Tensor>& virial, const ActionAtomistic* aa );
177 : };
178 :
179 : typedef ColvarShortcut<Distance> DistanceShortcut;
180 : PLUMED_REGISTER_ACTION(DistanceShortcut,"DISTANCE")
181 : PLUMED_REGISTER_ACTION(Distance,"DISTANCE_SCALAR")
182 : typedef MultiColvarTemplate<Distance> DistanceMulti;
183 : PLUMED_REGISTER_ACTION(DistanceMulti,"DISTANCE_VECTOR")
184 :
185 2380 : void Distance::registerKeywords( Keywords& keys ) {
186 2380 : Colvar::registerKeywords( keys );
187 4760 : keys.setDisplayName("DISTANCE");
188 : constexpr auto scalarOrVector = Keywords::componentType::scalar | Keywords::componentType::vector;
189 2380 : keys.add("atoms","ATOMS","the pair of atom that we are calculating the distance between");
190 2380 : 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");
191 2380 : 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");
192 4760 : keys.addOutputComponent("x","COMPONENTS",scalarOrVector,"the x-component of the vector connecting the two atoms");
193 4760 : keys.addOutputComponent("y","COMPONENTS",scalarOrVector,"the y-component of the vector connecting the two atoms");
194 4760 : keys.addOutputComponent("z","COMPONENTS",scalarOrVector,"the z-component of the vector connecting the two atoms");
195 4760 : keys.addOutputComponent("a","SCALED_COMPONENTS",scalarOrVector,"the normalized projection on the first lattice vector of the vector connecting the two atoms");
196 4760 : keys.addOutputComponent("b","SCALED_COMPONENTS",scalarOrVector,"the normalized projection on the second lattice vector of the vector connecting the two atoms");
197 4760 : keys.addOutputComponent("c","SCALED_COMPONENTS",scalarOrVector,"the normalized projection on the third lattice vector of the vector connecting the two atoms");
198 2380 : keys.add("hidden","NO_ACTION_LOG","suppresses printing from action on the log");
199 2380 : keys.setValueDescription(scalarOrVector,"the DISTANCE between this pair of atoms");
200 2380 : keys.addDOI("10.48550/arXiv.1812.08213");
201 2380 : keys.addDOI("10.1007/978-1-4939-9608-7_21");
202 2380 : }
203 :
204 623 : Distance::Distance(const ActionOptions&ao):
205 : PLUMED_COLVAR_INIT(ao),
206 623 : components(false),
207 623 : scaled_components(false),
208 623 : pbc(true),
209 623 : value(1),
210 625 : derivs(1),
211 1246 : virial(1) {
212 623 : derivs[0].resize(2);
213 : std::vector<AtomNumber> atoms;
214 623 : parseAtomList(-1,atoms,this);
215 623 : if(atoms.size()!=2) {
216 1 : error("Number of specified atoms should be 2");
217 : }
218 :
219 622 : bool nopbc=!pbc;
220 624 : parseFlag("NOPBC",nopbc);
221 622 : pbc=!nopbc;
222 :
223 622 : if(pbc) {
224 618 : log.printf(" using periodic boundary conditions\n");
225 : } else {
226 4 : log.printf(" without periodic boundary conditions\n");
227 : }
228 :
229 622 : unsigned mode = getModeAndSetupValues( this );
230 621 : if(mode==1) {
231 35 : components=true;
232 586 : } else if(mode==2) {
233 5 : scaled_components=true;
234 : }
235 621 : if( components || scaled_components ) {
236 40 : value.resize(3);
237 40 : derivs.resize(3);
238 40 : virial.resize(3);
239 160 : for(unsigned i=0; i<3; ++i) {
240 120 : derivs[i].resize(2);
241 : }
242 : }
243 621 : requestAtoms(atoms);
244 627 : }
245 :
246 30941 : void Distance::parseAtomList( const int& num, std::vector<AtomNumber>& t, ActionAtomistic* aa ) {
247 61882 : aa->parseAtomList("ATOMS",num,t);
248 30941 : if( t.size()==2 ) {
249 30832 : aa->log.printf(" between atoms %d %d\n",t[0].serial(),t[1].serial());
250 : }
251 30941 : }
252 :
253 730 : unsigned Distance::getModeAndSetupValues( ActionWithValue* av ) {
254 : bool c;
255 730 : av->parseFlag("COMPONENTS",c);
256 : bool sc;
257 730 : av->parseFlag("SCALED_COMPONENTS",sc);
258 730 : if( c && sc ) {
259 1 : av->error("COMPONENTS and SCALED_COMPONENTS are not compatible");
260 : }
261 :
262 729 : if(c) {
263 160 : av->addComponentWithDerivatives("x");
264 80 : av->componentIsNotPeriodic("x");
265 160 : av->addComponentWithDerivatives("y");
266 80 : av->componentIsNotPeriodic("y");
267 160 : av->addComponentWithDerivatives("z");
268 80 : av->componentIsNotPeriodic("z");
269 80 : av->log<<" WARNING: components will not have the proper periodicity - see manual\n";
270 80 : return 1;
271 649 : } else if(sc) {
272 12 : av->addComponentWithDerivatives("a");
273 12 : av->componentIsPeriodic("a","-0.5","+0.5");
274 12 : av->addComponentWithDerivatives("b");
275 12 : av->componentIsPeriodic("b","-0.5","+0.5");
276 12 : av->addComponentWithDerivatives("c");
277 12 : av->componentIsPeriodic("c","-0.5","+0.5");
278 6 : return 2;
279 : }
280 643 : av->addValueWithDerivatives();
281 643 : av->setNotPeriodic();
282 : return 0;
283 : }
284 :
285 : // calculator
286 74406 : void Distance::calculate() {
287 :
288 74406 : if(pbc) {
289 74386 : makeWhole();
290 : }
291 :
292 74406 : if( components ) {
293 409 : calculateCV( 1, masses, charges, getPositions(), value, derivs, virial, this );
294 409 : Value* valuex=getPntrToComponent("x");
295 409 : Value* valuey=getPntrToComponent("y");
296 409 : Value* valuez=getPntrToComponent("z");
297 :
298 1227 : for(unsigned i=0; i<2; ++i) {
299 818 : setAtomsDerivatives(valuex,i,derivs[0][i] );
300 : }
301 409 : setBoxDerivatives(valuex,virial[0]);
302 409 : valuex->set(value[0]);
303 :
304 1227 : for(unsigned i=0; i<2; ++i) {
305 818 : setAtomsDerivatives(valuey,i,derivs[1][i] );
306 : }
307 409 : setBoxDerivatives(valuey,virial[1]);
308 409 : valuey->set(value[1]);
309 :
310 1227 : for(unsigned i=0; i<2; ++i) {
311 818 : setAtomsDerivatives(valuez,i,derivs[2][i] );
312 : }
313 409 : setBoxDerivatives(valuez,virial[2]);
314 409 : valuez->set(value[2]);
315 73997 : } else if( scaled_components ) {
316 26 : calculateCV( 2, masses, charges, getPositions(), value, derivs, virial, this );
317 :
318 26 : Value* valuea=getPntrToComponent("a");
319 26 : Value* valueb=getPntrToComponent("b");
320 26 : Value* valuec=getPntrToComponent("c");
321 78 : for(unsigned i=0; i<2; ++i) {
322 52 : setAtomsDerivatives(valuea,i,derivs[0][i] );
323 : }
324 26 : valuea->set(value[0]);
325 78 : for(unsigned i=0; i<2; ++i) {
326 52 : setAtomsDerivatives(valueb,i,derivs[1][i] );
327 : }
328 26 : valueb->set(value[1]);
329 78 : for(unsigned i=0; i<2; ++i) {
330 52 : setAtomsDerivatives(valuec,i,derivs[2][i] );
331 : }
332 26 : valuec->set(value[2]);
333 : } else {
334 73971 : calculateCV( 0, masses, charges, getPositions(), value, derivs, virial, this );
335 221913 : for(unsigned i=0; i<2; ++i) {
336 147942 : setAtomsDerivatives(i,derivs[0][i] );
337 : }
338 73971 : setBoxDerivatives(virial[0]);
339 73971 : setValue (value[0]);
340 : }
341 74406 : }
342 :
343 406099 : void Distance::calculateCV( const unsigned& mode, const std::vector<double>& masses, const std::vector<double>& charges,
344 : const std::vector<Vector>& pos, std::vector<double>& vals, std::vector<std::vector<Vector> >& derivs,
345 : std::vector<Tensor>& virial, const ActionAtomistic* aa ) {
346 406099 : Vector distance=delta(pos[0],pos[1]);
347 406099 : const double value=distance.modulo();
348 406099 : const double invvalue=1.0/value;
349 :
350 406099 : if(mode==1) {
351 94377 : derivs[0][0] = Vector(-1,0,0);
352 94377 : derivs[0][1] = Vector(+1,0,0);
353 94377 : vals[0] = distance[0];
354 :
355 94377 : derivs[1][0] = Vector(0,-1,0);
356 94377 : derivs[1][1] = Vector(0,+1,0);
357 94377 : vals[1] = distance[1];
358 :
359 94377 : derivs[2][0] = Vector(0,0,-1);
360 94377 : derivs[2][1] = Vector(0,0,+1);
361 94377 : vals[2] = distance[2];
362 94377 : setBoxDerivativesNoPbc( pos, derivs, virial );
363 311722 : } else if(mode==2) {
364 41 : Vector d=aa->getPbc().realToScaled(distance);
365 41 : derivs[0][0] = matmul(aa->getPbc().getInvBox(),Vector(-1,0,0));
366 41 : derivs[0][1] = matmul(aa->getPbc().getInvBox(),Vector(+1,0,0));
367 41 : vals[0] = Tools::pbc(d[0]);
368 41 : derivs[1][0] = matmul(aa->getPbc().getInvBox(),Vector(0,-1,0));
369 41 : derivs[1][1] = matmul(aa->getPbc().getInvBox(),Vector(0,+1,0));
370 41 : vals[1] = Tools::pbc(d[1]);
371 41 : derivs[2][0] = matmul(aa->getPbc().getInvBox(),Vector(0,0,-1));
372 41 : derivs[2][1] = matmul(aa->getPbc().getInvBox(),Vector(0,0,+1));
373 41 : vals[2] = Tools::pbc(d[2]);
374 : } else {
375 311681 : derivs[0][0] = -invvalue*distance;
376 311681 : derivs[0][1] = invvalue*distance;
377 311681 : setBoxDerivativesNoPbc( pos, derivs, virial );
378 311681 : vals[0] = value;
379 : }
380 406099 : }
381 :
382 : }
383 : }
384 :
385 :
386 :
|
