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 ); keys.setDisplayName("QUATERNION");
123 172 : keys.add("atoms","ATOMS","the three atom that we are using to calculate the quaternion");
124 172 : keys.addOutputComponent("w","default","the real component of quaternion");
125 172 : keys.addOutputComponent("i","default","the i component of the quaternion");
126 172 : keys.addOutputComponent("j","default","the j component of the quaternion");
127 172 : keys.addOutputComponent("k","default","the k component of the quaternion");
128 172 : keys.add("hidden","NO_ACTION_LOG","suppresses printing from action on the log");
129 86 : }
130 :
131 5 : Quaternion::Quaternion(const ActionOptions&ao):
132 : PLUMED_COLVAR_INIT(ao),
133 5 : pbc(true),
134 5 : value(4),
135 5 : derivs(4),
136 10 : virial(4)
137 : {
138 25 : for(unsigned i=0; i<4; ++i) derivs[i].resize(3);
139 : std::vector<AtomNumber> atoms;
140 5 : parseAtomList(-1,atoms,this);
141 5 : if(atoms.size()!=3) error("Number of specified atoms should be 3");
142 : // please note, I do NO checking if these atoms are in the same molecule at all, so be careful in your input
143 :
144 5 : bool nopbc=!pbc;
145 5 : parseFlag("NOPBC",nopbc);
146 5 : pbc=!nopbc;
147 :
148 5 : unsigned mode = getModeAndSetupValues( this );
149 5 : requestAtoms(atoms);
150 5 : }
151 :
152 1318 : void Quaternion::parseAtomList( const int& num, std::vector<AtomNumber>& t, ActionAtomistic* aa ) {
153 2636 : aa->parseAtomList("ATOMS",num,t);
154 1318 : if( t.size()==3 ) aa->log.printf(" involving atoms %d %d %d\n",t[0].serial(),t[1].serial(),t[0].serial());
155 1318 : }
156 :
157 17 : unsigned Quaternion::getModeAndSetupValues( ActionWithValue* av ) {
158 : // This sets up values that we can pass around in PLUMED
159 34 : av->addComponentWithDerivatives("w"); av->componentIsNotPeriodic("w");
160 34 : av->addComponentWithDerivatives("i"); av->componentIsNotPeriodic("i");
161 34 : av->addComponentWithDerivatives("j"); av->componentIsNotPeriodic("j");
162 34 : av->addComponentWithDerivatives("k"); av->componentIsNotPeriodic("k");
163 17 : return 0;
164 : }
165 :
166 45 : void Quaternion::calculate() {
167 45 : if(pbc) makeWhole();
168 :
169 45 : calculateCV( 0, masses, charges, getPositions(), value, derivs, virial, this );
170 225 : for(unsigned j=0; j<4; ++j) {
171 180 : Value* valuej=getPntrToComponent(j);
172 720 : for(unsigned i=0; i<3; ++i) setAtomsDerivatives(valuej,i,derivs[j][i] );
173 180 : setBoxDerivatives(valuej,virial[j]);
174 180 : valuej->set(value[j]);
175 : }
176 45 : }
177 :
178 : // calculator
179 5718 : void Quaternion::calculateCV( const unsigned& mode, const std::vector<double>& masses, const std::vector<double>& charges,
180 : const std::vector<Vector>& pos, std::vector<double>& vals, std::vector<std::vector<Vector> >& derivs,
181 : std::vector<Tensor>& virial, const ActionAtomistic* aa ) {
182 : //declarations
183 5718 : Vector vec1_comp = delta( pos[0], pos[1] ); //components between atom 1 and 2
184 5718 : Vector vec2_comp = delta( pos[0], pos[2] ); //components between atom 1 and 3
185 :
186 : ////////x-vector calculations///////
187 5718 : double magx = vec1_comp.modulo(); Vector xt = vec1_comp / magx;
188 5718 : std::vector<Tensor> dx(3); double magx3= magx*magx*magx;
189 : //dx[i] - derivatives of atom i's coordinates
190 5718 : dx[0](0,0) = ( -(vec1_comp[1]*vec1_comp[1]+vec1_comp[2]*vec1_comp[2])/magx3 ); //dx[0]/dx0
191 5718 : dx[0](0,1) = ( vec1_comp[0]*vec1_comp[1]/magx3 ); // dx[0]/dy0
192 5718 : dx[0](0,2) = ( vec1_comp[0]*vec1_comp[2]/magx3 ); // dx[0]/dz0
193 5718 : dx[0](1,0) = ( vec1_comp[1]*vec1_comp[0]/magx3 ); // dx[1]/dx0
194 5718 : dx[0](1,1) = ( -(vec1_comp[0]*vec1_comp[0]+vec1_comp[2]*vec1_comp[2])/magx3 ); // dx[1]/dy0
195 5718 : dx[0](1,2) = ( vec1_comp[1]*vec1_comp[2]/magx3 ); //dx[1]/dz0
196 5718 : dx[0](2,0) = ( vec1_comp[2]*vec1_comp[0]/magx3 );//etc
197 5718 : dx[0](2,1) = ( vec1_comp[2]*vec1_comp[1]/magx3 );
198 5718 : dx[0](2,2) = ( -(vec1_comp[1]*vec1_comp[1]+vec1_comp[0]*vec1_comp[0])/magx3 );
199 :
200 5718 : dx[1](0,0) = ( (vec1_comp[1]*vec1_comp[1]+vec1_comp[2]*vec1_comp[2])/magx3 );//dx[0]/dx1
201 5718 : dx[1](0,1) = ( -vec1_comp[0]*vec1_comp[1]/magx3 );//dx[0]/dy1
202 5718 : dx[1](0,2) = ( -vec1_comp[0]*vec1_comp[2]/magx3 );
203 5718 : dx[1](1,0) = ( -vec1_comp[1]*vec1_comp[0]/magx3 );
204 5718 : dx[1](1,1) = ( (vec1_comp[0]*vec1_comp[0]+vec1_comp[2]*vec1_comp[2])/magx3 );
205 5718 : dx[1](1,2) = ( -vec1_comp[1]*vec1_comp[2]/magx3 );
206 5718 : dx[1](2,0) = ( -vec1_comp[2]*vec1_comp[0]/magx3 );
207 5718 : dx[1](2,1) = ( -vec1_comp[2]*vec1_comp[1]/magx3 );
208 5718 : dx[1](2,2) = ( (vec1_comp[1]*vec1_comp[1]+vec1_comp[0]*vec1_comp[0])/magx3 );
209 5718 : dx[2].zero();//not atom[2] terms present
210 :
211 : ////////y-vector calculations////////
212 : //project vec2_comp on to vec1_comp
213 : //first do dot product of unormalized x and unormed y, divided by magnitude of x^2
214 5718 : double dp = dotProduct( vec1_comp, vec2_comp ); double magx2=magx*magx;
215 5718 : std::vector<Vector> fac_derivs(3); double magx4=magx2*magx2, fac = dp/magx2; //fac meaning factor on front
216 5718 : fac_derivs[0] = (-vec2_comp - vec1_comp)/magx2 + 2*dp*vec1_comp / magx4;
217 5718 : fac_derivs[1] = (vec2_comp)/(magx2) - 2*dp*vec1_comp / magx4;
218 5718 : fac_derivs[2] = (vec1_comp)/(magx2); //atom 1, components x2,y2,z2
219 : //now multiply fac by unormed x, and subtract it from unormed y, then normalize
220 5718 : Vector yt = vec2_comp - fac*vec1_comp; std::vector<Tensor> dy(3);
221 5718 : dy[0](0,0) = -1 - fac_derivs[0][0]*vec1_comp[0] + fac; // dy[0]/dx0
222 5718 : dy[0](0,1) = -fac_derivs[0][1]*vec1_comp[0]; // dy[0]/dy0
223 5718 : dy[0](0,2) = -fac_derivs[0][2]*vec1_comp[0];
224 5718 : dy[0](1,0) = -fac_derivs[0][0]*vec1_comp[1];
225 5718 : dy[0](1,1) = -1 - fac_derivs[0][1]*vec1_comp[1] + fac;
226 5718 : dy[0](1,2) = -fac_derivs[0][2]*vec1_comp[1];
227 5718 : dy[0](2,0) = -fac_derivs[0][0]*vec1_comp[2];
228 5718 : dy[0](2,1) = -fac_derivs[0][1]*vec1_comp[2];
229 5718 : dy[0](2,2) = -1 - fac_derivs[0][2]*vec1_comp[2] + fac;
230 :
231 5718 : dy[1](0,0) = -fac_derivs[1][0]*vec1_comp[0] - fac; //dy[0]/dx0
232 5718 : dy[1](0,1) = -fac_derivs[1][1]*vec1_comp[0];
233 5718 : dy[1](0,2) = -fac_derivs[1][2]*vec1_comp[0];
234 5718 : dy[1](1,0) = -fac_derivs[1][0]*vec1_comp[1];
235 5718 : dy[1](1,1) = -fac_derivs[1][1]*vec1_comp[1] - fac;
236 5718 : dy[1](1,2) = -fac_derivs[1][2]*vec1_comp[1];
237 5718 : dy[1](2,0) = -fac_derivs[1][0]*vec1_comp[2];
238 5718 : dy[1](2,1) = -fac_derivs[1][1]*vec1_comp[2];
239 5718 : dy[1](2,2) = -fac_derivs[1][2]*vec1_comp[2] - fac;
240 :
241 5718 : dy[2](0,0) = 1 - fac_derivs[2][0]*vec1_comp[0];//dy[0]/dx2
242 5718 : dy[2](0,1) = -fac_derivs[2][1]*vec1_comp[0];
243 5718 : dy[2](0,2) = -fac_derivs[2][2]*vec1_comp[0];
244 5718 : dy[2](1,0) = -fac_derivs[2][0]*vec1_comp[1];
245 5718 : dy[2](1,1) = 1 - fac_derivs[2][1]*vec1_comp[1];
246 5718 : dy[2](1,2) = -fac_derivs[2][2]*vec1_comp[1];
247 5718 : dy[2](2,0) = -fac_derivs[2][0]*vec1_comp[2];
248 5718 : dy[2](2,1) = -fac_derivs[2][1]*vec1_comp[2];
249 5718 : dy[2](2,2) = 1 - fac_derivs[2][2]*vec1_comp[2];
250 : //now normalize, and we have our y vector
251 5718 : double magy = yt.modulo(); double imagy = 1/magy, magy3 = magy*magy*magy;
252 22872 : Tensor abc; for(unsigned i=0; i<3; ++i) abc.setRow(i, yt);
253 5718 : Tensor abc_diag; abc_diag.setRow(0, Vector(yt[0], 0, 0)); abc_diag.setRow(1, Vector(0, yt[1], 0)); abc_diag.setRow(2, Vector(0, 0, yt[2]));
254 5718 : Tensor abc_prod = matmul(abc_diag, abc);
255 22872 : for(unsigned i=0; i<3; ++i) dy[i] = dy[i]/magy - matmul(abc_prod, dy[i])/magy3;
256 : //normalize now, derivatives are with respect to un-normalized y vector
257 5718 : yt = yt / magy;
258 :
259 : ///////z-vector calculations/////////
260 : //comparatively simple
261 5718 : Vector zt = crossProduct(xt,yt); std::vector<Tensor> dz(3);
262 5718 : dz[0].setCol( 0, crossProduct( dx[0].getCol(0), yt ) + crossProduct( xt, dy[0].getCol(0) ) );
263 5718 : dz[0].setCol( 1, crossProduct( dx[0].getCol(1), yt ) + crossProduct( xt, dy[0].getCol(1) ) );
264 5718 : dz[0].setCol( 2, crossProduct( dx[0].getCol(2), yt ) + crossProduct( xt, dy[0].getCol(2) ) );
265 :
266 5718 : dz[1].setCol( 0, crossProduct( dx[1].getCol(0), yt ) + crossProduct( xt, dy[1].getCol(0) ) );
267 5718 : dz[1].setCol( 1, crossProduct( dx[1].getCol(1), yt ) + crossProduct( xt, dy[1].getCol(1) ) );
268 5718 : dz[1].setCol( 2, crossProduct( dx[1].getCol(2), yt ) + crossProduct( xt, dy[1].getCol(2) ) );
269 :
270 5718 : dz[2].setCol( 0, crossProduct( xt, dy[2].getCol(0) ) );
271 5718 : dz[2].setCol( 1, crossProduct( xt, dy[2].getCol(1) ) );
272 5718 : dz[2].setCol( 2, crossProduct( xt, dy[2].getCol(2) ) );
273 :
274 : //for debugging frame values
275 : //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]);
276 :
277 : //for bebuffing derivatives
278 : //aa->log.printf("x1 x2 x3 y1 y2 y3 z1 z2 z3\n");
279 : //for (int i=0; i<3; i++){
280 : //for (int j=0;j<3;j++){
281 : //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));
282 : //}
283 : //}
284 : //
285 :
286 : //the above 9 components form an orthonormal basis, centered on the molecule in question
287 : //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
288 : //the inverse is simply the transpose
289 :
290 :
291 : //[x[0] x[1] x[2]]
292 : //[y[0] y[1] y[2]]
293 : //[z[0] z[1] z[2]]
294 : //QUICKFIX to transpose basis
295 5718 : Vector x(xt[0],yt[0],zt[0]);
296 5718 : Vector y(xt[1],yt[1],zt[1]);
297 5718 : Vector z(xt[2],yt[2],zt[2]);
298 :
299 : //likewise transposing the tensors into proper form
300 5718 : std::vector<Tensor> tdx(3);
301 5718 : std::vector<Tensor> tdy(3);
302 5718 : std::vector<Tensor> tdz(3);
303 22872 : for (int i=0; i<3; ++i) {
304 17154 : tdx[i].setRow(0, dx[i].getRow(0));
305 17154 : tdx[i].setRow(1, dy[i].getRow(0));
306 17154 : tdx[i].setRow(2, dz[i].getRow(0));
307 :
308 17154 : tdy[i].setRow(0, dx[i].getRow(1));
309 17154 : tdy[i].setRow(1, dy[i].getRow(1));
310 17154 : tdy[i].setRow(2, dz[i].getRow(1));
311 :
312 17154 : tdz[i].setRow(0, dx[i].getRow(2));
313 17154 : tdz[i].setRow(1, dy[i].getRow(2));
314 17154 : tdz[i].setRow(2, dz[i].getRow(2));
315 : }
316 :
317 : //convert to quaternion
318 5718 : double tr = x[0] + y[1] + z[2] + 1; //trace of the rotation matrix + 1
319 5718 : std::vector<Vector> dS(3);
320 5718 : if (tr > 1.0E-8) { //to avoid numerical instability
321 5718 : double S = 1/(sqrt(tr) * 2); // S=4*qw
322 22872 : for(unsigned i=0; i<3; ++i) dS[i] = (-2*S*S*S)*(tdx[i].getRow(0) + tdy[i].getRow(1) + tdz[i].getRow(2));
323 :
324 5718 : vals[0] = 0.25 / S;
325 22872 : for(unsigned i=0; i<3; ++i) derivs[0][i] =-0.25*dS[i]/(S*S);
326 :
327 5718 : vals[1] = (z[1] - y[2]) * S;
328 22872 : for(unsigned i=0; i<3; ++i) derivs[1][i] = (S)*(tdz[i].getRow(1) - tdy[i].getRow(2)) + (z[1]-y[2])*dS[i];
329 :
330 5718 : vals[2] = (x[2] - z[0]) * S;
331 22872 : for(unsigned i=0; i<3; ++i) derivs[2][i] = (S)*(tdx[i].getRow(2) - tdz[i].getRow(0)) + (x[2]-z[0])*dS[i];
332 :
333 5718 : vals[3] = (y[0] - x[1]) * S;
334 22872 : for(unsigned i=0; i<3; ++i) derivs[3][i] = (S)*(tdy[i].getRow(0) - tdx[i].getRow(1)) + (y[0]-x[1])*dS[i];
335 : }
336 0 : else if ((x[0] > y[1])&(x[0] > z[2])) {
337 0 : float S = sqrt(1.0 + x[0] - y[1] - z[2]) * 2; // S=4*qx
338 0 : for(unsigned i=0; i<3; ++i) dS[i] = (2/S)*(tdx[i].getRow(0) - tdy[i].getRow(1) - tdz[i].getRow(2));
339 :
340 0 : vals[0] = (z[1] - y[2]) / S;
341 0 : for(unsigned i=0; i<3; ++i) derivs[0][i] = (1/S)*(tdz[i].getRow(1) - tdy[i].getRow(2)) - (vals[0]/S)*dS[i];
342 :
343 0 : vals[1] = 0.25 * S;
344 0 : for(unsigned i=0; i<3; ++i) derivs[1][i] =0.25*dS[i];
345 :
346 0 : vals[2] = (x[1] + y[0]) / S;
347 0 : for(unsigned i=0; i<3; ++i) derivs[2][i] = (1/S)*(tdx[i].getRow(1) + tdy[i].getRow(0)) - (vals[2]/S)*dS[i];
348 :
349 0 : vals[3] = (x[2] + z[0]) / S;
350 0 : for(unsigned i=0; i<3; ++i) derivs[3][i] = (1/S)*(tdx[i].getRow(2) + tdz[i].getRow(0)) - (vals[3]/S)*dS[i];
351 : }
352 0 : else if (y[1] > z[2]) {
353 0 : float S = sqrt(1.0 + y[1] - x[0] - z[2]) * 2; // S=4*qy
354 0 : for(unsigned i=0; i<3; ++i) dS[i] = (2/S)*( -tdx[i].getRow(0) + tdy[i].getRow(1) - tdz[i].getRow(2));
355 :
356 :
357 0 : vals[0] = (x[2] - z[0]) / S;
358 0 : for(unsigned i=0; i<3; ++i) derivs[0][i] = (1/S)*(tdx[i].getRow(2) - tdz[i].getRow(0)) - (vals[0]/S)*dS[i];
359 :
360 0 : vals[1] = (x[1] + y[0]) / S;
361 0 : for(unsigned i=0; i<3; ++i) derivs[1][i] = (1/S)*(tdx[i].getRow(1) + tdy[i].getRow(0)) - (vals[1]/S)*dS[i];
362 :
363 0 : vals[2] = 0.25 * S;
364 0 : for(unsigned i=0; i<3; ++i) derivs[2][i] =0.25*dS[i];
365 :
366 0 : vals[3] = (y[2] + z[1]) / S;
367 0 : for(unsigned i=0; i<3; ++i) derivs[3][i] = (1/S)*(tdy[i].getRow(2) + tdz[i].getRow(1)) - (vals[3]/S)*dS[i];
368 : }
369 : else {
370 0 : float S = sqrt(1.0 + z[2] - x[0] - y[1]) * 2; // S=4*qz
371 0 : for(unsigned i=0; i<3; ++i) dS[i] = (2/S)*(-tdx[i].getRow(0) - tdy[i].getRow(1) + tdz[i].getRow(2));
372 :
373 :
374 0 : vals[0] = (y[0] - x[1]) / S;
375 0 : for(unsigned i=0; i<3; ++i) derivs[0][i] = (1/S)*(tdy[i].getRow(0) - tdx[i].getRow(1)) - (vals[0]/S)*dS[i];
376 :
377 0 : vals[1] = (x[2] + z[0]) / S;
378 0 : for(unsigned i=0; i<3; ++i) derivs[1][i] = (1/S)*(tdx[i].getRow(2) + tdz[i].getRow(0)) - (vals[1]/S)*dS[i];
379 :
380 0 : vals[2] = (y[2] + z[1]) / S;
381 0 : for(unsigned i=0; i<3; ++i) derivs[2][i] = (1/S)*(tdy[i].getRow(2) + tdz[i].getRow(1)) - (vals[2]/S)*dS[i];
382 :
383 0 : vals[3] = 0.25 * S;
384 0 : for(unsigned i=0; i<3; ++i) derivs[3][i] =0.25*dS[i];
385 : }
386 5718 : setBoxDerivativesNoPbc( pos, derivs, virial );
387 :
388 5718 : }
389 :
390 : }
391 : }
392 :
393 :
394 :
|
