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 "ActionRegister.h"
24 : #include "tools/Pbc.h"
25 :
26 : namespace PLMD {
27 : namespace colvar {
28 :
29 : //+PLUMEDOC COLVAR DISTANCE
30 : /*
31 : Calculate the distance between a pair of atoms.
32 :
33 : By default the distance is computed taking into account periodic
34 : boundary conditions. This behavior can be changed with the NOPBC flag.
35 : Moreover, single components in Cartesian space (x,y, and z, with COMPONENTS)
36 : or single components projected to the three lattice vectors (a,b, and c, with SCALED_COMPONENTS)
37 : can be also computed.
38 :
39 : Notice that Cartesian components will not have the proper periodicity!
40 : If you have to study e.g. the permeation of a molecule across a membrane,
41 : better to use SCALED_COMPONENTS.
42 :
43 : \par Examples
44 :
45 : The following input tells plumed to print the distance between atoms 3 and 5,
46 : the distance between atoms 2 and 4 and the x component of the distance between atoms 2 and 4.
47 : \plumedfile
48 : d1: DISTANCE ATOMS=3,5
49 : d2: DISTANCE ATOMS=2,4
50 : d2c: DISTANCE ATOMS=2,4 COMPONENTS
51 : PRINT ARG=d1,d2,d2c.x
52 : \endplumedfile
53 :
54 : The following input computes the end-to-end distance for a polymer
55 : of 100 atoms and keeps it at a value around 5.
56 : \plumedfile
57 : WHOLEMOLECULES ENTITY0=1-100
58 : e2e: DISTANCE ATOMS=1,100 NOPBC
59 : RESTRAINT ARG=e2e KAPPA=1 AT=5
60 : \endplumedfile
61 :
62 : Notice that NOPBC is used
63 : to be sure that if the end-to-end distance is larger than half the simulation
64 : box the distance is compute properly. Also notice that, since many MD
65 : codes break molecules across cell boundary, it might be necessary to
66 : use the \ref WHOLEMOLECULES keyword (also notice that it should be
67 : _before_ distance). The list of atoms provided to \ref WHOLEMOLECULES
68 : here contains all the atoms between 1 and 100. Strictly speaking, this
69 : is not necessary. If you know for sure that atoms with difference in
70 : the index say equal to 10 are _not_ going to be farther than half cell
71 : you can e.g. use
72 : \plumedfile
73 : WHOLEMOLECULES ENTITY0=1,10,20,30,40,50,60,70,80,90,100
74 : e2e: DISTANCE ATOMS=1,100 NOPBC
75 : RESTRAINT ARG=e2e KAPPA=1 AT=5
76 : \endplumedfile
77 : Just be sure that the ordered list provide to \ref WHOLEMOLECULES has the following
78 : properties:
79 : - Consecutive atoms should be closer than half-cell throughout the entire simulation.
80 : - Atoms required later for the distance (e.g. 1 and 100) should be included in the list
81 :
82 : The following example shows how to take into account periodicity e.g.
83 : in z-component of a distance
84 : \plumedfile
85 : # this is a center of mass of a large group
86 : c: COM ATOMS=1-100
87 : # this is the distance between atom 101 and the group
88 : d: DISTANCE ATOMS=c,101 COMPONENTS
89 : # this makes a new variable, dd, equal to d and periodic, with domain -10,10
90 : # this is the right choise if e.g. the cell is orthorombic and its size in
91 : # z direction is 20.
92 : dz: COMBINE ARG=d.z PERIODIC=-10,10
93 : # metadynamics on dd
94 : METAD ARG=dz SIGMA=0.1 HEIGHT=0.1 PACE=200
95 : \endplumedfile
96 :
97 : Using SCALED_COMPONENTS this problem should not arise because they are always periodic
98 : with domain (-0.5,+0.5).
99 :
100 :
101 :
102 :
103 : */
104 : //+ENDPLUMEDOC
105 :
106 : class Distance : public Colvar {
107 : bool components;
108 : bool scaled_components;
109 : bool pbc;
110 :
111 : public:
112 : static void registerKeywords( Keywords& keys );
113 : explicit Distance(const ActionOptions&);
114 : // active methods:
115 : void calculate() override;
116 : };
117 :
118 14721 : PLUMED_REGISTER_ACTION(Distance,"DISTANCE")
119 :
120 473 : void Distance::registerKeywords( Keywords& keys ) {
121 473 : Colvar::registerKeywords( keys );
122 946 : keys.add("atoms","ATOMS","the pair of atom that we are calculating the distance between");
123 946 : 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");
124 946 : 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");
125 946 : keys.addOutputComponent("x","COMPONENTS","the x-component of the vector connecting the two atoms");
126 946 : keys.addOutputComponent("y","COMPONENTS","the y-component of the vector connecting the two atoms");
127 946 : keys.addOutputComponent("z","COMPONENTS","the z-component of the vector connecting the two atoms");
128 946 : keys.addOutputComponent("a","SCALED_COMPONENTS","the normalized projection on the first lattice vector of the vector connecting the two atoms");
129 946 : keys.addOutputComponent("b","SCALED_COMPONENTS","the normalized projection on the second lattice vector of the vector connecting the two atoms");
130 946 : keys.addOutputComponent("c","SCALED_COMPONENTS","the normalized projection on the third lattice vector of the vector connecting the two atoms");
131 473 : }
132 :
133 469 : Distance::Distance(const ActionOptions&ao):
134 : PLUMED_COLVAR_INIT(ao),
135 469 : components(false),
136 469 : scaled_components(false),
137 469 : pbc(true) {
138 : std::vector<AtomNumber> atoms;
139 938 : parseAtomList("ATOMS",atoms);
140 469 : if(atoms.size()!=2) {
141 1 : error("Number of specified atoms should be 2");
142 : }
143 468 : parseFlag("COMPONENTS",components);
144 468 : parseFlag("SCALED_COMPONENTS",scaled_components);
145 468 : bool nopbc=!pbc;
146 468 : parseFlag("NOPBC",nopbc);
147 468 : pbc=!nopbc;
148 468 : checkRead();
149 :
150 468 : log.printf(" between atoms %d %d\n",atoms[0].serial(),atoms[1].serial());
151 468 : if(pbc) {
152 464 : log.printf(" using periodic boundary conditions\n");
153 : } else {
154 4 : log.printf(" without periodic boundary conditions\n");
155 : }
156 :
157 468 : if(components && scaled_components) {
158 1 : error("COMPONENTS and SCALED_COMPONENTS are not compatible");
159 : }
160 :
161 467 : if(components) {
162 31 : addComponentWithDerivatives("x");
163 31 : componentIsNotPeriodic("x");
164 31 : addComponentWithDerivatives("y");
165 31 : componentIsNotPeriodic("y");
166 31 : addComponentWithDerivatives("z");
167 31 : componentIsNotPeriodic("z");
168 31 : log<<" WARNING: components will not have the proper periodicity - see manual\n";
169 436 : } else if(scaled_components) {
170 2 : addComponentWithDerivatives("a");
171 4 : componentIsPeriodic("a","-0.5","+0.5");
172 2 : addComponentWithDerivatives("b");
173 4 : componentIsPeriodic("b","-0.5","+0.5");
174 2 : addComponentWithDerivatives("c");
175 6 : componentIsPeriodic("c","-0.5","+0.5");
176 : } else {
177 434 : addValueWithDerivatives();
178 434 : setNotPeriodic();
179 : }
180 :
181 :
182 467 : requestAtoms(atoms);
183 471 : }
184 :
185 :
186 : // calculator
187 42399 : void Distance::calculate() {
188 :
189 42399 : if(pbc) {
190 42379 : makeWhole();
191 : }
192 :
193 42399 : Vector distance=delta(getPosition(0),getPosition(1));
194 42399 : const double value=distance.modulo();
195 42399 : const double invvalue=1.0/value;
196 :
197 42399 : if(components) {
198 394 : Value* valuex=getPntrToComponent("x");
199 394 : Value* valuey=getPntrToComponent("y");
200 394 : Value* valuez=getPntrToComponent("z");
201 :
202 394 : setAtomsDerivatives (valuex,0,Vector(-1,0,0));
203 394 : setAtomsDerivatives (valuex,1,Vector(+1,0,0));
204 394 : setBoxDerivativesNoPbc(valuex);
205 394 : valuex->set(distance[0]);
206 :
207 394 : setAtomsDerivatives (valuey,0,Vector(0,-1,0));
208 394 : setAtomsDerivatives (valuey,1,Vector(0,+1,0));
209 394 : setBoxDerivativesNoPbc(valuey);
210 394 : valuey->set(distance[1]);
211 :
212 394 : setAtomsDerivatives (valuez,0,Vector(0,0,-1));
213 394 : setAtomsDerivatives (valuez,1,Vector(0,0,+1));
214 394 : setBoxDerivativesNoPbc(valuez);
215 394 : valuez->set(distance[2]);
216 42005 : } else if(scaled_components) {
217 11 : Value* valuea=getPntrToComponent("a");
218 11 : Value* valueb=getPntrToComponent("b");
219 22 : Value* valuec=getPntrToComponent("c");
220 11 : Vector d=getPbc().realToScaled(distance);
221 11 : setAtomsDerivatives (valuea,0,matmul(getPbc().getInvBox(),Vector(-1,0,0)));
222 11 : setAtomsDerivatives (valuea,1,matmul(getPbc().getInvBox(),Vector(+1,0,0)));
223 11 : valuea->set(Tools::pbc(d[0]));
224 11 : setAtomsDerivatives (valueb,0,matmul(getPbc().getInvBox(),Vector(0,-1,0)));
225 11 : setAtomsDerivatives (valueb,1,matmul(getPbc().getInvBox(),Vector(0,+1,0)));
226 11 : valueb->set(Tools::pbc(d[1]));
227 11 : setAtomsDerivatives (valuec,0,matmul(getPbc().getInvBox(),Vector(0,0,-1)));
228 11 : setAtomsDerivatives (valuec,1,matmul(getPbc().getInvBox(),Vector(0,0,+1)));
229 11 : valuec->set(Tools::pbc(d[2]));
230 : } else {
231 41994 : setAtomsDerivatives(0,-invvalue*distance);
232 41994 : setAtomsDerivatives(1,invvalue*distance);
233 41994 : setBoxDerivativesNoPbc();
234 41994 : setValue (value);
235 : }
236 :
237 42399 : }
238 :
239 : }
240 : }
241 :
242 :
243 :
|
