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 quaternions for molecules.
32 :
33 : The reference frame for the molecule is defined using the positions of three user selected atoms. From the positions of these atoms,
34 : \f$\mathbf{x}_1\f$, \f$\mathbf{x}_2\f$ and \f$\mathbf{x}_3\f$, we define the vectors of the reference frame as:
35 :
36 : \f[
37 : \begin{aligned}
38 : \mathbf{x} & = \mathbf{x}_2 - \mathbf{x}_1 \\
39 : \mathbf{y} & = (\mathbf{x}_2 - \mathbf{x}_1) \times (\mathbf{x}_3 - \mathbf{x}_1) \\
40 : \mathbf{z} & \mathbf{x} \times \mathbf{y}
41 : \f]
42 :
43 : \par Examples
44 :
45 : This calculates the quaternions for a molecule with 10 atoms
46 :
47 : \plumedfile
48 : q1: QUATERNION ATOMS1=1,2,3
49 : PRINT ARG=q1.w,q1.i,q1.j,q1.k FILE=colvar
50 : \endplumedfile
51 :
52 : This calculate the quaternions for two molecules with 10 atoms
53 :
54 : \plumedfile
55 : q1: QUATERNION ATOMS1=1,2,3 ATOMS=4,5,6
56 : PRINT ARG=q1.w,q1.i,q1.j,q1.k FILE=colvar
57 : \endplumedfile
58 :
59 : */
60 : //+ENDPLUMEDOC
61 :
62 : //+PLUMEDOC COLVAR QUATERNION_SCALAR
63 : /*
64 : Calculate a single quaternion
65 :
66 : See \ref QUATERNION for more details
67 :
68 : \par Examples
69 :
70 : */
71 : //+ENDPLUMEDOC
72 :
73 : //+PLUMEDOC COLVAR QUATERNION_VECTOR
74 : /*
75 : Calculate multiple quaternions
76 :
77 : See \ref QUATERNION for more details
78 :
79 : \par Examples
80 :
81 : */
82 : //+ENDPLUMEDOC
83 :
84 : //simple hamilton/quaternion product, which is just expanding the two quats as expected, then applying the rules for i j k
85 : //passing 12 references might be a bit silly
86 : //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) {
87 : //
88 : // w = a1*a2 - b1*b2 - c1*c2 - d1*d2;
89 : // i = a1*b2 + b1*a2 + c1*d2 - d1*c2;
90 : // j = a1*c2 - b1*d2 + c1*a2 + d1*b2;
91 : // k = a1*d2 + b1*c2 - c1*b2 + d1*a2;
92 : //}
93 : //
94 : //
95 : //
96 :
97 : class Quaternion : public Colvar {
98 : private:
99 : bool pbc;
100 : std::vector<double> value, masses, charges;
101 : std::vector<std::vector<Vector> > derivs;
102 : std::vector<Tensor> virial;
103 : public:
104 : static void registerKeywords( Keywords& keys );
105 : explicit Quaternion(const ActionOptions&);
106 : static void parseAtomList( const int& num, std::vector<AtomNumber>& t, ActionAtomistic* aa );
107 : static unsigned getModeAndSetupValues( ActionWithValue* av );
108 : // active methods:
109 : void calculate() override;
110 : static void calculateCV( const unsigned& mode, const std::vector<double>& masses, const std::vector<double>& charges,
111 : const std::vector<Vector>& pos, std::vector<double>& vals, std::vector<std::vector<Vector> >& derivs,
112 : std::vector<Tensor>& virial, const ActionAtomistic* aa );
113 : };
114 :
115 : typedef colvar::ColvarShortcut<Quaternion> QuaternionShortcut;
116 : PLUMED_REGISTER_ACTION(QuaternionShortcut,"QUATERNION")
117 : PLUMED_REGISTER_ACTION(Quaternion,"QUATERNION_SCALAR")
118 : typedef colvar::MultiColvarTemplate<Quaternion> QuaternionMulti;
119 : PLUMED_REGISTER_ACTION(QuaternionMulti,"QUATERNION_VECTOR")
120 :
121 86 : void Quaternion::registerKeywords( Keywords& keys ) {
122 86 : Colvar::registerKeywords( keys );
123 86 : keys.setDisplayName("QUATERNION");
124 172 : keys.add("atoms","ATOMS","the three atom that we are using to calculate the quaternion");
125 172 : keys.addOutputComponent("w","default","the real component of quaternion");
126 172 : keys.addOutputComponent("i","default","the i component of the quaternion");
127 172 : keys.addOutputComponent("j","default","the j component of the quaternion");
128 172 : keys.addOutputComponent("k","default","the k component of the quaternion");
129 172 : keys.add("hidden","NO_ACTION_LOG","suppresses printing from action on the log");
130 86 : }
131 :
132 5 : Quaternion::Quaternion(const ActionOptions&ao):
133 : PLUMED_COLVAR_INIT(ao),
134 5 : pbc(true),
135 5 : value(4),
136 5 : derivs(4),
137 10 : virial(4) {
138 25 : for(unsigned i=0; i<4; ++i) {
139 20 : derivs[i].resize(3);
140 : }
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 1318 : void Quaternion::parseAtomList( const int& num, std::vector<AtomNumber>& t, ActionAtomistic* aa ) {
157 2636 : aa->parseAtomList("ATOMS",num,t);
158 1318 : if( t.size()==3 ) {
159 1306 : aa->log.printf(" involving atoms %d %d %d\n",t[0].serial(),t[1].serial(),t[0].serial());
160 : }
161 1318 : }
162 :
163 17 : unsigned Quaternion::getModeAndSetupValues( ActionWithValue* av ) {
164 : // This sets up values that we can pass around in PLUMED
165 34 : av->addComponentWithDerivatives("w");
166 17 : av->componentIsNotPeriodic("w");
167 34 : av->addComponentWithDerivatives("i");
168 17 : av->componentIsNotPeriodic("i");
169 34 : av->addComponentWithDerivatives("j");
170 17 : av->componentIsNotPeriodic("j");
171 34 : av->addComponentWithDerivatives("k");
172 17 : av->componentIsNotPeriodic("k");
173 17 : return 0;
174 : }
175 :
176 45 : void Quaternion::calculate() {
177 45 : if(pbc) {
178 45 : makeWhole();
179 : }
180 :
181 45 : calculateCV( 0, masses, charges, getPositions(), value, derivs, virial, this );
182 225 : for(unsigned j=0; j<4; ++j) {
183 180 : Value* valuej=getPntrToComponent(j);
184 720 : for(unsigned i=0; i<3; ++i) {
185 540 : setAtomsDerivatives(valuej,i,derivs[j][i] );
186 : }
187 180 : setBoxDerivatives(valuej,virial[j]);
188 180 : valuej->set(value[j]);
189 : }
190 45 : }
191 :
192 : // calculator
193 5718 : void Quaternion::calculateCV( const unsigned& mode, const std::vector<double>& masses, const std::vector<double>& charges,
194 : const std::vector<Vector>& pos, std::vector<double>& vals, std::vector<std::vector<Vector> >& derivs,
195 : std::vector<Tensor>& virial, const ActionAtomistic* aa ) {
196 : //declarations
197 5718 : Vector vec1_comp = delta( pos[0], pos[1] ); //components between atom 1 and 2
198 5718 : Vector vec2_comp = delta( pos[0], pos[2] ); //components between atom 1 and 3
199 :
200 : ////////x-vector calculations///////
201 5718 : double magx = vec1_comp.modulo();
202 5718 : Vector xt = vec1_comp / magx;
203 5718 : std::vector<Tensor> dx(3);
204 5718 : double magx3= magx*magx*magx;
205 : //dx[i] - derivatives of atom i's coordinates
206 5718 : dx[0](0,0) = ( -(vec1_comp[1]*vec1_comp[1]+vec1_comp[2]*vec1_comp[2])/magx3 ); //dx[0]/dx0
207 5718 : dx[0](0,1) = ( vec1_comp[0]*vec1_comp[1]/magx3 ); // dx[0]/dy0
208 5718 : dx[0](0,2) = ( vec1_comp[0]*vec1_comp[2]/magx3 ); // dx[0]/dz0
209 5718 : dx[0](1,0) = ( vec1_comp[1]*vec1_comp[0]/magx3 ); // dx[1]/dx0
210 5718 : dx[0](1,1) = ( -(vec1_comp[0]*vec1_comp[0]+vec1_comp[2]*vec1_comp[2])/magx3 ); // dx[1]/dy0
211 5718 : dx[0](1,2) = ( vec1_comp[1]*vec1_comp[2]/magx3 ); //dx[1]/dz0
212 5718 : dx[0](2,0) = ( vec1_comp[2]*vec1_comp[0]/magx3 );//etc
213 5718 : dx[0](2,1) = ( vec1_comp[2]*vec1_comp[1]/magx3 );
214 5718 : dx[0](2,2) = ( -(vec1_comp[1]*vec1_comp[1]+vec1_comp[0]*vec1_comp[0])/magx3 );
215 :
216 5718 : dx[1](0,0) = ( (vec1_comp[1]*vec1_comp[1]+vec1_comp[2]*vec1_comp[2])/magx3 );//dx[0]/dx1
217 5718 : dx[1](0,1) = ( -vec1_comp[0]*vec1_comp[1]/magx3 );//dx[0]/dy1
218 5718 : dx[1](0,2) = ( -vec1_comp[0]*vec1_comp[2]/magx3 );
219 5718 : dx[1](1,0) = ( -vec1_comp[1]*vec1_comp[0]/magx3 );
220 5718 : dx[1](1,1) = ( (vec1_comp[0]*vec1_comp[0]+vec1_comp[2]*vec1_comp[2])/magx3 );
221 5718 : dx[1](1,2) = ( -vec1_comp[1]*vec1_comp[2]/magx3 );
222 5718 : dx[1](2,0) = ( -vec1_comp[2]*vec1_comp[0]/magx3 );
223 5718 : dx[1](2,1) = ( -vec1_comp[2]*vec1_comp[1]/magx3 );
224 5718 : dx[1](2,2) = ( (vec1_comp[1]*vec1_comp[1]+vec1_comp[0]*vec1_comp[0])/magx3 );
225 5718 : dx[2].zero();//not atom[2] terms present
226 :
227 : ////////y-vector calculations////////
228 : //project vec2_comp on to vec1_comp
229 : //first do dot product of unormalized x and unormed y, divided by magnitude of x^2
230 5718 : double dp = dotProduct( vec1_comp, vec2_comp );
231 : double magx2=magx*magx;
232 5718 : std::vector<Vector> fac_derivs(3);
233 5718 : double magx4=magx2*magx2, fac = dp/magx2; //fac meaning factor on front
234 5718 : fac_derivs[0] = (-vec2_comp - vec1_comp)/magx2 + 2*dp*vec1_comp / magx4;
235 5718 : fac_derivs[1] = (vec2_comp)/(magx2) - 2*dp*vec1_comp / magx4;
236 5718 : fac_derivs[2] = (vec1_comp)/(magx2); //atom 1, components x2,y2,z2
237 : //now multiply fac by unormed x, and subtract it from unormed y, then normalize
238 5718 : Vector yt = vec2_comp - fac*vec1_comp;
239 5718 : std::vector<Tensor> dy(3);
240 5718 : dy[0](0,0) = -1 - fac_derivs[0][0]*vec1_comp[0] + fac; // dy[0]/dx0
241 5718 : dy[0](0,1) = -fac_derivs[0][1]*vec1_comp[0]; // dy[0]/dy0
242 5718 : dy[0](0,2) = -fac_derivs[0][2]*vec1_comp[0];
243 5718 : dy[0](1,0) = -fac_derivs[0][0]*vec1_comp[1];
244 5718 : dy[0](1,1) = -1 - fac_derivs[0][1]*vec1_comp[1] + fac;
245 5718 : dy[0](1,2) = -fac_derivs[0][2]*vec1_comp[1];
246 5718 : dy[0](2,0) = -fac_derivs[0][0]*vec1_comp[2];
247 5718 : dy[0](2,1) = -fac_derivs[0][1]*vec1_comp[2];
248 5718 : dy[0](2,2) = -1 - fac_derivs[0][2]*vec1_comp[2] + fac;
249 :
250 5718 : dy[1](0,0) = -fac_derivs[1][0]*vec1_comp[0] - fac; //dy[0]/dx0
251 5718 : dy[1](0,1) = -fac_derivs[1][1]*vec1_comp[0];
252 5718 : dy[1](0,2) = -fac_derivs[1][2]*vec1_comp[0];
253 5718 : dy[1](1,0) = -fac_derivs[1][0]*vec1_comp[1];
254 5718 : dy[1](1,1) = -fac_derivs[1][1]*vec1_comp[1] - fac;
255 5718 : dy[1](1,2) = -fac_derivs[1][2]*vec1_comp[1];
256 5718 : dy[1](2,0) = -fac_derivs[1][0]*vec1_comp[2];
257 5718 : dy[1](2,1) = -fac_derivs[1][1]*vec1_comp[2];
258 5718 : dy[1](2,2) = -fac_derivs[1][2]*vec1_comp[2] - fac;
259 :
260 5718 : dy[2](0,0) = 1 - fac_derivs[2][0]*vec1_comp[0];//dy[0]/dx2
261 5718 : dy[2](0,1) = -fac_derivs[2][1]*vec1_comp[0];
262 5718 : dy[2](0,2) = -fac_derivs[2][2]*vec1_comp[0];
263 5718 : dy[2](1,0) = -fac_derivs[2][0]*vec1_comp[1];
264 5718 : dy[2](1,1) = 1 - fac_derivs[2][1]*vec1_comp[1];
265 5718 : dy[2](1,2) = -fac_derivs[2][2]*vec1_comp[1];
266 5718 : dy[2](2,0) = -fac_derivs[2][0]*vec1_comp[2];
267 5718 : dy[2](2,1) = -fac_derivs[2][1]*vec1_comp[2];
268 5718 : dy[2](2,2) = 1 - fac_derivs[2][2]*vec1_comp[2];
269 : //now normalize, and we have our y vector
270 5718 : double magy = yt.modulo();
271 5718 : double imagy = 1/magy, magy3 = magy*magy*magy;
272 5718 : Tensor abc;
273 22872 : for(unsigned i=0; i<3; ++i) {
274 17154 : abc.setRow(i, yt);
275 : }
276 5718 : Tensor abc_diag;
277 5718 : abc_diag.setRow(0, Vector(yt[0], 0, 0));
278 5718 : abc_diag.setRow(1, Vector(0, yt[1], 0));
279 5718 : abc_diag.setRow(2, Vector(0, 0, yt[2]));
280 5718 : Tensor abc_prod = matmul(abc_diag, abc);
281 22872 : for(unsigned i=0; i<3; ++i) {
282 17154 : dy[i] = dy[i]/magy - matmul(abc_prod, dy[i])/magy3;
283 : }
284 : //normalize now, derivatives are with respect to un-normalized y vector
285 5718 : yt = yt / magy;
286 :
287 : ///////z-vector calculations/////////
288 : //comparatively simple
289 5718 : Vector zt = crossProduct(xt,yt);
290 5718 : std::vector<Tensor> dz(3);
291 5718 : dz[0].setCol( 0, crossProduct( dx[0].getCol(0), yt ) + crossProduct( xt, dy[0].getCol(0) ) );
292 5718 : dz[0].setCol( 1, crossProduct( dx[0].getCol(1), yt ) + crossProduct( xt, dy[0].getCol(1) ) );
293 5718 : dz[0].setCol( 2, crossProduct( dx[0].getCol(2), yt ) + crossProduct( xt, dy[0].getCol(2) ) );
294 :
295 5718 : dz[1].setCol( 0, crossProduct( dx[1].getCol(0), yt ) + crossProduct( xt, dy[1].getCol(0) ) );
296 5718 : dz[1].setCol( 1, crossProduct( dx[1].getCol(1), yt ) + crossProduct( xt, dy[1].getCol(1) ) );
297 5718 : dz[1].setCol( 2, crossProduct( dx[1].getCol(2), yt ) + crossProduct( xt, dy[1].getCol(2) ) );
298 :
299 5718 : dz[2].setCol( 0, crossProduct( xt, dy[2].getCol(0) ) );
300 5718 : dz[2].setCol( 1, crossProduct( xt, dy[2].getCol(1) ) );
301 5718 : dz[2].setCol( 2, crossProduct( xt, dy[2].getCol(2) ) );
302 :
303 : //for debugging frame values
304 : //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]);
305 :
306 : //for bebuffing derivatives
307 : //aa->log.printf("x1 x2 x3 y1 y2 y3 z1 z2 z3\n");
308 : //for (int i=0; i<3; i++){
309 : //for (int j=0;j<3;j++){
310 : //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));
311 : //}
312 : //}
313 : //
314 :
315 : //the above 9 components form an orthonormal basis, centered on the molecule in question
316 : //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
317 : //the inverse is simply the transpose
318 :
319 :
320 : //[x[0] x[1] x[2]]
321 : //[y[0] y[1] y[2]]
322 : //[z[0] z[1] z[2]]
323 : //QUICKFIX to transpose basis
324 5718 : Vector x(xt[0],yt[0],zt[0]);
325 5718 : Vector y(xt[1],yt[1],zt[1]);
326 5718 : Vector z(xt[2],yt[2],zt[2]);
327 :
328 : //likewise transposing the tensors into proper form
329 5718 : std::vector<Tensor> tdx(3);
330 5718 : std::vector<Tensor> tdy(3);
331 5718 : std::vector<Tensor> tdz(3);
332 22872 : for (int i=0; i<3; ++i) {
333 17154 : tdx[i].setRow(0, dx[i].getRow(0));
334 17154 : tdx[i].setRow(1, dy[i].getRow(0));
335 17154 : tdx[i].setRow(2, dz[i].getRow(0));
336 :
337 17154 : tdy[i].setRow(0, dx[i].getRow(1));
338 17154 : tdy[i].setRow(1, dy[i].getRow(1));
339 17154 : tdy[i].setRow(2, dz[i].getRow(1));
340 :
341 17154 : tdz[i].setRow(0, dx[i].getRow(2));
342 17154 : tdz[i].setRow(1, dy[i].getRow(2));
343 17154 : tdz[i].setRow(2, dz[i].getRow(2));
344 : }
345 :
346 : //convert to quaternion
347 5718 : double tr = x[0] + y[1] + z[2] + 1; //trace of the rotation matrix + 1
348 5718 : std::vector<Vector> dS(3);
349 5718 : if (tr > 1.0E-8) { //to avoid numerical instability
350 5718 : double S = 1/(sqrt(tr) * 2); // S=4*qw
351 22872 : for(unsigned i=0; i<3; ++i) {
352 17154 : dS[i] = (-2*S*S*S)*(tdx[i].getRow(0) + tdy[i].getRow(1) + tdz[i].getRow(2));
353 : }
354 :
355 5718 : vals[0] = 0.25 / S;
356 22872 : for(unsigned i=0; i<3; ++i) {
357 17154 : derivs[0][i] =-0.25*dS[i]/(S*S);
358 : }
359 :
360 5718 : vals[1] = (z[1] - y[2]) * S;
361 22872 : for(unsigned i=0; i<3; ++i) {
362 17154 : derivs[1][i] = (S)*(tdz[i].getRow(1) - tdy[i].getRow(2)) + (z[1]-y[2])*dS[i];
363 : }
364 :
365 5718 : vals[2] = (x[2] - z[0]) * S;
366 22872 : for(unsigned i=0; i<3; ++i) {
367 17154 : derivs[2][i] = (S)*(tdx[i].getRow(2) - tdz[i].getRow(0)) + (x[2]-z[0])*dS[i];
368 : }
369 :
370 5718 : vals[3] = (y[0] - x[1]) * S;
371 22872 : for(unsigned i=0; i<3; ++i) {
372 17154 : derivs[3][i] = (S)*(tdy[i].getRow(0) - tdx[i].getRow(1)) + (y[0]-x[1])*dS[i];
373 : }
374 0 : } else if ((x[0] > y[1])&(x[0] > z[2])) {
375 0 : double S = sqrt(1.0 + x[0] - y[1] - z[2]) * 2; // S=4*qx
376 0 : for(unsigned i=0; i<3; ++i) {
377 0 : dS[i] = (2/S)*(tdx[i].getRow(0) - tdy[i].getRow(1) - tdz[i].getRow(2));
378 : }
379 :
380 0 : vals[0] = (z[1] - y[2]) / S;
381 0 : for(unsigned i=0; i<3; ++i) {
382 0 : derivs[0][i] = (1/S)*(tdz[i].getRow(1) - tdy[i].getRow(2)) - (vals[0]/S)*dS[i];
383 : }
384 :
385 0 : vals[1] = 0.25 * S;
386 0 : for(unsigned i=0; i<3; ++i) {
387 0 : derivs[1][i] =0.25*dS[i];
388 : }
389 :
390 0 : vals[2] = (x[1] + y[0]) / S;
391 0 : for(unsigned i=0; i<3; ++i) {
392 0 : derivs[2][i] = (1/S)*(tdx[i].getRow(1) + tdy[i].getRow(0)) - (vals[2]/S)*dS[i];
393 : }
394 :
395 0 : vals[3] = (x[2] + z[0]) / S;
396 0 : for(unsigned i=0; i<3; ++i) {
397 0 : derivs[3][i] = (1/S)*(tdx[i].getRow(2) + tdz[i].getRow(0)) - (vals[3]/S)*dS[i];
398 : }
399 0 : } else if (y[1] > z[2]) {
400 0 : double S = sqrt(1.0 + y[1] - x[0] - z[2]) * 2; // S=4*qy
401 0 : for(unsigned i=0; i<3; ++i) {
402 0 : dS[i] = (2/S)*( -tdx[i].getRow(0) + tdy[i].getRow(1) - tdz[i].getRow(2));
403 : }
404 :
405 :
406 0 : vals[0] = (x[2] - z[0]) / S;
407 0 : for(unsigned i=0; i<3; ++i) {
408 0 : derivs[0][i] = (1/S)*(tdx[i].getRow(2) - tdz[i].getRow(0)) - (vals[0]/S)*dS[i];
409 : }
410 :
411 0 : vals[1] = (x[1] + y[0]) / S;
412 0 : for(unsigned i=0; i<3; ++i) {
413 0 : derivs[1][i] = (1/S)*(tdx[i].getRow(1) + tdy[i].getRow(0)) - (vals[1]/S)*dS[i];
414 : }
415 :
416 0 : vals[2] = 0.25 * S;
417 0 : for(unsigned i=0; i<3; ++i) {
418 0 : derivs[2][i] =0.25*dS[i];
419 : }
420 :
421 0 : vals[3] = (y[2] + z[1]) / S;
422 0 : for(unsigned i=0; i<3; ++i) {
423 0 : derivs[3][i] = (1/S)*(tdy[i].getRow(2) + tdz[i].getRow(1)) - (vals[3]/S)*dS[i];
424 : }
425 : } else {
426 0 : double S = sqrt(1.0 + z[2] - x[0] - y[1]) * 2; // S=4*qz
427 0 : for(unsigned i=0; i<3; ++i) {
428 0 : dS[i] = (2/S)*(-tdx[i].getRow(0) - tdy[i].getRow(1) + tdz[i].getRow(2));
429 : }
430 :
431 :
432 0 : vals[0] = (y[0] - x[1]) / S;
433 0 : for(unsigned i=0; i<3; ++i) {
434 0 : derivs[0][i] = (1/S)*(tdy[i].getRow(0) - tdx[i].getRow(1)) - (vals[0]/S)*dS[i];
435 : }
436 :
437 0 : vals[1] = (x[2] + z[0]) / S;
438 0 : for(unsigned i=0; i<3; ++i) {
439 0 : derivs[1][i] = (1/S)*(tdx[i].getRow(2) + tdz[i].getRow(0)) - (vals[1]/S)*dS[i];
440 : }
441 :
442 0 : vals[2] = (y[2] + z[1]) / S;
443 0 : for(unsigned i=0; i<3; ++i) {
444 0 : derivs[2][i] = (1/S)*(tdy[i].getRow(2) + tdz[i].getRow(1)) - (vals[2]/S)*dS[i];
445 : }
446 :
447 0 : vals[3] = 0.25 * S;
448 0 : for(unsigned i=0; i<3; ++i) {
449 0 : derivs[3][i] =0.25*dS[i];
450 : }
451 : }
452 5718 : setBoxDerivativesNoPbc( pos, derivs, virial );
453 :
454 5718 : }
455 :
456 : }
457 : }
458 :
459 :
460 :
|
