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