1 /* ----------------------------------------------------------------------
2 This is the
3
4 ██╗ ██╗ ██████╗ ██████╗ ██████╗ ██╗ ██╗████████╗███████╗
5 ██║ ██║██╔════╝ ██╔════╝ ██╔════╝ ██║ ██║╚══██╔══╝██╔════╝
6 ██║ ██║██║ ███╗██║ ███╗██║ ███╗███████║ ██║ ███████╗
7 ██║ ██║██║ ██║██║ ██║██║ ██║██╔══██║ ██║ ╚════██║
8 ███████╗██║╚██████╔╝╚██████╔╝╚██████╔╝██║ ██║ ██║ ███████║
9 ╚══════╝╚═╝ ╚═════╝ ╚═════╝ ╚═════╝ ╚═╝ ╚═╝ ╚═╝ ╚══════╝®
10
11 DEM simulation engine, released by
12 DCS Computing Gmbh, Linz, Austria
13 http://www.dcs-computing.com, office@dcs-computing.com
14
15 LIGGGHTS® is part of CFDEM®project:
16 http://www.liggghts.com | http://www.cfdem.com
17
18 Core developer and main author:
19 Christoph Kloss, christoph.kloss@dcs-computing.com
20
21 LIGGGHTS® is open-source, distributed under the terms of the GNU Public
22 License, version 2 or later. It is distributed in the hope that it will
23 be useful, but WITHOUT ANY WARRANTY; without even the implied warranty
24 of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. You should have
25 received a copy of the GNU General Public License along with LIGGGHTS®.
26 If not, see http://www.gnu.org/licenses . See also top-level README
27 and LICENSE files.
28
29 LIGGGHTS® and CFDEM® are registered trade marks of DCS Computing GmbH,
30 the producer of the LIGGGHTS® software and the CFDEM®coupling software
31 See http://www.cfdem.com/terms-trademark-policy for details.
32
33 -------------------------------------------------------------------------
34 Contributing author and copyright for this file:
35 Alexander Podlozhnyuk, DCS Computing GmbH, Linz
36 Christoph Kloss, DCS Computing GmbH, Linz
37
38 Copyright 2015- DCS Computing GmbH, Linz
39 ------------------------------------------------------------------------- */
40
41 #include <cmath>
42 #include <stdio.h>
43 #include <string.h>
44 #include "fix_nve_asphere_base.h"
45 #include "atom.h"
46 #include "atom_vec.h"
47 #include "update.h"
48 #include "respa.h"
49 #include "force.h"
50 #include "error.h"
51 #include "domain.h"
52 #include "math_extra_liggghts_nonspherical.h"
53 #include "fix_property_atom.h"
54
55 /* ---------------------------------------------------------------------- */
56
FixNVEAsphereBase(LAMMPS * lmp,int narg,char ** arg)57 FixNVEAsphereBase::FixNVEAsphereBase(LAMMPS *lmp, int narg, char **arg) :
58 FixNVE(lmp, narg, arg)
59 {
60 if (narg < 3) error->all(FLERR,"Illegal fix nve/superquadric command");
61
62 time_integrate = 1;
63
64 // process extra keywords
65
66 integration_scheme = 1;
67
68 int iarg = 3;
69 while (iarg < narg) {
70 if (strcmp(arg[iarg],"integration_scheme") == 0) {
71 integration_scheme = force->numeric(FLERR,arg[iarg+1]);
72 iarg += 2;
73 }
74 else error->fix_error(FLERR,this,"unknown keyword");
75 }
76
77 hdtorque = NULL;
78 orientation = NULL;
79 ksl_rotation = NULL;
80 couple_fix_id = -1;
81 for(int ifix = 0; ifix < modify->nfix; ifix++) {
82 if(strcmp(modify->fix[ifix]->style, "couple/cfd/force/implicit")==0) {
83 couple_fix_id = ifix;
84 break;
85 }
86 }
87 }
88
89 /* ---------------------------------------------------------------------- */
90
init()91 void FixNVEAsphereBase::init()
92 {
93 FixNVE::init();
94
95 // error checks might go here
96 }
97
98 /* ---------------------------------------------------------------------- */
99
100 //dynamic Euler equations for angular velocity in body's principal axes
dynamic_euler(double * wbody,double * tbody,double * inertia,double * result)101 void FixNVEAsphereBase::dynamic_euler(double *wbody, double *tbody, double *inertia, double *result)
102 {
103 result[0] = tbody[0] / inertia[0] + wbody[1]*wbody[2]*((inertia[1] - inertia[2]) / inertia[0]);
104 result[1] = tbody[1] / inertia[1] + wbody[2]*wbody[0]*((inertia[2] - inertia[0]) / inertia[1]);
105 result[2] = tbody[2] / inertia[2] + wbody[0]*wbody[1]*((inertia[0] - inertia[1]) / inertia[2]);
106 }
107
108 /* ---------------------------------------------------------------------- */
109
integrate_dynamic_euler(double dt,double * wbody,double * tbody,double * inertia)110 void FixNVEAsphereBase::integrate_dynamic_euler(double dt, double *wbody, double *tbody, double *inertia)
111 {
112 double omega_der[3];
113 double tol = 1e-12;
114 if(LAMMPS_NS::vectorMag3D(wbody)*dt > 1.0)
115 error->one(FLERR, "Timestep is too big for rotation integration!");
116 dynamic_euler(wbody, tbody, inertia, omega_der);
117
118 double omega_half_prev[3], delta[3];
119 double omega_half[] = {0.0, 0.0, 0.0};
120 while(1) {
121 LAMMPS_NS::vectorCopy3D(omega_half, omega_half_prev);
122 omega_half[0] = wbody[0] + 0.5*dt*omega_der[0];
123 omega_half[1] = wbody[1] + 0.5*dt*omega_der[1];
124 omega_half[2] = wbody[2] + 0.5*dt*omega_der[2];
125 LAMMPS_NS::vectorSubtract3D(omega_half_prev, omega_half, delta);
126 double omega_half_mag = LAMMPS_NS::vectorMag3D(omega_half);
127 if(omega_half_mag > 0.0) {
128 double eps = LAMMPS_NS::vectorMag3D(delta) / omega_half_mag;
129 if(eps < tol)
130 break;
131 dynamic_euler(omega_half, tbody, inertia, omega_der);
132 } else
133 break;
134 }
135 wbody[0] += dt*omega_der[0];
136 wbody[1] += dt*omega_der[1];
137 wbody[2] += dt*omega_der[2];
138 }
139
140 /* ---------------------------------------------------------------------- */
141
integrate_quaternion(double dtq,double * quat,double * wbody)142 void FixNVEAsphereBase::integrate_quaternion(double dtq, double *quat, double *wbody)
143 {
144 double q2[4];
145 q2[0] = 1.0;
146 q2[1] = -wbody[0] * dtq;
147 q2[2] = -wbody[1] * dtq;
148 q2[3] = -wbody[2] * dtq;
149 MathExtraLiggghtsNonspherical::invquat(q2); //using implicit scheme
150
151 double q_temp[4];
152 MathExtra::quatquat(quat, q2, q_temp);
153 vectorCopy4D(q_temp, quat);
154 MathExtra::qnormalize(quat);
155 }
156
157 /* ---------------------------------------------------------------------- */
158
update_hdtorque(int i,double * rotation_matrix,double * omegaOld,double * omegaNew)159 void FixNVEAsphereBase::update_hdtorque(int i, double *rotation_matrix, double *omegaOld, double *omegaNew)
160 {
161 if(ksl_rotation and hdtorque) {
162 double deltaHydrotorquePrime[3], deltaHydrotorque[3];
163 for(int dirI = 0; dirI < 3; dirI++)
164 deltaHydrotorquePrime[dirI] = ksl_rotation[i][dirI]*(omegaOld[dirI]-omegaNew[dirI]);
165 MathExtraLiggghtsNonspherical::matvec(rotation_matrix, deltaHydrotorquePrime, deltaHydrotorque);
166 for(int dirI=0; dirI<3; dirI++)
167 hdtorque[i][dirI] += deltaHydrotorque[dirI];
168 }
169 }
170
171 /* ---------------------------------------------------------------------- */
172
implicitRotationUpdate(double deltaT,double * inertia,double * angMom,double * torque,double * KslRot,double * omegaNew,double * deltaHydrotorquePrime)173 void FixNVEAsphereBase::implicitRotationUpdate
174 (
175 double deltaT, double* inertia,
176 double *angMom, double *torque, double* KslRot,
177 double *omegaNew,
178 double *deltaHydrotorquePrime
179 )
180 {
181 int index1_[3], index2_[3];
182 index1_[0]=2;index1_[1]=0;index1_[2]=1;
183 index2_[0]=1;index2_[1]=2;index2_[2]=0;
184 double omegaOld[3], dtfm, KslMDeltaT;
185 double omegaAngMomTerm;
186 for(int dirI=0;dirI<3;dirI++)
187 {
188 if (inertia[dirI] == 0.0)
189 {
190 omegaOld[dirI] = 0.0;
191 continue;
192 }
193 omegaOld[dirI] = angMom[dirI] / inertia[dirI];
194 }
195 for(int dirI=0;dirI<3;dirI++)
196 {
197 omegaAngMomTerm = omegaOld[index2_[dirI]]*angMom[index1_[dirI]]
198 - omegaOld[index1_[dirI]]*angMom[index2_[dirI]];
199
200 dtfm = deltaT
201 / ( inertia[dirI]);
202 KslMDeltaT = KslRot[dirI] * dtfm;
203
204 //calculate new rotation velocity
205 omegaNew[dirI] = ( omegaOld[dirI]
206 + dtfm
207 * ( torque[dirI] //this is the TOTAL torque!
208 - omegaAngMomTerm
209 + KslRot[dirI]*omegaOld[dirI]
210 )
211 )
212 /
213 ( 1.0+KslMDeltaT);
214
215 //save increment in hdtorque to update torque since particle is now rotation with omegaNew
216 deltaHydrotorquePrime[dirI] = KslRot[dirI]*(omegaOld[dirI]-omegaNew[dirI]);
217
218 //update the angular momentum
219 angMom[dirI] = omegaNew[dirI] * inertia[dirI]; //update velocity for a half step!
220 }
221 }
222
223 /* ---------------------------------------------------------------------- */
224
rotationUpdate(bool updateQuaternion)225 void FixNVEAsphereBase::rotationUpdate(bool updateQuaternion)
226 {
227 double omegaNew[3];
228 double exone[3],eyone[3],ezone[3];
229 double angMomPrime[3];
230 double torquePrime[3];
231 double omegaNewPrime[3];
232 double deltaHydrotorquePrime[3];
233 double deltaHydrotorque[3];
234 double rot[3][3];
235
236 // This is okey only for SPHERICAL particles
237 // update 1/2 step for omega
238
239 double **angmom = atom->angmom;
240 double **quat = atom->quaternion;
241 double **inertia = atom->inertia;
242 double **omega = atom->omega;
243 double **torque = atom->torque;
244 int *mask = atom->mask;
245 int nlocal = atom->nlocal;
246
247 double dtq = 0.50 * dtv;
248
249 //save rotation rate to array if necessary
250 orientation = NULL;
251 ksl_rotation = NULL;
252 hdtorque = NULL;
253
254 if(couple_fix_id > -1) {
255 ksl_rotation = ((FixCfdCouplingForceImplicit*)modify->fix[couple_fix_id])->fix_KslRotation_->array_atom;
256 hdtorque = ((FixCfdCouplingForceImplicit*)modify->fix[couple_fix_id])->fix_hdtorque_->array_atom;
257 orientation = ((FixCfdCouplingForceImplicit*)modify->fix[couple_fix_id])->fix_ex_->array_atom;
258 }
259
260 for (int i = 0; i < nlocal; i++) {
261 if (mask[i] & groupbit)
262 {
263 MathExtra::quat_to_mat(quat[i],rot); //rotation matrix
264
265 //Compute all relevant quantities in the particle's coordinate system called 'prime'
266 MathExtra::transpose_matvec(rot,angmom[i],angMomPrime);
267 MathExtra::transpose_matvec(rot,torque[i],torquePrime);
268
269 /* printf("\n***rotationUpdate(): angmom/prime: %g %g %g/%g %g %g, torque/prime: %g %g %g/%g %g %g \n\n",
270 angmom[i][0],angmom[i][1],angmom[i][2],
271 angMomPrime[0],angMomPrime[1],angMomPrime[2],
272 torque[i][0],torque[i][1],torque[i][2],
273 torquePrime[0],torquePrime[1],torquePrime[2]
274 );
275
276 printf("\n***rotationUpdate(): KslRotation: %g %g %g\n\n",
277 KslRotation[i][0],KslRotation[i][1],KslRotation[i][2]
278 );
279 */
280 //Implicit angmom and omegaNew update in the 'prime' coordinate system
281 implicitRotationUpdate
282 (
283 dtf, inertia[i],
284 angMomPrime, torquePrime, ksl_rotation[i],
285 omegaNewPrime,
286 deltaHydrotorquePrime
287 );
288
289 //Transform result to global coordinate system and update hydro torque (since it is the total!)
290 MathExtra::matvec(rot,angMomPrime,angmom[i]);
291 MathExtra::matvec(rot,omegaNewPrime,omegaNew);
292 MathExtra::matvec(rot,deltaHydrotorquePrime,deltaHydrotorque); //save torque
293 for(int dirI=0; dirI<3; dirI++)
294 hdtorque[i][dirI] += deltaHydrotorque[dirI];
295
296 // compute omega at 1/2 step from angmom at 1/2 step and current q
297 // update quaternion a full step via Richardson iteration
298 // returns new normalized quaternion
299 if(updateQuaternion)
300 MathExtra::richardson(quat[i],angmom[i],omegaNew,inertia[i],dtq);
301
302 omega[i][0]=omegaNew[0];
303 omega[i][1]=omegaNew[1];
304 omega[i][2]=omegaNew[2];
305
306 if(ksl_rotation && updateQuaternion)
307 {
308 MathExtra::q_to_exyz(quat[i],exone,eyone,ezone);
309 orientation[i][0] = exone[0];
310 orientation[i][1] = exone[1];
311 orientation[i][2] = exone[2];
312 }
313 } //end of particle
314 }
315 }
316
317 /* ---------------------------------------------------------------------- */
318
initial_integrate(int vflag)319 void FixNVEAsphereBase::initial_integrate(int vflag)
320 {
321 double **x = atom->x;
322 double **v = atom->v;
323 double **f = atom->f;
324 double **omega = atom->omega;
325 double **torque = atom->torque;
326 double *rmass = atom->rmass;
327 double **angmom = atom->angmom;
328 double **quat = atom->quaternion;
329 double **inertia = atom->inertia;
330 int *mask = atom->mask;
331 int nlocal = atom->nlocal;
332 if (igroup == atom->firstgroup) nlocal = atom->nfirst;
333
334 orientation = NULL;
335 ksl_rotation = NULL;
336 hdtorque = NULL;
337
338 if(couple_fix_id > -1) {
339 ksl_rotation = ((FixCfdCouplingForceImplicit*)modify->fix[couple_fix_id])->fix_KslRotation_->array_atom;
340 hdtorque = ((FixCfdCouplingForceImplicit*)modify->fix[couple_fix_id])->fix_hdtorque_->array_atom;
341 orientation = ((FixCfdCouplingForceImplicit*)modify->fix[couple_fix_id])->fix_ex_->array_atom;
342 }
343
344 double tbody[3], rotation_matrix[9];
345 double fquat[4], mbody[3], wbody[3], conjqm[4];
346
347 double angMomPrime[3];
348 double omegaNewPrime[3], omegaOldPrime[3];
349 const double dtf2 = 2.0 * dtf;
350
351 const double dtq = 0.5 * dtv;
352 for (int i = 0; i < nlocal; i++)
353 {
354 if (mask[i] & groupbit)
355 {
356 #ifdef LIGGGHTS_DEBUG
357 if(std::isnan(LAMMPS_NS::vectorMag3D(x[i])))
358 error->fix_error(FLERR,this,"x[i] is NaN!");
359 if(std::isnan(LAMMPS_NS::vectorMag3D(v[i])))
360 error->fix_error(FLERR,this,"v[i] is NaN!");
361 if(std::isnan(LAMMPS_NS::vectorMag3D(f[i])))
362 error->fix_error(FLERR,this,"f[i] is NaN!");
363 if(std::isnan(LAMMPS_NS::vectorMag4D(quat[i])))
364 error->fix_error(FLERR,this,"quat[i] is NaN!");
365 if(std::isnan(LAMMPS_NS::vectorMag3D(angmom[i])))
366 error->fix_error(FLERR,this,"angmom[i] is NaN!");
367 if(std::isnan(LAMMPS_NS::vectorMag3D(omega[i])))
368 error->fix_error(FLERR,this,"omega[i] is NaN!");
369 if(std::isnan(LAMMPS_NS::vectorMag3D(torque[i])))
370 error->fix_error(FLERR,this,"torque[i] is NaN!");
371 #endif
372 MathExtraLiggghtsNonspherical::quat_to_mat(quat[i], rotation_matrix);
373 const double dtfm = dtf / rmass[i];
374 MathExtraLiggghtsNonspherical::transpose_matvec(rotation_matrix, angmom[i], angMomPrime);
375 omegaOldPrime[0] = angMomPrime[0] / inertia[i][0];
376 omegaOldPrime[1] = angMomPrime[1] / inertia[i][1];
377 omegaOldPrime[2] = angMomPrime[2] / inertia[i][2];
378
379 //update velocity by step t+1/2dt
380
381 v[i][0] += dtfm * f[i][0];
382 v[i][1] += dtfm * f[i][1];
383 v[i][2] += dtfm * f[i][2];
384
385 #ifdef LIGGGHTS_DEBUG
386 if(std::isnan(LAMMPS_NS::vectorMag3D(v[i])))
387 error->fix_error(FLERR,this,"v[i] is NaN!");
388 #endif
389 //update position by step t+dt
390 x[i][0] += dtv * v[i][0];
391 x[i][1] += dtv * v[i][1];
392 x[i][2] += dtv * v[i][2];
393 #ifdef LIGGGHTS_DEBUG
394 if(std::isnan(LAMMPS_NS::vectorMag3D(x[i])))
395 error->fix_error(FLERR,this,"x[i] is NaN!");
396 #endif
397
398 if(integration_scheme == 0)
399 {
400 //update angular moment by step t+1/2dt
401 angmom[i][0] += dtf * torque[i][0];
402 angmom[i][1] += dtf * torque[i][1];
403 angmom[i][2] += dtf * torque[i][2];
404 //update angular velocity by step t+1/2dt
405 MathExtra::mq_to_omega(angmom[i],quat[i],inertia[i],omega[i]);
406 MathExtraLiggghtsNonspherical::transpose_matvec(rotation_matrix, omega[i], omegaNewPrime);
407 update_hdtorque(i, rotation_matrix, omegaOldPrime, omegaNewPrime);
408 //update quaternion by step t+dt
409 MathExtra::richardson(quat[i],angmom[i],omega[i],inertia[i],dtq);
410 }
411 //symplectic scheme
412 else if(integration_scheme == 1)
413 {
414 MathExtraLiggghtsNonspherical::transpose_matvec(rotation_matrix, angmom[i],mbody);
415 MathExtraLiggghtsNonspherical::calc_conjqm(quat[i],mbody,conjqm);
416
417 MathExtraLiggghtsNonspherical::transpose_matvec(rotation_matrix, torque[i], tbody);
418 MathExtra::quatvec(quat[i], tbody, fquat);
419
420 conjqm[0] += dtf2 * fquat[0];
421 conjqm[1] += dtf2 * fquat[1];
422 conjqm[2] += dtf2 * fquat[2];
423 conjqm[3] += dtf2 * fquat[3];
424
425 MathExtraLiggghtsNonspherical::no_squish_rotate(3,conjqm,quat[i],inertia[i],dtq);
426 MathExtraLiggghtsNonspherical::no_squish_rotate(2,conjqm,quat[i],inertia[i],dtq);
427 MathExtraLiggghtsNonspherical::no_squish_rotate(1,conjqm,quat[i],inertia[i],dtv);
428 MathExtraLiggghtsNonspherical::no_squish_rotate(2,conjqm,quat[i],inertia[i],dtq);
429 MathExtraLiggghtsNonspherical::no_squish_rotate(3,conjqm,quat[i],inertia[i],dtq);
430
431 MathExtra::qnormalize(quat[i]);
432 MathExtraLiggghtsNonspherical::quat_to_mat(quat[i], rotation_matrix);
433
434 MathExtra::invquatvec(quat[i],conjqm,mbody);
435 mbody[0] *= 0.5;
436 mbody[1] *= 0.5;
437 mbody[2] *= 0.5;
438
439 wbody[0] = mbody[0] / inertia[i][0];
440 wbody[1] = mbody[1] / inertia[i][1];
441 wbody[2] = mbody[2] / inertia[i][2];
442
443 MathExtraLiggghtsNonspherical::matvec(rotation_matrix, mbody, angmom[i]);
444 MathExtraLiggghtsNonspherical::matvec(rotation_matrix, wbody, omega[i]);
445 } else if (integration_scheme == 2) { //direct integration, 2nd order predictor-corrector
446 double omega_half[3];
447 MathExtraLiggghtsNonspherical::quat_to_mat(quat[i], rotation_matrix); //calculate rotation matrix from quaternion
448 MathExtraLiggghtsNonspherical::transpose_matvec(rotation_matrix, omega[i],wbody); //angular velocity in body principal axes
449 MathExtraLiggghtsNonspherical::transpose_matvec(rotation_matrix, torque[i], tbody); //torque in body principal axes
450 omega_half[0] = wbody[0];
451 omega_half[1] = wbody[1];
452 omega_half[2] = wbody[2];
453 integrate_dynamic_euler(dtf, wbody, tbody, inertia[i]); //calculate angular velocity at step t+dt
454 omega_half[0] = 0.5*(wbody[0] + omega_half[0]); //angular velocity at step t+dt/2
455 omega_half[1] = 0.5*(wbody[1] + omega_half[1]);
456 omega_half[2] = 0.5*(wbody[2] + omega_half[2]);
457
458 integrate_quaternion(dtq, quat[i], wbody); //calculate quaternion at step t+dt
459 update_hdtorque(i, rotation_matrix, omegaOldPrime, wbody);
460
461 mbody[0] = inertia[i][0]*wbody[0]; //calculate angular momentum at step t+dt from angular velocity in body's principal axes
462 mbody[1] = inertia[i][1]*wbody[1];
463 mbody[2] = inertia[i][2]*wbody[2];
464
465 MathExtraLiggghtsNonspherical::matvec(rotation_matrix, mbody, angmom[i]); // angular momentum to global frame
466 MathExtraLiggghtsNonspherical::matvec(rotation_matrix, wbody, omega[i]); // angular velocity to global frame
467
468 } else if(integration_scheme == 3) { //woodem scheme
469 double angmom_half[3], angmom_half_local[3], angmom_next_local[3];
470
471 for(int j = 0; j < 3; j++)
472 angmom_half[j] = angmom[i][j] + torque[i][j]*dtf;
473 MathExtraLiggghtsNonspherical::transpose_matvec(rotation_matrix, angmom_half, angmom_half_local);
474 for(int j = 0; j < 3; j++)
475 angmom[i][j] += torque[i][j]*dtf2;
476 MathExtraLiggghtsNonspherical::transpose_matvec(rotation_matrix, angmom[i], angmom_next_local);
477
478 double omega_half_local[3], omega_next_local[3], q_der[4], quat_half[4];
479 for(int j = 0; j < 3; j++) {
480 omega_half_local[j] = angmom_half_local[j] / inertia[i][j];
481 omega_next_local[j] = angmom_next_local[j] / inertia[i][j];
482 }
483 update_hdtorque(i, rotation_matrix, omegaOldPrime, omega_next_local);
484
485 MathExtra::quatvec(quat[i], omega_half_local, q_der);
486 for(int j = 0; j < 4; j++)
487 quat_half[j] = quat[i][j] + q_der[j]*dtq*0.5;
488 MathExtra::qnormalize(quat_half);
489
490 MathExtra::quatvec(quat_half, omega_next_local, q_der);
491 for(int j = 0; j < 4; j++)
492 quat[i][j] += q_der[j]*dtq;
493 MathExtra::qnormalize(quat[i]);
494
495 MathExtraLiggghtsNonspherical::quat_to_mat(quat[i], rotation_matrix);
496 MathExtraLiggghtsNonspherical::matvec(rotation_matrix, omega_next_local, omega[i]);
497 } else if(integration_scheme == 4) {
498 rotationUpdate(true);
499 }
500 else
501 error->one(FLERR,"Invalid integration scheme! Please choose between 0, 1 (default), 2 or 3!");
502 #ifdef LIGGGHTS_DEBUG
503 if(std::isnan(LAMMPS_NS::vectorMag4D(quat[i])))
504 error->fix_error(FLERR,this,"quat[i] is NaN!");
505 if(std::isnan(LAMMPS_NS::vectorMag3D(angmom[i])))
506 error->fix_error(FLERR,this,"angmom[i] is NaN!");
507 if(std::isnan(LAMMPS_NS::vectorMag3D(omega[i])))
508 error->fix_error(FLERR,this,"omega[i] is NaN!");
509 #endif
510 if(orientation) {
511 double exone[3], eyone[3], ezone[3];
512 MathExtra::q_to_exyz(quat[i],exone,eyone,ezone);
513 orientation[i][0] = exone[0];
514 orientation[i][1] = exone[1];
515 orientation[i][2] = exone[2];
516 }
517 }
518 }
519 }
520
521 /* ---------------------------------------------------------------------- */
522
final_integrate()523 void FixNVEAsphereBase::final_integrate()
524 {
525 double **v = atom->v;
526 double **f = atom->f;
527 double **omega = atom->omega;
528 double **torque = atom->torque;
529 double *rmass = atom->rmass;
530 double **angmom = atom->angmom;
531 double **quat = atom->quaternion;
532 double **inertia = atom->inertia;
533 int *mask = atom->mask;
534 int nlocal = atom->nlocal;
535 if (igroup == atom->firstgroup) nlocal = atom->nfirst;
536
537 double tbody[3], rotation_matrix[9];
538 double fquat[4], mbody[3], wbody[3];
539 double dtf2 = dtf * 2.0;
540 double conjqm[4];
541 double angMomPrime[3], omegaNewPrime[3], omegaOldPrime[3];
542
543 ksl_rotation = 0;
544 hdtorque = 0;
545 orientation = 0;
546 if(couple_fix_id > -1) {
547 ksl_rotation = ((FixCfdCouplingForceImplicit*)modify->fix[couple_fix_id])->fix_KslRotation_->array_atom;
548 hdtorque = ((FixCfdCouplingForceImplicit*)modify->fix[couple_fix_id])->fix_hdtorque_->array_atom;
549 orientation = ((FixCfdCouplingForceImplicit*)modify->fix[couple_fix_id])->fix_ex_->array_atom;
550 }
551
552 for (int i = 0; i < nlocal; i++)
553 if (mask[i] & groupbit) {
554 #ifdef LIGGGHTS_DEBUG
555 if(std::isnan(LAMMPS_NS::vectorMag3D(v[i])))
556 error->fix_error(FLERR,this,"v[i] is NaN!");
557 if(std::isnan(LAMMPS_NS::vectorMag3D(f[i])))
558 error->fix_error(FLERR,this,"f[i] is NaN!");
559 if(std::isnan(LAMMPS_NS::vectorMag4D(quat[i])))
560 error->fix_error(FLERR,this,"quat[i] is NaN!");
561 if(std::isnan(LAMMPS_NS::vectorMag3D(angmom[i])))
562 error->fix_error(FLERR,this,"angmom[i] is NaN!");
563 if(std::isnan(LAMMPS_NS::vectorMag3D(omega[i])))
564 error->fix_error(FLERR,this,"omega[i] is NaN!");
565 if(std::isnan(LAMMPS_NS::vectorMag3D(torque[i])))
566 error->fix_error(FLERR,this,"torque[i] is NaN!");
567 #endif
568 const double dtfm = dtf / rmass[i];
569 MathExtraLiggghtsNonspherical::quat_to_mat(quat[i], rotation_matrix);
570 MathExtraLiggghtsNonspherical::transpose_matvec(rotation_matrix, angmom[i], angMomPrime);
571 omegaOldPrime[0] = angMomPrime[0] / inertia[i][0];
572 omegaOldPrime[1] = angMomPrime[1] / inertia[i][1];
573 omegaOldPrime[2] = angMomPrime[2] / inertia[i][2];
574 if(hdtorque) {
575 vectorAdd3D(torque[i], hdtorque[i], torque[i]);
576 vectorAdd3D(torque[i], hdtorque[i], torque[i]);
577 }
578
579 //update velocity by step t+dt
580 v[i][0] += dtfm * f[i][0];
581 v[i][1] += dtfm * f[i][1];
582 v[i][2] += dtfm * f[i][2];
583
584 #ifdef LIGGGHTS_DEBUG
585 if(std::isnan(LAMMPS_NS::vectorMag3D(v[i])))
586 error->fix_error(FLERR,this,"v[i] is NaN!");
587 #endif
588
589 if(integration_scheme == 0) {
590 //update angular moment by step t+dt
591 angmom[i][0] += dtf * torque[i][0];
592 angmom[i][1] += dtf * torque[i][1];
593 angmom[i][2] += dtf * torque[i][2];
594 //update angular velocity by step t+dt
595 MathExtra::mq_to_omega(angmom[i],quat[i],inertia[i],omega[i]);
596 MathExtraLiggghtsNonspherical::transpose_matvec(rotation_matrix, omega[i], omegaNewPrime);
597 update_hdtorque(i, rotation_matrix, omegaOldPrime, omegaNewPrime);
598 } else if(integration_scheme == 1) {
599 MathExtraLiggghtsNonspherical::transpose_matvec(rotation_matrix, torque[i], tbody);
600 MathExtraLiggghtsNonspherical::transpose_matvec(rotation_matrix, angmom[i], mbody);
601 MathExtraLiggghtsNonspherical::calc_conjqm(quat[i], mbody, conjqm);
602
603 MathExtra::quatvec(quat[i], tbody, fquat);
604
605 conjqm[0] += dtf2 * fquat[0];
606 conjqm[1] += dtf2 * fquat[1];
607 conjqm[2] += dtf2 * fquat[2];
608 conjqm[3] += dtf2 * fquat[3];
609
610 MathExtra::invquatvec(quat[i], conjqm, mbody);
611 mbody[0] *= 0.5;
612 mbody[1] *= 0.5;
613 mbody[2] *= 0.5;
614
615 wbody[0] = mbody[0] / inertia[i][0];
616 wbody[1] = mbody[1] / inertia[i][1];
617 wbody[2] = mbody[2] / inertia[i][2];
618
619 MathExtraLiggghtsNonspherical::matvec(rotation_matrix, mbody, angmom[i]);
620 MathExtraLiggghtsNonspherical::matvec(rotation_matrix, wbody, omega[i]);
621 update_hdtorque(i, rotation_matrix, omegaOldPrime, wbody);
622 } else if (integration_scheme == 2) {
623 MathExtraLiggghtsNonspherical::quat_to_mat(quat[i], rotation_matrix); //calculate rotation matrix from quaternion
624 MathExtraLiggghtsNonspherical::transpose_matvec(rotation_matrix, omega[i],wbody); //angular velocity in body principal axes
625 MathExtraLiggghtsNonspherical::transpose_matvec(rotation_matrix, torque[i], tbody); //torque in body principal axes
626
627 integrate_dynamic_euler(dtf, wbody, tbody, inertia[i]); //calculate angular velocity at step t+dt
628
629 update_hdtorque(i, rotation_matrix, omegaOldPrime, wbody);
630
631 mbody[0] = inertia[i][0]*wbody[0]; //calculate angular momentum at step t+dt from angular velocity in body's principal axes
632 mbody[1] = inertia[i][1]*wbody[1];
633 mbody[2] = inertia[i][2]*wbody[2];
634
635 MathExtraLiggghtsNonspherical::matvec(rotation_matrix, mbody, angmom[i]); // angular momentum to global frame
636 MathExtraLiggghtsNonspherical::matvec(rotation_matrix, wbody, omega[i]); // angular velocity to global frame
637 } else if(integration_scheme == 4) {
638 rotationUpdate(false);
639 }
640
641 #ifdef LIGGGHTS_DEBUG
642 if(std::isnan(LAMMPS_NS::vectorMag4D(quat[i])))
643 error->fix_error(FLERR,this,"quat[i] is NaN!");
644 if(std::isnan(LAMMPS_NS::vectorMag3D(angmom[i])))
645 error->fix_error(FLERR,this,"angmom[i] is NaN!");
646 if(std::isnan(LAMMPS_NS::vectorMag3D(omega[i])))
647 error->fix_error(FLERR,this,"omega[i] is NaN!");
648 #endif
649 }
650 }
651
652