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 : #ifndef __PLUMED_colvar_Distance_h
23 : #define __PLUMED_colvar_Distance_h
24 : #include "Colvar.h"
25 : #include "ColvarInput.h"
26 : #include "MultiColvarTemplate.h"
27 : #include "tools/Pbc.h"
28 :
29 : namespace PLMD {
30 : namespace colvar {
31 :
32 : template <typename T=double>
33 : class Distance : public Colvar {
34 : enum components:unsigned {squared=3,scaled=2,standard=1,scalar=0} mode{scalar};
35 : bool pbc;
36 :
37 : //not Typep, the float version is only used for the PTM
38 : std::vector<double> value;
39 : std::vector<double> derivs;
40 : public:
41 : using precision = T;
42 : static void registerKeywords( Keywords& keys );
43 : explicit Distance(const ActionOptions&);
44 : static void parseAtomList( const int& num,
45 : std::vector<AtomNumber>& t,
46 : ActionAtomistic* aa );
47 : static unsigned getModeAndSetupValues( ActionWithValue* av );
48 : // active methods:
49 : void calculate() override;
50 : static void calculateCV( const ColvarInput<T>& cvin,
51 : ColvarOutput<T>& cvout );
52 : };
53 :
54 : template <typename T>
55 5296 : void Distance<T>::registerKeywords( Keywords& keys ) {
56 5296 : Colvar::registerKeywords( keys );
57 10592 : keys.setDisplayName("DISTANCE");
58 : constexpr auto scalarOrVector = Keywords::componentType::scalar
59 : | Keywords::componentType::vector;
60 5296 : keys.add("atoms","ATOMS","the pair of atom that we are calculating the distance between");
61 5296 : keys.addFlag("SQUARED",false,"calculate and store the square of the distance");
62 5296 : 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");
63 5296 : 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");
64 10592 : keys.addOutputComponent("x","COMPONENTS",scalarOrVector,"the x-component of the vector connecting the two atoms");
65 10592 : keys.addOutputComponent("y","COMPONENTS",scalarOrVector,"the y-component of the vector connecting the two atoms");
66 10592 : keys.addOutputComponent("z","COMPONENTS",scalarOrVector,"the z-component of the vector connecting the two atoms");
67 10592 : keys.addOutputComponent("a","SCALED_COMPONENTS",scalarOrVector,"the normalized projection on the first lattice vector of the vector connecting the two atoms");
68 10592 : keys.addOutputComponent("b","SCALED_COMPONENTS",scalarOrVector,"the normalized projection on the second lattice vector of the vector connecting the two atoms");
69 10592 : keys.addOutputComponent("c","SCALED_COMPONENTS",scalarOrVector,"the normalized projection on the third lattice vector of the vector connecting the two atoms");
70 5296 : keys.add("hidden","NO_ACTION_LOG","suppresses printing from action on the log");
71 5296 : keys.setValueDescription(scalarOrVector,"the DISTANCE between this pair of atoms");
72 5296 : keys.addDOI("10.1007/978-1-4939-9608-7_21");
73 10592 : keys.reset_style("NUMERICAL_DERIVATIVES","hidden");
74 5296 : }
75 :
76 : template <typename T>
77 1441 : Distance<T>::Distance(const ActionOptions&ao):
78 : PLUMED_COLVAR_INIT(ao),
79 1441 : pbc(true),
80 1441 : value(1) {
81 1441 : std::vector<AtomNumber> atoms;
82 1441 : parseAtomList(-1,atoms,this);
83 1441 : if(atoms.size()!=2) {
84 1 : error("Number of specified atoms should be 2");
85 : }
86 :
87 1440 : bool nopbc=!pbc;
88 1440 : parseFlag("NOPBC",nopbc);
89 1440 : pbc=!nopbc;
90 :
91 1440 : if(pbc) {
92 763 : log.printf(" using periodic boundary conditions\n");
93 : } else {
94 677 : log.printf(" without periodic boundary conditions\n");
95 : }
96 :
97 1440 : mode = static_cast<components>(getModeAndSetupValues(this));
98 1439 : if( mode == components::standard || mode == components::scaled ) {
99 275 : value.resize(3);
100 : }
101 1439 : requestAtoms(atoms);
102 1445 : }
103 :
104 : template <typename T>
105 231507 : void Distance<T>::parseAtomList( const int& num,
106 : std::vector<AtomNumber>& t,
107 : ActionAtomistic* aa ) {
108 463014 : aa->parseAtomList("ATOMS",num,t);
109 231507 : if( t.size()==2 ) {
110 231287 : aa->log.printf(" between atoms %d %d\n",t[0].serial(),t[1].serial());
111 : }
112 231507 : }
113 :
114 : template <typename T>
115 1659 : unsigned Distance<T>::getModeAndSetupValues( ActionWithValue* av ) {
116 : bool s;
117 1659 : av->parseFlag("SQUARED",s);
118 : bool c;
119 1659 : av->parseFlag("COMPONENTS",c);
120 : bool sc;
121 1659 : av->parseFlag("SCALED_COMPONENTS",sc);
122 1659 : if( c && sc ) {
123 1 : av->error("COMPONENTS and SCALED_COMPONENTS are not compatible");
124 : }
125 1658 : if ( s & c ) {
126 0 : av->error("SQUARED and COMPONENTS are not compatible");
127 : }
128 1658 : if ( s & sc ) {
129 0 : av->error("SQUARED and SCALED_COMPONENTS are not compatible");
130 : }
131 1658 : if(c) {
132 678 : av->addComponentWithDerivatives("x");
133 339 : av->componentIsNotPeriodic("x");
134 678 : av->addComponentWithDerivatives("y");
135 339 : av->componentIsNotPeriodic("y");
136 678 : av->addComponentWithDerivatives("z");
137 339 : av->componentIsNotPeriodic("z");
138 339 : av->log<<" WARNING: components will not have the proper periodicity - see manual\n";
139 339 : return components::standard;
140 1319 : } else if(sc) {
141 12 : av->addComponentWithDerivatives("a");
142 12 : av->componentIsPeriodic("a","-0.5","+0.5");
143 12 : av->addComponentWithDerivatives("b");
144 12 : av->componentIsPeriodic("b","-0.5","+0.5");
145 12 : av->addComponentWithDerivatives("c");
146 12 : av->componentIsPeriodic("c","-0.5","+0.5");
147 6 : return components::scaled;
148 : }
149 1313 : av->addValueWithDerivatives();
150 1313 : av->setNotPeriodic();
151 1313 : if(s) {
152 1 : return components::squared;
153 : }
154 : return components::scalar;
155 : }
156 :
157 : // calculator
158 : template <typename T>
159 218599 : void Distance<T>::calculate() {
160 :
161 218599 : if(pbc) {
162 83334 : makeWhole();
163 : }
164 : using CVInput= ColvarInput<double>;
165 : using CVOutput=ColvarOutput<double>;
166 218599 : auto cvout = CVOutput::createColvarOutput(value, derivs, this);
167 218599 : Distance<double>::calculateCV( CVInput::createColvarInput( mode, getPositions(), this ), cvout );
168 218599 : if( mode==components::standard ) {
169 2001 : Value* valuex=getPntrToComponent("x");
170 2001 : Value* valuey=getPntrToComponent("y");
171 2001 : Value* valuez=getPntrToComponent("z");
172 :
173 6003 : for(unsigned i=0; i<2; ++i) {
174 4002 : setAtomsDerivatives(valuex,i,cvout.getAtomDerivatives(0,i) );
175 : }
176 2001 : setBoxDerivatives(valuex,cvout.virial[0]);
177 2001 : valuex->set(value[0]);
178 :
179 6003 : for(unsigned i=0; i<2; ++i) {
180 4002 : setAtomsDerivatives(valuey,i,cvout.getAtomDerivatives(1,i) );
181 : }
182 2001 : setBoxDerivatives(valuey,cvout.virial[1]);
183 2001 : valuey->set(value[1]);
184 :
185 6003 : for(unsigned i=0; i<2; ++i) {
186 4002 : setAtomsDerivatives(valuez,i,cvout.getAtomDerivatives(2,i) );
187 : }
188 2001 : setBoxDerivatives(valuez,cvout.virial[2]);
189 2001 : valuez->set(value[2]);
190 216598 : } else if( mode==components::scaled ) {
191 :
192 26 : Value* valuea=getPntrToComponent("a");
193 26 : Value* valueb=getPntrToComponent("b");
194 26 : Value* valuec=getPntrToComponent("c");
195 78 : for(unsigned i=0; i<2; ++i) {
196 52 : setAtomsDerivatives(valuea,i,cvout.getAtomDerivatives(0,i) );
197 : }
198 26 : valuea->set(value[0]);
199 78 : for(unsigned i=0; i<2; ++i) {
200 52 : setAtomsDerivatives(valueb,i,cvout.getAtomDerivatives(1,i) );
201 : }
202 26 : valueb->set(value[1]);
203 78 : for(unsigned i=0; i<2; ++i) {
204 52 : setAtomsDerivatives(valuec,i,cvout.getAtomDerivatives(2,i) );
205 : }
206 26 : valuec->set(value[2]);
207 : } else {
208 649716 : for(unsigned i=0; i<2; ++i) {
209 433144 : setAtomsDerivatives(i,cvout.getAtomDerivatives(0,i) );
210 : }
211 216572 : setBoxDerivatives(cvout.virial[0]);
212 216572 : setValue (value[0]);
213 : }
214 218599 : }
215 :
216 : template <typename T>
217 2739991 : void Distance<T>::calculateCV( const ColvarInput<T>& cvin, ColvarOutput<T>& cvout ) {
218 2739991 : auto distance=delta(cvin.pos[0],cvin.pos[1]);
219 : //these six can be compacted in a scoped `using namespac3 PLMD::Versors` but plumedcheck won't approve (even in a comment, hence the 3)
220 : using PLMD::Versors::xp;
221 : using PLMD::Versors::xm;
222 : using PLMD::Versors::yp;
223 : using PLMD::Versors::ym;
224 : using PLMD::Versors::zp;
225 : using PLMD::Versors::zm;
226 2739991 : if(cvin.mode==components::standard) {
227 :
228 : cvout.derivs[0][0] = xm<T>;
229 : cvout.derivs[0][1] = xp<T>;
230 106175 : cvout.values[0] = distance[0];
231 :
232 : cvout.derivs[1][0] = ym<T>;
233 : cvout.derivs[1][1] = yp<T>;
234 106175 : cvout.values[1] = distance[1];
235 :
236 : cvout.derivs[2][0] = zm<T>;
237 : cvout.derivs[2][1] = zp<T>;
238 106175 : cvout.values[2] = distance[2];
239 106175 : ColvarInput<T>::setBoxDerivativesNoPbc( cvin, cvout );
240 2633816 : } else if(cvin.mode==components::scaled) {
241 41 : auto d=cvin.pbc.realToScaled(distance);
242 41 : cvout.derivs[0][0].copyConv(matmul(cvin.pbc.getInvBox(),xm<T>));
243 41 : cvout.derivs[0][1].copyConv(matmul(cvin.pbc.getInvBox(),xp<T>));
244 41 : cvout.values[0] = Tools::pbc(d[0]);
245 41 : cvout.derivs[1][0].copyConv(matmul(cvin.pbc.getInvBox(),ym<T>));
246 41 : cvout.derivs[1][1].copyConv(matmul(cvin.pbc.getInvBox(),yp<T>));
247 41 : cvout.values[1] = Tools::pbc(d[1]);
248 41 : cvout.derivs[2][0].copyConv(matmul(cvin.pbc.getInvBox(),zm<T>));
249 41 : cvout.derivs[2][1].copyConv(matmul(cvin.pbc.getInvBox(),zp<T>));
250 41 : cvout.values[2] = Tools::pbc(d[2]);
251 2633775 : } else if(cvin.mode==components::squared) {
252 : const double value=distance.modulo2();
253 : cvout.derivs[0][0] = -2*distance;
254 : cvout.derivs[0][1] = 2*distance;
255 13400 : ColvarInput<T>::setBoxDerivativesNoPbc( cvin, cvout );
256 13400 : cvout.values[0] = value;
257 : } else {
258 2620375 : const double value=distance.modulo();
259 2620375 : const double invvalue=1.0/value;
260 2620375 : cvout.derivs[0][0] = -invvalue*distance;
261 : cvout.derivs[0][1] = invvalue*distance;
262 2620375 : ColvarInput<T>::setBoxDerivativesNoPbc( cvin, cvout );
263 2620375 : cvout.values[0] = value;
264 : }
265 2739991 : }
266 :
267 : } //namespace colvar
268 : } //namespace PLMD
269 :
270 : #endif //__PLUMED_colvar_Distance_h
271 :
272 :
|
