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 );
154 2380 : keys.setDisplayName("DISTANCE");
155 4760 : keys.add("atoms","ATOMS","the pair of atom that we are calculating the distance between");
156 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");
157 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");
158 4760 : keys.addOutputComponent("x","COMPONENTS","the x-component of the vector connecting the two atoms");
159 4760 : keys.addOutputComponent("y","COMPONENTS","the y-component of the vector connecting the two atoms");
160 4760 : keys.addOutputComponent("z","COMPONENTS","the z-component of the vector connecting the two atoms");
161 4760 : keys.addOutputComponent("a","SCALED_COMPONENTS","the normalized projection on the first lattice vector of the vector connecting the two atoms");
162 4760 : keys.addOutputComponent("b","SCALED_COMPONENTS","the normalized projection on the second lattice vector of the vector connecting the two atoms");
163 4760 : keys.addOutputComponent("c","SCALED_COMPONENTS","the normalized projection on the third lattice vector of the vector connecting the two atoms");
164 4760 : keys.add("hidden","NO_ACTION_LOG","suppresses printing from action on the log");
165 2380 : keys.setValueDescription("the DISTANCE between this pair of atoms");
166 2380 : }
167 :
168 623 : Distance::Distance(const ActionOptions&ao):
169 : PLUMED_COLVAR_INIT(ao),
170 623 : components(false),
171 623 : scaled_components(false),
172 623 : pbc(true),
173 623 : value(1),
174 625 : derivs(1),
175 1246 : virial(1) {
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 :
183 622 : bool nopbc=!pbc;
184 624 : parseFlag("NOPBC",nopbc);
185 622 : pbc=!nopbc;
186 :
187 622 : if(pbc) {
188 618 : log.printf(" using periodic boundary conditions\n");
189 : } else {
190 4 : log.printf(" without periodic boundary conditions\n");
191 : }
192 :
193 622 : unsigned mode = getModeAndSetupValues( this );
194 621 : if(mode==1) {
195 35 : components=true;
196 586 : } else if(mode==2) {
197 5 : scaled_components=true;
198 : }
199 621 : if( components || scaled_components ) {
200 40 : value.resize(3);
201 40 : derivs.resize(3);
202 40 : virial.resize(3);
203 160 : for(unsigned i=0; i<3; ++i) {
204 120 : derivs[i].resize(2);
205 : }
206 : }
207 621 : requestAtoms(atoms);
208 627 : }
209 :
210 30941 : void Distance::parseAtomList( const int& num, std::vector<AtomNumber>& t, ActionAtomistic* aa ) {
211 61882 : aa->parseAtomList("ATOMS",num,t);
212 30941 : if( t.size()==2 ) {
213 30832 : aa->log.printf(" between atoms %d %d\n",t[0].serial(),t[1].serial());
214 : }
215 30941 : }
216 :
217 730 : unsigned Distance::getModeAndSetupValues( ActionWithValue* av ) {
218 : bool c;
219 730 : av->parseFlag("COMPONENTS",c);
220 : bool sc;
221 730 : av->parseFlag("SCALED_COMPONENTS",sc);
222 730 : if( c && sc ) {
223 1 : av->error("COMPONENTS and SCALED_COMPONENTS are not compatible");
224 : }
225 :
226 729 : if(c) {
227 160 : av->addComponentWithDerivatives("x");
228 80 : av->componentIsNotPeriodic("x");
229 160 : av->addComponentWithDerivatives("y");
230 80 : av->componentIsNotPeriodic("y");
231 160 : av->addComponentWithDerivatives("z");
232 80 : av->componentIsNotPeriodic("z");
233 80 : av->log<<" WARNING: components will not have the proper periodicity - see manual\n";
234 80 : return 1;
235 649 : } else if(sc) {
236 12 : av->addComponentWithDerivatives("a");
237 12 : av->componentIsPeriodic("a","-0.5","+0.5");
238 12 : av->addComponentWithDerivatives("b");
239 12 : av->componentIsPeriodic("b","-0.5","+0.5");
240 12 : av->addComponentWithDerivatives("c");
241 12 : av->componentIsPeriodic("c","-0.5","+0.5");
242 6 : return 2;
243 : }
244 643 : av->addValueWithDerivatives();
245 643 : av->setNotPeriodic();
246 : return 0;
247 : }
248 :
249 : // calculator
250 74406 : void Distance::calculate() {
251 :
252 74406 : if(pbc) {
253 74386 : makeWhole();
254 : }
255 :
256 74406 : if( components ) {
257 409 : calculateCV( 1, masses, charges, getPositions(), value, derivs, virial, this );
258 409 : Value* valuex=getPntrToComponent("x");
259 409 : Value* valuey=getPntrToComponent("y");
260 409 : Value* valuez=getPntrToComponent("z");
261 :
262 1227 : for(unsigned i=0; i<2; ++i) {
263 818 : setAtomsDerivatives(valuex,i,derivs[0][i] );
264 : }
265 409 : setBoxDerivatives(valuex,virial[0]);
266 409 : valuex->set(value[0]);
267 :
268 1227 : for(unsigned i=0; i<2; ++i) {
269 818 : setAtomsDerivatives(valuey,i,derivs[1][i] );
270 : }
271 409 : setBoxDerivatives(valuey,virial[1]);
272 409 : valuey->set(value[1]);
273 :
274 1227 : for(unsigned i=0; i<2; ++i) {
275 818 : setAtomsDerivatives(valuez,i,derivs[2][i] );
276 : }
277 409 : setBoxDerivatives(valuez,virial[2]);
278 409 : valuez->set(value[2]);
279 73997 : } else if( scaled_components ) {
280 26 : calculateCV( 2, masses, charges, getPositions(), value, derivs, virial, this );
281 :
282 26 : Value* valuea=getPntrToComponent("a");
283 26 : Value* valueb=getPntrToComponent("b");
284 26 : Value* valuec=getPntrToComponent("c");
285 78 : for(unsigned i=0; i<2; ++i) {
286 52 : setAtomsDerivatives(valuea,i,derivs[0][i] );
287 : }
288 26 : valuea->set(value[0]);
289 78 : for(unsigned i=0; i<2; ++i) {
290 52 : setAtomsDerivatives(valueb,i,derivs[1][i] );
291 : }
292 26 : valueb->set(value[1]);
293 78 : for(unsigned i=0; i<2; ++i) {
294 52 : setAtomsDerivatives(valuec,i,derivs[2][i] );
295 : }
296 26 : valuec->set(value[2]);
297 : } else {
298 73971 : calculateCV( 0, masses, charges, getPositions(), value, derivs, virial, this );
299 221913 : for(unsigned i=0; i<2; ++i) {
300 147942 : setAtomsDerivatives(i,derivs[0][i] );
301 : }
302 73971 : setBoxDerivatives(virial[0]);
303 73971 : setValue (value[0]);
304 : }
305 74406 : }
306 :
307 406099 : void Distance::calculateCV( const unsigned& mode, const std::vector<double>& masses, const std::vector<double>& charges,
308 : const std::vector<Vector>& pos, std::vector<double>& vals, std::vector<std::vector<Vector> >& derivs,
309 : std::vector<Tensor>& virial, const ActionAtomistic* aa ) {
310 406099 : Vector distance=delta(pos[0],pos[1]);
311 406099 : const double value=distance.modulo();
312 406099 : const double invvalue=1.0/value;
313 :
314 406099 : if(mode==1) {
315 94377 : derivs[0][0] = Vector(-1,0,0);
316 94377 : derivs[0][1] = Vector(+1,0,0);
317 94377 : vals[0] = distance[0];
318 :
319 94377 : derivs[1][0] = Vector(0,-1,0);
320 94377 : derivs[1][1] = Vector(0,+1,0);
321 94377 : vals[1] = distance[1];
322 :
323 94377 : derivs[2][0] = Vector(0,0,-1);
324 94377 : derivs[2][1] = Vector(0,0,+1);
325 94377 : vals[2] = distance[2];
326 94377 : setBoxDerivativesNoPbc( pos, derivs, virial );
327 311722 : } else if(mode==2) {
328 41 : Vector d=aa->getPbc().realToScaled(distance);
329 41 : derivs[0][0] = matmul(aa->getPbc().getInvBox(),Vector(-1,0,0));
330 41 : derivs[0][1] = matmul(aa->getPbc().getInvBox(),Vector(+1,0,0));
331 41 : vals[0] = Tools::pbc(d[0]);
332 41 : derivs[1][0] = matmul(aa->getPbc().getInvBox(),Vector(0,-1,0));
333 41 : derivs[1][1] = matmul(aa->getPbc().getInvBox(),Vector(0,+1,0));
334 41 : vals[1] = Tools::pbc(d[1]);
335 41 : derivs[2][0] = matmul(aa->getPbc().getInvBox(),Vector(0,0,-1));
336 41 : derivs[2][1] = matmul(aa->getPbc().getInvBox(),Vector(0,0,+1));
337 41 : vals[2] = Tools::pbc(d[2]);
338 : } else {
339 311681 : derivs[0][0] = -invvalue*distance;
340 311681 : derivs[0][1] = invvalue*distance;
341 311681 : setBoxDerivativesNoPbc( pos, derivs, virial );
342 311681 : vals[0] = value;
343 : }
344 406099 : }
345 :
346 : }
347 : }
348 :
349 :
350 :
|
