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