Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) crystdistrib 2023-2023 The code team
3 : (see the PEOPLE-crystdistrib file at the root of this folder for a list of names)
4 :
5 : This file is part of crystdistrib code module.
6 :
7 : The crystdistrib code module is free software: you can redistribute it and/or modify
8 : it under the terms of the GNU Lesser General Public License as published by
9 : the Free Software Foundation, either version 3 of the License, or
10 : (at your option) any later version.
11 :
12 : The crystdistrib code module is distributed in the hope that it will be useful,
13 : but WITHOUT ANY WARRANTY; without even the implied warranty of
14 : MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
15 : GNU Lesser General Public License for more details.
16 :
17 : You should have received a copy of the GNU Lesser General Public License
18 : along with the crystdistrib code module. If not, see <http://www.gnu.org/licenses/>.
19 : +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
20 : #include "core/Colvar.h"
21 : #include "colvar/ColvarShortcut.h"
22 : #include "colvar/MultiColvarTemplate.h"
23 : #include "core/ActionRegister.h"
24 : #include "tools/Pbc.h"
25 :
26 : namespace PLMD {
27 : namespace crystdistrib {
28 :
29 : //+PLUMEDOC COLVAR QUATERNION
30 : /*
31 : Calculate unit quaternions for molecules.
32 :
33 : This action calculates a unit quaternion to define an internal coordinate frame for a molecule. The reference frame for the molecule is user-defined using the positions of three (non-collinear) atoms.
34 : The vectors which will define the frame are calculated as follows from atomic coordinates, all vectors in $\mathbb{R}^3$, $x_1, x_2, x_3$:
35 :
36 : The first axis is the normalized difference of $x_1$ and $x_2$:
37 :
38 : $$
39 : \mathbf{\hat{v_1}} = x_2 - x_1 / \|x_2 - x_1\|
40 : $$
41 :
42 : In general, the vector $\mathbf{v'_2} = x_3 - x_1$ will not be orthogonal to $\mathbf{\hat{v_1}}$. This is fixed by taking the difference between the projection of $\mathbf{v'_2}$
43 : onto $\mathbf{\hat{v_1}}$ and $\mathbf{\hat{v_1}}$.
44 :
45 : $$
46 : \mathbf{v_2} = \mathbf{v'_2} - Proj_{\mathbf{\hat{v}_1}}(\mathbf{v'_2}) = \mathbf{v'_2} - \mathbf{\hat{v}_1} \cdot \mathbf{v'_2}
47 : $$
48 :
49 : This is then normalized to form the second axis.
50 :
51 : $$
52 : \mathbf{\hat{v_2}} = \mathbf{v_2} / \|\mathbf{v_2}\|
53 : $$
54 :
55 : Finally, the third axis is the cross product between the first two.
56 :
57 : $$
58 : \mathbf{\hat{v_3}} = \mathbf{\hat{v_1}} \times \mathbf{\hat{v_2}}
59 : $$
60 :
61 : The above 9 components form an orthonormal basis, centered on the molecule provided. The rotation matrix is generally the inverse of this matrix, and
62 : in this case since the matrix is orthogonal and its determinant is 1, the inverse is simply the transpose. The rotation matrix is then converted to a quaternion.
63 : The resulting quaternion has 4 real numbers attached to it, and they can be called as w, i , j ,and k. Note the quaternions are not unique e.g. q and -q perform the same rotation,
64 : so take care when using the results. Take care that the components are simply 4 real numbers, and the usual non-commutativity of quaternions, and any other algebraic difference
65 : will need to be accounted for manually in later usage. No checks are made for co-linearity, or if the atoms are a part of the same molecule.
66 :
67 : An example input file, in a system of 12 atoms. It calculates four molecular frames, then uses them in an order parameter.
68 :
69 : ```plumed
70 : q1: QUATERNION ATOMS1=1,3,2 ATOMS2=4,6,5
71 : q2: QUATERNION ATOMS1=7,9,8 ATOMS2=10,12,11
72 : #There are no checks to make sure the atoms belong to the same molecule
73 :
74 : fake: CUSTOM ARG=q1.w,q1.i,q1.j,q1.k,q2.w,q2.i,q2.j,q2.k VAR=w1,i1,j1,k1,w2,i2,j2,k2 FUNC=w1*w2+i1*i2+j1*j2+k1*k2 PERIODIC=NO
75 : #this isn’t a real order parameter, for the record
76 : PRINT ARG=fake FILE=fakeout
77 : ```
78 :
79 : The last thing to note is that by default a procedure akin to that used in [WHOLEMOLECULES](WHOLEMOLECULES.md)
80 : is used to ensure that the sets of atoms that are specified to each ATOMS keyword are not broken by the periodic
81 : boundary conditions. If you would like to turn this off for any reason you add the NOPBC in your input file as shown
82 : below:
83 :
84 : ```plumed
85 : d: QUATERNION ATOMS=1,2,3 NOPBC
86 : ```
87 :
88 : */
89 : //+ENDPLUMEDOC
90 :
91 : //simple hamilton/quaternion product, which is just expanding the two quats as expected, then applying the rules for i j k
92 : //passing 12 references might be a bit silly
93 : //void QuatProduct(double &a1, double &b1, double &c1, double &d1, double &a2, double &b2, double &c2, double &d2, double &w, double &i, double &j, double &k) {
94 : //
95 : // w = a1*a2 - b1*b2 - c1*c2 - d1*d2;
96 : // i = a1*b2 + b1*a2 + c1*d2 - d1*c2;
97 : // j = a1*c2 - b1*d2 + c1*a2 + d1*b2;
98 : // k = a1*d2 + b1*c2 - c1*b2 + d1*a2;
99 : //}
100 : //
101 : //
102 : //
103 :
104 : class Quaternion : public Colvar {
105 : private:
106 : bool pbc;
107 : std::vector<double> value;
108 : std::vector<double> derivs;
109 : public:
110 : static void registerKeywords( Keywords& keys );
111 : explicit Quaternion(const ActionOptions&);
112 : static void parseAtomList( const int& num, std::vector<AtomNumber>& t, ActionAtomistic* aa );
113 : static unsigned getModeAndSetupValues( ActionWithValue* av );
114 : // active methods:
115 : void calculate() override;
116 : static void calculateCV( const colvar::ColvarInput& cvin, ColvarOutput& cvout );
117 : };
118 :
119 : typedef colvar::ColvarShortcut<Quaternion> QuaternionShortcut;
120 : PLUMED_REGISTER_ACTION(QuaternionShortcut,"QUATERNION")
121 : PLUMED_REGISTER_ACTION(Quaternion,"QUATERNION_SCALAR")
122 : typedef colvar::MultiColvarTemplate<Quaternion> QuaternionMulti;
123 : PLUMED_REGISTER_ACTION(QuaternionMulti,"QUATERNION_VECTOR")
124 :
125 90 : void Quaternion::registerKeywords( Keywords& keys ) {
126 90 : Colvar::registerKeywords( keys );
127 180 : keys.setDisplayName("QUATERNION");
128 90 : keys.add("atoms","ATOMS","the three atom that we are using to calculate the quaternion");
129 180 : keys.addOutputComponent("w","default","scalar/vector","the real component of quaternion");
130 180 : keys.addOutputComponent("i","default","scalar/vector","the i component of the quaternion");
131 180 : keys.addOutputComponent("j","default","scalar/vector","the j component of the quaternion");
132 180 : keys.addOutputComponent("k","default","scalar/vector","the k component of the quaternion");
133 90 : keys.add("hidden","NO_ACTION_LOG","suppresses printing from action on the log");
134 180 : keys.reset_style("NUMERICAL_DERIVATIVES","hidden");
135 90 : }
136 :
137 5 : Quaternion::Quaternion(const ActionOptions&ao):
138 : PLUMED_COLVAR_INIT(ao),
139 5 : pbc(true),
140 5 : value(4) {
141 : std::vector<AtomNumber> atoms;
142 5 : parseAtomList(-1,atoms,this);
143 5 : if(atoms.size()!=3) {
144 0 : error("Number of specified atoms should be 3");
145 : }
146 : // please note, I do NO checking if these atoms are in the same molecule at all, so be careful in your input
147 :
148 5 : bool nopbc=!pbc;
149 5 : parseFlag("NOPBC",nopbc);
150 5 : pbc=!nopbc;
151 :
152 5 : unsigned mode = getModeAndSetupValues( this );
153 5 : requestAtoms(atoms);
154 5 : }
155 :
156 1322 : void Quaternion::parseAtomList( const int& num, std::vector<AtomNumber>& t, ActionAtomistic* aa ) {
157 2644 : aa->parseAtomList("ATOMS",num,t);
158 1322 : if( t.size()==3 ) {
159 1309 : aa->log.printf(" involving atoms %d %d %d\n",t[0].serial(),t[1].serial(),t[0].serial());
160 : }
161 1322 : }
162 :
163 18 : unsigned Quaternion::getModeAndSetupValues( ActionWithValue* av ) {
164 : // This sets up values that we can pass around in PLUMED
165 36 : av->addComponentWithDerivatives("w");
166 18 : av->componentIsNotPeriodic("w");
167 36 : av->addComponentWithDerivatives("i");
168 18 : av->componentIsNotPeriodic("i");
169 36 : av->addComponentWithDerivatives("j");
170 18 : av->componentIsNotPeriodic("j");
171 36 : av->addComponentWithDerivatives("k");
172 18 : av->componentIsNotPeriodic("k");
173 18 : return 0;
174 : }
175 :
176 45 : void Quaternion::calculate() {
177 45 : if(pbc) {
178 45 : makeWhole();
179 : }
180 :
181 45 : ColvarOutput cvout = ColvarOutput::createColvarOutput(value,derivs,this);
182 45 : calculateCV( colvar::ColvarInput::createColvarInput( 0, getPositions(), this ), cvout );
183 225 : for(unsigned j=0; j<4; ++j) {
184 180 : Value* valuej=getPntrToComponent(j);
185 720 : for(unsigned i=0; i<3; ++i) {
186 540 : setAtomsDerivatives(valuej,i,cvout.getAtomDerivatives(j,i) );
187 : }
188 180 : setBoxDerivatives(valuej,cvout.virial[j]);
189 180 : valuej->set(value[j]);
190 : }
191 45 : }
192 :
193 : // calculator
194 7203 : void Quaternion::calculateCV( const colvar::ColvarInput& cvin, ColvarOutput& cvout ) {
195 : //declarations
196 7203 : Vector vec1_comp = delta( cvin.pos[0], cvin.pos[1] ); //components between atom 1 and 2
197 7203 : Vector vec2_comp = delta( cvin.pos[0], cvin.pos[2] ); //components between atom 1 and 3
198 :
199 : ////////x-vector calculations///////
200 7203 : double magx = vec1_comp.modulo();
201 : Vector xt = vec1_comp / magx;
202 7203 : std::vector<Tensor> dx(3);
203 7203 : double magx3= magx*magx*magx;
204 : //dx[i] - derivatives of atom i's coordinates
205 7203 : dx[0](0,0) = ( -(vec1_comp[1]*vec1_comp[1]+vec1_comp[2]*vec1_comp[2])/magx3 ); //dx[0]/dx0
206 7203 : dx[0](0,1) = ( vec1_comp[0]*vec1_comp[1]/magx3 ); // dx[0]/dy0
207 7203 : dx[0](0,2) = ( vec1_comp[0]*vec1_comp[2]/magx3 ); // dx[0]/dz0
208 7203 : dx[0](1,0) = ( vec1_comp[1]*vec1_comp[0]/magx3 ); // dx[1]/dx0
209 7203 : dx[0](1,1) = ( -(vec1_comp[0]*vec1_comp[0]+vec1_comp[2]*vec1_comp[2])/magx3 ); // dx[1]/dy0
210 7203 : dx[0](1,2) = ( vec1_comp[1]*vec1_comp[2]/magx3 ); //dx[1]/dz0
211 7203 : dx[0](2,0) = ( vec1_comp[2]*vec1_comp[0]/magx3 );//etc
212 7203 : dx[0](2,1) = ( vec1_comp[2]*vec1_comp[1]/magx3 );
213 7203 : dx[0](2,2) = ( -(vec1_comp[1]*vec1_comp[1]+vec1_comp[0]*vec1_comp[0])/magx3 );
214 :
215 7203 : dx[1](0,0) = ( (vec1_comp[1]*vec1_comp[1]+vec1_comp[2]*vec1_comp[2])/magx3 );//dx[0]/dx1
216 7203 : dx[1](0,1) = ( -vec1_comp[0]*vec1_comp[1]/magx3 );//dx[0]/dy1
217 7203 : dx[1](0,2) = ( -vec1_comp[0]*vec1_comp[2]/magx3 );
218 7203 : dx[1](1,0) = ( -vec1_comp[1]*vec1_comp[0]/magx3 );
219 7203 : dx[1](1,1) = ( (vec1_comp[0]*vec1_comp[0]+vec1_comp[2]*vec1_comp[2])/magx3 );
220 7203 : dx[1](1,2) = ( -vec1_comp[1]*vec1_comp[2]/magx3 );
221 7203 : dx[1](2,0) = ( -vec1_comp[2]*vec1_comp[0]/magx3 );
222 7203 : dx[1](2,1) = ( -vec1_comp[2]*vec1_comp[1]/magx3 );
223 7203 : dx[1](2,2) = ( (vec1_comp[1]*vec1_comp[1]+vec1_comp[0]*vec1_comp[0])/magx3 );
224 : dx[2].zero();//not atom[2] terms present
225 :
226 : ////////y-vector calculations////////
227 : //project vec2_comp on to vec1_comp
228 : //first do dot product of unormalized x and unormed y, divided by magnitude of x^2
229 : double dp = dotProduct( vec1_comp, vec2_comp );
230 : double magx2=magx*magx;
231 7203 : std::vector<Vector> fac_derivs(3);
232 7203 : double magx4=magx2*magx2, fac = dp/magx2; //fac meaning factor on front
233 7203 : fac_derivs[0] = (-vec2_comp - vec1_comp)/magx2 + 2*dp*vec1_comp / magx4;
234 7203 : fac_derivs[1] = (vec2_comp)/(magx2) - 2*dp*vec1_comp / magx4;
235 7203 : fac_derivs[2] = (vec1_comp)/(magx2); //atom 1, components x2,y2,z2
236 : //now multiply fac by unormed x, and subtract it from unormed y, then normalize
237 : Vector yt = vec2_comp - fac*vec1_comp;
238 7203 : std::vector<Tensor> dy(3);
239 7203 : dy[0](0,0) = -1 - fac_derivs[0][0]*vec1_comp[0] + fac; // dy[0]/dx0
240 7203 : dy[0](0,1) = -fac_derivs[0][1]*vec1_comp[0]; // dy[0]/dy0
241 7203 : dy[0](0,2) = -fac_derivs[0][2]*vec1_comp[0];
242 7203 : dy[0](1,0) = -fac_derivs[0][0]*vec1_comp[1];
243 7203 : dy[0](1,1) = -1 - fac_derivs[0][1]*vec1_comp[1] + fac;
244 7203 : dy[0](1,2) = -fac_derivs[0][2]*vec1_comp[1];
245 7203 : dy[0](2,0) = -fac_derivs[0][0]*vec1_comp[2];
246 7203 : dy[0](2,1) = -fac_derivs[0][1]*vec1_comp[2];
247 7203 : dy[0](2,2) = -1 - fac_derivs[0][2]*vec1_comp[2] + fac;
248 :
249 7203 : dy[1](0,0) = -fac_derivs[1][0]*vec1_comp[0] - fac; //dy[0]/dx0
250 7203 : dy[1](0,1) = -fac_derivs[1][1]*vec1_comp[0];
251 7203 : dy[1](0,2) = -fac_derivs[1][2]*vec1_comp[0];
252 7203 : dy[1](1,0) = -fac_derivs[1][0]*vec1_comp[1];
253 7203 : dy[1](1,1) = -fac_derivs[1][1]*vec1_comp[1] - fac;
254 7203 : dy[1](1,2) = -fac_derivs[1][2]*vec1_comp[1];
255 7203 : dy[1](2,0) = -fac_derivs[1][0]*vec1_comp[2];
256 7203 : dy[1](2,1) = -fac_derivs[1][1]*vec1_comp[2];
257 7203 : dy[1](2,2) = -fac_derivs[1][2]*vec1_comp[2] - fac;
258 :
259 7203 : dy[2](0,0) = 1 - fac_derivs[2][0]*vec1_comp[0];//dy[0]/dx2
260 7203 : dy[2](0,1) = -fac_derivs[2][1]*vec1_comp[0];
261 7203 : dy[2](0,2) = -fac_derivs[2][2]*vec1_comp[0];
262 7203 : dy[2](1,0) = -fac_derivs[2][0]*vec1_comp[1];
263 7203 : dy[2](1,1) = 1 - fac_derivs[2][1]*vec1_comp[1];
264 7203 : dy[2](1,2) = -fac_derivs[2][2]*vec1_comp[1];
265 7203 : dy[2](2,0) = -fac_derivs[2][0]*vec1_comp[2];
266 7203 : dy[2](2,1) = -fac_derivs[2][1]*vec1_comp[2];
267 7203 : dy[2](2,2) = 1 - fac_derivs[2][2]*vec1_comp[2];
268 : //now normalize, and we have our y vector
269 7203 : double magy = yt.modulo();
270 7203 : double imagy = 1/magy, magy3 = magy*magy*magy;
271 : Tensor abc;
272 28812 : for(unsigned i=0; i<3; ++i) {
273 : abc.setRow(i, yt);
274 : }
275 : Tensor abc_diag;
276 7203 : abc_diag.setRow(0, Vector(yt[0], 0, 0));
277 7203 : abc_diag.setRow(1, Vector(0, yt[1], 0));
278 7203 : abc_diag.setRow(2, Vector(0, 0, yt[2]));
279 7203 : Tensor abc_prod = matmul(abc_diag, abc);
280 28812 : for(unsigned i=0; i<3; ++i) {
281 43218 : dy[i] = dy[i]/magy - matmul(abc_prod, dy[i])/magy3;
282 : }
283 : //normalize now, derivatives are with respect to un-normalized y vector
284 7203 : yt = yt / magy;
285 :
286 : ///////z-vector calculations/////////
287 : //comparatively simple
288 7203 : Vector zt = crossProduct(xt,yt);
289 7203 : std::vector<Tensor> dz(3);
290 21609 : dz[0].setCol( 0, crossProduct( dx[0].getCol(0), yt ) + crossProduct( xt, dy[0].getCol(0) ) );
291 21609 : dz[0].setCol( 1, crossProduct( dx[0].getCol(1), yt ) + crossProduct( xt, dy[0].getCol(1) ) );
292 21609 : dz[0].setCol( 2, crossProduct( dx[0].getCol(2), yt ) + crossProduct( xt, dy[0].getCol(2) ) );
293 :
294 21609 : dz[1].setCol( 0, crossProduct( dx[1].getCol(0), yt ) + crossProduct( xt, dy[1].getCol(0) ) );
295 21609 : dz[1].setCol( 1, crossProduct( dx[1].getCol(1), yt ) + crossProduct( xt, dy[1].getCol(1) ) );
296 21609 : dz[1].setCol( 2, crossProduct( dx[1].getCol(2), yt ) + crossProduct( xt, dy[1].getCol(2) ) );
297 :
298 14406 : dz[2].setCol( 0, crossProduct( xt, dy[2].getCol(0) ) );
299 14406 : dz[2].setCol( 1, crossProduct( xt, dy[2].getCol(1) ) );
300 7203 : dz[2].setCol( 2, crossProduct( xt, dy[2].getCol(2) ) );
301 :
302 : //for debugging frame values
303 : //aa->log.printf("%8.6f %8.6f %8.6f\n%8.6f %8.6f %8.6f\n%8.6f %8.6f %8.6f\n",xt[0],xt[1],xt[2],yt[0],yt[1],yt[2],zt[0],zt[1],zt[2]);
304 :
305 : //for bebuffing derivatives
306 : //aa->log.printf("x1 x2 x3 y1 y2 y3 z1 z2 z3\n");
307 : //for (int i=0; i<3; i++){
308 : //for (int j=0;j<3;j++){
309 : //aa->log.printf("%8.4f %8.4f %8.4f\n%8.4f %8.4f %8.4f\n%8.4f %8.4f %8.4f\n",dx[i](0,j), dx[i](1,j), dx[i](2,j), dy[i](0,j), dy[i](1,j), dy[i](2,j), dz[i](0,j), dz[i](1,j), dz[i](2,j));
310 : //}
311 : //}
312 : //
313 :
314 : //the above 9 components form an orthonormal basis, centered on the molecule in question
315 : //the rotation matrix is generally the inverse of this matrix, and in this case since it is 1) orthogonal and 2) its determinant is 1
316 : //the inverse is simply the transpose
317 :
318 :
319 : //[x[0] x[1] x[2]]
320 : //[y[0] y[1] y[2]]
321 : //[z[0] z[1] z[2]]
322 : //QUICKFIX to transpose basis
323 7203 : Vector x(xt[0],yt[0],zt[0]);
324 7203 : Vector y(xt[1],yt[1],zt[1]);
325 7203 : Vector z(xt[2],yt[2],zt[2]);
326 :
327 : //likewise transposing the tensors into proper form
328 7203 : std::vector<Tensor> tdx(3);
329 7203 : std::vector<Tensor> tdy(3);
330 7203 : std::vector<Tensor> tdz(3);
331 28812 : for (int i=0; i<3; ++i) {
332 43218 : tdx[i].setRow(0, dx[i].getRow(0));
333 21609 : tdx[i].setRow(1, dy[i].getRow(0));
334 21609 : tdx[i].setRow(2, dz[i].getRow(0));
335 :
336 21609 : tdy[i].setRow(0, dx[i].getRow(1));
337 21609 : tdy[i].setRow(1, dy[i].getRow(1));
338 21609 : tdy[i].setRow(2, dz[i].getRow(1));
339 :
340 21609 : tdz[i].setRow(0, dx[i].getRow(2));
341 21609 : tdz[i].setRow(1, dy[i].getRow(2));
342 21609 : tdz[i].setRow(2, dz[i].getRow(2));
343 : }
344 :
345 : //convert to quaternion
346 7203 : double tr = x[0] + y[1] + z[2] + 1; //trace of the rotation matrix + 1
347 7203 : std::vector<Vector> dS(3);
348 7203 : if (tr > 1.0E-8) { //to avoid numerical instability
349 7203 : double S = 1/(sqrt(tr) * 2); // S=4*qw
350 28812 : for(unsigned i=0; i<3; ++i) {
351 43218 : dS[i] = (-2*S*S*S)*(tdx[i].getRow(0) + tdy[i].getRow(1) + tdz[i].getRow(2));
352 : }
353 :
354 7203 : cvout.values[0] = 0.25 / S;
355 28812 : for(unsigned i=0; i<3; ++i) {
356 21609 : cvout.derivs[0][i] =-0.25*dS[i]/(S*S);
357 : }
358 :
359 7203 : cvout.values[1] = (z[1] - y[2]) * S;
360 28812 : for(unsigned i=0; i<3; ++i) {
361 43218 : cvout.derivs[1][i] = (S)*(tdz[i].getRow(1) - tdy[i].getRow(2)) + (z[1]-y[2])*dS[i];
362 : }
363 :
364 7203 : cvout.values[2] = (x[2] - z[0]) * S;
365 28812 : for(unsigned i=0; i<3; ++i) {
366 43218 : cvout.derivs[2][i] = (S)*(tdx[i].getRow(2) - tdz[i].getRow(0)) + (x[2]-z[0])*dS[i];
367 : }
368 :
369 7203 : cvout.values[3] = (y[0] - x[1]) * S;
370 28812 : for(unsigned i=0; i<3; ++i) {
371 43218 : cvout.derivs[3][i] = (S)*(tdy[i].getRow(0) - tdx[i].getRow(1)) + (y[0]-x[1])*dS[i];
372 : }
373 0 : } else if ((x[0] > y[1])&(x[0] > z[2])) {
374 0 : double S = sqrt(1.0 + x[0] - y[1] - z[2]) * 2; // S=4*qx
375 0 : for(unsigned i=0; i<3; ++i) {
376 0 : dS[i] = (2/S)*(tdx[i].getRow(0) - tdy[i].getRow(1) - tdz[i].getRow(2));
377 : }
378 :
379 0 : cvout.values[0] = (z[1] - y[2]) / S;
380 0 : for(unsigned i=0; i<3; ++i) {
381 0 : cvout.derivs[0][i] = (1/S)*(tdz[i].getRow(1) - tdy[i].getRow(2)) - (cvout.values[0]/S)*dS[i];
382 : }
383 :
384 0 : cvout.values[1] = 0.25 * S;
385 0 : for(unsigned i=0; i<3; ++i) {
386 0 : cvout.derivs[1][i] =0.25*dS[i];
387 : }
388 :
389 0 : cvout.values[2] = (x[1] + y[0]) / S;
390 0 : for(unsigned i=0; i<3; ++i) {
391 0 : cvout.derivs[2][i] = (1/S)*(tdx[i].getRow(1) + tdy[i].getRow(0)) - (cvout.values[2]/S)*dS[i];
392 : }
393 :
394 0 : cvout.values[3] = (x[2] + z[0]) / S;
395 0 : for(unsigned i=0; i<3; ++i) {
396 0 : cvout.derivs[3][i] = (1/S)*(tdx[i].getRow(2) + tdz[i].getRow(0)) - (cvout.values[3]/S)*dS[i];
397 : }
398 0 : } else if (y[1] > z[2]) {
399 0 : double S = sqrt(1.0 + y[1] - x[0] - z[2]) * 2; // S=4*qy
400 0 : for(unsigned i=0; i<3; ++i) {
401 0 : dS[i] = (2/S)*( -tdx[i].getRow(0) + tdy[i].getRow(1) - tdz[i].getRow(2));
402 : }
403 :
404 :
405 0 : cvout.values[0] = (x[2] - z[0]) / S;
406 0 : for(unsigned i=0; i<3; ++i) {
407 0 : cvout.derivs[0][i] = (1/S)*(tdx[i].getRow(2) - tdz[i].getRow(0)) - (cvout.values[0]/S)*dS[i];
408 : }
409 :
410 0 : cvout.values[1] = (x[1] + y[0]) / S;
411 0 : for(unsigned i=0; i<3; ++i) {
412 0 : cvout.derivs[1][i] = (1/S)*(tdx[i].getRow(1) + tdy[i].getRow(0)) - (cvout.values[1]/S)*dS[i];
413 : }
414 :
415 0 : cvout.values[2] = 0.25 * S;
416 0 : for(unsigned i=0; i<3; ++i) {
417 0 : cvout.derivs[2][i] =0.25*dS[i];
418 : }
419 :
420 0 : cvout.values[3] = (y[2] + z[1]) / S;
421 0 : for(unsigned i=0; i<3; ++i) {
422 0 : cvout.derivs[3][i] = (1/S)*(tdy[i].getRow(2) + tdz[i].getRow(1)) - (cvout.values[3]/S)*dS[i];
423 : }
424 : } else {
425 0 : double S = sqrt(1.0 + z[2] - x[0] - y[1]) * 2; // S=4*qz
426 0 : for(unsigned i=0; i<3; ++i) {
427 0 : dS[i] = (2/S)*(-tdx[i].getRow(0) - tdy[i].getRow(1) + tdz[i].getRow(2));
428 : }
429 :
430 :
431 0 : cvout.values[0] = (y[0] - x[1]) / S;
432 0 : for(unsigned i=0; i<3; ++i) {
433 0 : cvout.derivs[0][i] = (1/S)*(tdy[i].getRow(0) - tdx[i].getRow(1)) - (cvout.values[0]/S)*dS[i];
434 : }
435 :
436 0 : cvout.values[1] = (x[2] + z[0]) / S;
437 0 : for(unsigned i=0; i<3; ++i) {
438 0 : cvout.derivs[1][i] = (1/S)*(tdx[i].getRow(2) + tdz[i].getRow(0)) - (cvout.values[1]/S)*dS[i];
439 : }
440 :
441 0 : cvout.values[2] = (y[2] + z[1]) / S;
442 0 : for(unsigned i=0; i<3; ++i) {
443 0 : cvout.derivs[2][i] = (1/S)*(tdy[i].getRow(2) + tdz[i].getRow(1)) - (cvout.values[2]/S)*dS[i];
444 : }
445 :
446 0 : cvout.values[3] = 0.25 * S;
447 0 : for(unsigned i=0; i<3; ++i) {
448 0 : cvout.derivs[3][i] =0.25*dS[i];
449 : }
450 : }
451 7203 : colvar::ColvarInput::setBoxDerivativesNoPbc( cvin, cvout );
452 :
453 7203 : }
454 :
455 : }
456 : }
457 :
458 :
459 :
|
