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