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