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 11327 : PLUMED_REGISTER_ACTION(Distance,"DISTANCE")
119 :
120 456 : void Distance::registerKeywords( Keywords& keys ) {
121 456 : Colvar::registerKeywords( keys );
122 912 : keys.add("atoms","ATOMS","the pair of atom that we are calculating the distance between");
123 912 : 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 912 : 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 912 : keys.addOutputComponent("x","COMPONENTS","the x-component of the vector connecting the two atoms");
126 912 : keys.addOutputComponent("y","COMPONENTS","the y-component of the vector connecting the two atoms");
127 912 : keys.addOutputComponent("z","COMPONENTS","the z-component of the vector connecting the two atoms");
128 1368 : keys.addOutputComponent("a","SCALED_COMPONENTS","the normalized projection on the first lattice vector of the vector connecting the two atoms");
129 1368 : keys.addOutputComponent("b","SCALED_COMPONENTS","the normalized projection on the second lattice vector of the vector connecting the two atoms");
130 1368 : keys.addOutputComponent("c","SCALED_COMPONENTS","the normalized projection on the third lattice vector of the vector connecting the two atoms");
131 456 : }
132 :
133 455 : Distance::Distance(const ActionOptions&ao):
134 : PLUMED_COLVAR_INIT(ao),
135 455 : components(false),
136 455 : scaled_components(false),
137 455 : pbc(true)
138 : {
139 : std::vector<AtomNumber> atoms;
140 910 : parseAtomList("ATOMS",atoms);
141 455 : if(atoms.size()!=2)
142 1 : error("Number of specified atoms should be 2");
143 454 : parseFlag("COMPONENTS",components);
144 454 : parseFlag("SCALED_COMPONENTS",scaled_components);
145 454 : bool nopbc=!pbc;
146 454 : parseFlag("NOPBC",nopbc);
147 454 : pbc=!nopbc;
148 454 : checkRead();
149 :
150 454 : log.printf(" between atoms %d %d\n",atoms[0].serial(),atoms[1].serial());
151 454 : if(pbc) log.printf(" using periodic boundary conditions\n");
152 4 : else log.printf(" without periodic boundary conditions\n");
153 :
154 454 : if(components && scaled_components) error("COMPONENTS and SCALED_COMPONENTS are not compatible");
155 :
156 453 : if(components) {
157 62 : addComponentWithDerivatives("x"); componentIsNotPeriodic("x");
158 62 : addComponentWithDerivatives("y"); componentIsNotPeriodic("y");
159 62 : addComponentWithDerivatives("z"); componentIsNotPeriodic("z");
160 31 : log<<" WARNING: components will not have the proper periodicity - see manual\n";
161 422 : } else if(scaled_components) {
162 6 : addComponentWithDerivatives("a"); componentIsPeriodic("a","-0.5","+0.5");
163 6 : addComponentWithDerivatives("b"); componentIsPeriodic("b","-0.5","+0.5");
164 8 : addComponentWithDerivatives("c"); componentIsPeriodic("c","-0.5","+0.5");
165 : } else {
166 420 : addValueWithDerivatives(); setNotPeriodic();
167 : }
168 :
169 :
170 453 : requestAtoms(atoms);
171 457 : }
172 :
173 :
174 : // calculator
175 42339 : void Distance::calculate() {
176 :
177 42339 : if(pbc) makeWhole();
178 :
179 42339 : Vector distance=delta(getPosition(0),getPosition(1));
180 42339 : const double value=distance.modulo();
181 42339 : const double invvalue=1.0/value;
182 :
183 42339 : if(components) {
184 394 : Value* valuex=getPntrToComponent("x");
185 394 : Value* valuey=getPntrToComponent("y");
186 394 : Value* valuez=getPntrToComponent("z");
187 :
188 394 : setAtomsDerivatives (valuex,0,Vector(-1,0,0));
189 394 : setAtomsDerivatives (valuex,1,Vector(+1,0,0));
190 394 : setBoxDerivativesNoPbc(valuex);
191 394 : valuex->set(distance[0]);
192 :
193 394 : setAtomsDerivatives (valuey,0,Vector(0,-1,0));
194 394 : setAtomsDerivatives (valuey,1,Vector(0,+1,0));
195 394 : setBoxDerivativesNoPbc(valuey);
196 394 : valuey->set(distance[1]);
197 :
198 394 : setAtomsDerivatives (valuez,0,Vector(0,0,-1));
199 394 : setAtomsDerivatives (valuez,1,Vector(0,0,+1));
200 394 : setBoxDerivativesNoPbc(valuez);
201 394 : valuez->set(distance[2]);
202 41945 : } else if(scaled_components) {
203 11 : Value* valuea=getPntrToComponent("a");
204 11 : Value* valueb=getPntrToComponent("b");
205 22 : Value* valuec=getPntrToComponent("c");
206 11 : Vector d=getPbc().realToScaled(distance);
207 11 : setAtomsDerivatives (valuea,0,matmul(getPbc().getInvBox(),Vector(-1,0,0)));
208 11 : setAtomsDerivatives (valuea,1,matmul(getPbc().getInvBox(),Vector(+1,0,0)));
209 11 : valuea->set(Tools::pbc(d[0]));
210 11 : setAtomsDerivatives (valueb,0,matmul(getPbc().getInvBox(),Vector(0,-1,0)));
211 11 : setAtomsDerivatives (valueb,1,matmul(getPbc().getInvBox(),Vector(0,+1,0)));
212 11 : valueb->set(Tools::pbc(d[1]));
213 11 : setAtomsDerivatives (valuec,0,matmul(getPbc().getInvBox(),Vector(0,0,-1)));
214 11 : setAtomsDerivatives (valuec,1,matmul(getPbc().getInvBox(),Vector(0,0,+1)));
215 11 : valuec->set(Tools::pbc(d[2]));
216 : } else {
217 41934 : setAtomsDerivatives(0,-invvalue*distance);
218 41934 : setAtomsDerivatives(1,invvalue*distance);
219 41934 : setBoxDerivativesNoPbc();
220 41934 : setValue (value);
221 : }
222 :
223 42339 : }
224 :
225 : }
226 : }
227 :
228 :
229 :
|
