git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@5978 f3b2605a-c512-4ea7-a41b-209d697bcdaa

src/ASPHERE/fix_nh_asphere.cpp

@@ -57,83 +57,6 @@ void FixNHAsphere::init()
   FixNH::init();
 }
 
-/* ----------------------------------------------------------------------
-   Richardson iteration to update quaternion accurately
-------------------------------------------------------------------------- */
-
-void FixNHAsphere::richardson(double *q, double *m, double *moments)
-{
-  // compute omega at 1/2 step from m at 1/2 step and q at 0
-
-  double w[3];
-  omega_from_mq(q,m,moments,w);
-
-  // full update from dq/dt = 1/2 w q
-
-  double wq[4];
-  MathExtra::vecquat(w,q,wq);
-
-  double qfull[4];
-  qfull[0] = q[0] + dtq * wq[0];
-  qfull[1] = q[1] + dtq * wq[1];
-  qfull[2] = q[2] + dtq * wq[2];
-  qfull[3] = q[3] + dtq * wq[3];
-  MathExtra::qnormalize(qfull);
-
-  // 1st half of update from dq/dt = 1/2 w q
-
-  double qhalf[4];
-  qhalf[0] = q[0] + 0.5*dtq * wq[0];
-  qhalf[1] = q[1] + 0.5*dtq * wq[1];
-  qhalf[2] = q[2] + 0.5*dtq * wq[2];
-  qhalf[3] = q[3] + 0.5*dtq * wq[3];
-  MathExtra::qnormalize(qhalf);
-
-  // re-compute omega at 1/2 step from m at 1/2 step and q at 1/2 step
-  // recompute wq
-
-  omega_from_mq(qhalf,m,moments,w);
-  MathExtra::vecquat(w,qhalf,wq);
-
-  // 2nd half of update from dq/dt = 1/2 w q
-
-  qhalf[0] += 0.5*dtq * wq[0];
-  qhalf[1] += 0.5*dtq * wq[1];
-  qhalf[2] += 0.5*dtq * wq[2];
-  qhalf[3] += 0.5*dtq * wq[3];
-  MathExtra::qnormalize(qhalf);
-
-  // corrected Richardson update
-
-  q[0] = 2.0*qhalf[0] - qfull[0];
-  q[1] = 2.0*qhalf[1] - qfull[1];
-  q[2] = 2.0*qhalf[2] - qfull[2];
-  q[3] = 2.0*qhalf[3] - qfull[3];
-  MathExtra::qnormalize(q);
-}
-
-/* ----------------------------------------------------------------------
-   compute omega from angular momentum
-   w = omega = angular velocity in space frame
-   wbody = angular velocity in body frame
-   project space-frame angular momentum onto body axes
-     and divide by principal moments
-------------------------------------------------------------------------- */
-
-void FixNHAsphere::omega_from_mq(double *q, double *m, double *inertia,
-				  double *w)
-{
-  double rot[3][3];
-  MathExtra::quat_to_mat(q,rot);
-  
-  double wbody[3];
-  MathExtra::transpose_times_column3(rot,m,wbody);
-  wbody[0] /= inertia[0];
-  wbody[1] /= inertia[1];
-  wbody[2] /= inertia[2];
-  MathExtra::times_column3(rot,wbody,w);
-}
-
 /* ----------------------------------------------------------------------
    perform half-step update of angular momentum
 -----------------------------------------------------------------------*/
@@ -167,6 +90,8 @@ void FixNHAsphere::nve_v()
 
 void FixNHAsphere::nve_x()
 {
+  double omega[3];
+
   // standard nve_x position update
 
   FixNH::nve_x();
@@ -192,6 +117,9 @@ void FixNHAsphere::nve_x()
 
   for (int i = 0; i < nlocal; i++)
     if (mask[i] & groupbit) {
+
+      // principal moments of inertia
+
       shape = bonus[ellipsoid[i]].shape;
       quat = bonus[ellipsoid[i]].quat;
 
@@ -199,7 +127,12 @@ void FixNHAsphere::nve_x()
       inertia[1] = rmass[i] * (shape[0]*shape[0]+shape[2]*shape[2]) / 5.0;
       inertia[2] = rmass[i] * (shape[0]*shape[0]+shape[1]*shape[1]) / 5.0;
       
-      richardson(quat,angmom[i],inertia);
+      // compute omega at 1/2 step from angmom at 1/2 step and current q
+      // update quaternion a full step via Richardson iteration
+      // returns new normalized quaternion
+
+      MathExtra::mq_to_omega(angmom[i],quat,inertia,omega);
+      MathExtra::richardson(quat,angmom[i],omega,inertia,dtq);
     }
 }
 

src/ASPHERE/fix_nh_asphere.h

@@ -28,9 +28,6 @@ class FixNHAsphere : public FixNH {
   double dtq;
   class AtomVecEllipsoid *avec;
 
-  void richardson(double *, double *, double *);
-  void omega_from_mq(double *, double *, double *, double *);
-  void calculate_inertia();
   void nve_v();
   void nve_x();
   void nh_v_temp();

src/ASPHERE/fix_nve_asphere.cpp

@@ -63,7 +63,7 @@ void FixNVEAsphere::init()
 void FixNVEAsphere::initial_integrate(int vflag)
 {
   double dtfm;
-  double inertia[3];
+  double inertia[3],omega[3];
   double *shape,*quat;
 
   AtomVecEllipsoid::Bonus *bonus = avec->bonus;
@@ -93,8 +93,6 @@ void FixNVEAsphere::initial_integrate(int vflag)
       x[i][2] += dtv * v[i][2];
 
       // update angular momentum by 1/2 step
-      // update quaternion a full step via Richardson iteration
-      // returns new normalized quaternion
       
       angmom[i][0] += dtf * torque[i][0];
       angmom[i][1] += dtf * torque[i][1];
@@ -109,7 +107,12 @@ void FixNVEAsphere::initial_integrate(int vflag)
       inertia[1] = rmass[i] * (shape[0]*shape[0]+shape[2]*shape[2]) / 5.0;
       inertia[2] = rmass[i] * (shape[0]*shape[0]+shape[1]*shape[1]) / 5.0;
 
-      richardson(quat,angmom[i],inertia);
+      // compute omega at 1/2 step from angmom at 1/2 step and current q
+      // update quaternion a full step via Richardson iteration
+      // returns new normalized quaternion
+
+      MathExtra::mq_to_omega(angmom[i],quat,inertia,omega);
+      MathExtra::richardson(quat,angmom[i],omega,inertia,dtq);
     }
 }
 
@@ -140,80 +143,3 @@ void FixNVEAsphere::final_integrate()
       angmom[i][2] += dtf * torque[i][2];
     }
 }
-
-/* ----------------------------------------------------------------------
-   Richardson iteration to update quaternion accurately
-------------------------------------------------------------------------- */
-
-void FixNVEAsphere::richardson(double *q, double *m, double *moments)
-{
-  // compute omega at 1/2 step from m at 1/2 step and q at 0
-
-  double w[3];
-  omega_from_mq(m,q,moments,w);
-
-  // full update from dq/dt = 1/2 w q
-
-  double wq[4];
-  MathExtra::vecquat(w,q,wq);
-
-  double qfull[4];
-  qfull[0] = q[0] + dtq * wq[0];
-  qfull[1] = q[1] + dtq * wq[1];
-  qfull[2] = q[2] + dtq * wq[2];
-  qfull[3] = q[3] + dtq * wq[3];
-  MathExtra::qnormalize(qfull);
-
-  // 1st half of update from dq/dt = 1/2 w q
-
-  double qhalf[4];
-  qhalf[0] = q[0] + 0.5*dtq * wq[0];
-  qhalf[1] = q[1] + 0.5*dtq * wq[1];
-  qhalf[2] = q[2] + 0.5*dtq * wq[2];
-  qhalf[3] = q[3] + 0.5*dtq * wq[3];
-  MathExtra::qnormalize(qhalf);
-
-  // re-compute omega at 1/2 step from m at 1/2 step and q at 1/2 step
-  // recompute wq
-
-  omega_from_mq(m,qhalf,moments,w);
-  MathExtra::vecquat(w,qhalf,wq);
-
-  // 2nd half of update from dq/dt = 1/2 w q
-
-  qhalf[0] += 0.5*dtq * wq[0];
-  qhalf[1] += 0.5*dtq * wq[1];
-  qhalf[2] += 0.5*dtq * wq[2];
-  qhalf[3] += 0.5*dtq * wq[3];
-  MathExtra::qnormalize(qhalf);
-
-  // corrected Richardson update
-
-  q[0] = 2.0*qhalf[0] - qfull[0];
-  q[1] = 2.0*qhalf[1] - qfull[1];
-  q[2] = 2.0*qhalf[2] - qfull[2];
-  q[3] = 2.0*qhalf[3] - qfull[3];
-  MathExtra::qnormalize(q);
-}
-
-/* ----------------------------------------------------------------------
-   compute omega from angular momentum
-   w = omega = angular velocity in space frame
-   wbody = angular velocity in body frame
-   project space-frame angular momentum onto body axes
-     and divide by principal moments
-------------------------------------------------------------------------- */
-
-void FixNVEAsphere::omega_from_mq(double *m, double *q, double *moments,
-				  double *w)
-{
-  double rot[3][3];
-  MathExtra::quat_to_mat(q,rot);
-  
-  double wbody[3];
-  MathExtra::transpose_times_column3(rot,m,wbody);
-  wbody[0] /= moments[0];
-  wbody[1] /= moments[1];
-  wbody[2] /= moments[2];
-  MathExtra::times_column3(rot,wbody,w);
-}

src/ASPHERE/pair_gayberne.cpp

@@ -545,7 +545,7 @@ double PairGayBerne::gayberne_analytic(const int i,const int j,double a1[3][3],
   MathExtra::plus3(g1,g2,g12);
   double kappa[3];
   int ierror = MathExtra::mldivide3(g12,r12,kappa);
-  if (error) error->all("Bad matrix inversion in mldivide3");
+  if (ierror) error->all("Bad matrix inversion in mldivide3");
 
   // tempv = G12^-1*r12hat
 
@@ -576,7 +576,7 @@ double PairGayBerne::gayberne_analytic(const int i,const int j,double a1[3][3],
   double iota[3];
   MathExtra::plus3(b1,b2,b12);
   ierror = MathExtra::mldivide3(b12,r12,iota);
-  if (error) error->all("Bad matrix inversion in mldivide3");
+  if (ierror) error->all("Bad matrix inversion in mldivide3");
 
   // tempv = G12^-1*r12hat
 
@@ -726,7 +726,7 @@ double PairGayBerne::gayberne_lj(const int i,const int j,double a1[3][3],
   g12[1][2] = g1[1][2]; g12[2][1] = g1[2][1];
   double kappa[3];
   int ierror = MathExtra::mldivide3(g12,r12,kappa);
-  if (error) error->all("Bad matrix inversion in mldivide3");
+  if (ierror) error->all("Bad matrix inversion in mldivide3");
 
   // tempv = G12^-1*r12hat
 
@@ -762,7 +762,7 @@ double PairGayBerne::gayberne_lj(const int i,const int j,double a1[3][3],
   b12[0][2] = b1[0][2]; b12[2][0] = b1[2][0];
   b12[1][2] = b1[1][2]; b12[2][1] = b1[2][1];
   ierror = MathExtra::mldivide3(b12,r12,iota);
-  if (error) error->all("Bad matrix inversion in mldivide3");
+  if (ierror) error->all("Bad matrix inversion in mldivide3");
 
   // tempv = G12^-1*r12hat
 

src/ASPHERE/pair_resquared.cpp

@@ -612,7 +612,7 @@ double PairRESquared::resquared_analytic(const int i, const int j,
   double temp[3][3];
   MathExtra::plus3(wi.gamma,wj.gamma,temp);
   int ierror = MathExtra::mldivide3(temp,rhat,s);
-  if (error) error->all("Bad matrix inversion in mldivide3");
+  if (ierror) error->all("Bad matrix inversion in mldivide3");
 
   sigma12 = 1.0/sqrt(0.5*MathExtra::dot3(s,rhat));
   MathExtra::times_column3(wi.A,rhat,z1);
@@ -644,7 +644,7 @@ double PairRESquared::resquared_analytic(const int i, const int j,
   MathExtra::times3(wj.aTe,wj.A,temp2);
   MathExtra::plus3(temp,temp2,temp);
   ierror = MathExtra::mldivide3(temp,rhat,w);
-  if (error) error->all("Bad matrix inversion in mldivide3");
+  if (ierror) error->all("Bad matrix inversion in mldivide3");
 
   h12 = rnorm-sigma12;
   eta = lambda/nu;
@@ -898,7 +898,7 @@ double PairRESquared::resquared_lj(const int i, const int j,
   // energy
 
   int ierror = MathExtra::mldivide3(gamma,rhat,s);
-  if (error) error->all("Bad matrix inversion in mldivide3");
+  if (ierror) error->all("Bad matrix inversion in mldivide3");
 
   sigma12 = 1.0/sqrt(0.5*MathExtra::dot3(s,rhat));
   double temp[3][3];
@@ -907,7 +907,7 @@ double PairRESquared::resquared_lj(const int i, const int j,
   temp[1][1] += 1.0;
   temp[2][2] += 1.0;
   ierror = MathExtra::mldivide3(temp,rhat,w);
-  if (error) error->all("Bad matrix inversion in mldivide3");
+  if (ierror) error->all("Bad matrix inversion in mldivide3");
 
   h12 = rnorm-sigma12;
   chi = 2.0*MathExtra::dot3(rhat,w);

src/fix_rigid.cpp

@@ -1070,13 +1070,17 @@ void FixRigid::initial_integrate(int vflag)
     angmom[ibody][1] += dtf * torque[ibody][1] * tflag[ibody][1];
     angmom[ibody][2] += dtf * torque[ibody][2] * tflag[ibody][2];
 
-    // update quaternion a full step
+    // compute omega at 1/2 step from angmom at 1/2 step and current q
-    // returns new normalized quat
+    // update quaternion a full step via Richardson iteration
-    // returns ex_space,ey_space,ez_space for new quat
+    // returns new normalized quaternion, also updated omega at 1/2 step
-    // returns omega at 1/2 step (depends on angmom and quat)
+    // update ex,ey,ez to reflect new quaternion
 
-    richardson(quat[ibody],omega[ibody],angmom[ibody],inertia[ibody],
+    MathExtra::angmom_to_omega(angmom[ibody],ex_space[ibody],ey_space[ibody],
-	       ex_space[ibody],ey_space[ibody],ez_space[ibody]);
+			       ez_space[ibody],inertia[ibody],omega[ibody]);
+    MathExtra::richardson(quat[ibody],angmom[ibody],omega[ibody],
+			  inertia[ibody],dtq);
+    MathExtra::q_to_exyz(quat[ibody],
+			 ex_space[ibody],ey_space[ibody],ez_space[ibody]);
   }
 
   // virial setup before call to set_xv
@@ -1199,67 +1203,6 @@ void FixRigid::final_integrate()
   set_v();
 }
 
-/* ---------------------------------------------------------------------- */
-
-void FixRigid::richardson(double *q, double *w,
-			  double *m, double *moments,
-			  double *ex, double *ey, double *ez)
-{
-  // compute omega at 1/2 step from m at 1/2 step and q at 0
-  
-  MathExtra::angmom_to_omega(m,ex,ey,ez,moments,w);
-
-  // full update from dq/dt = 1/2 w q
-
-  double wq[4];
-  MathExtra::vecquat(w,q,wq);
-
-  double qfull[4];
-  qfull[0] = q[0] + dtq * wq[0];
-  qfull[1] = q[1] + dtq * wq[1];
-  qfull[2] = q[2] + dtq * wq[2];
-  qfull[3] = q[3] + dtq * wq[3];
-
-  MathExtra::qnormalize(qfull);
-
-  // 1st half update from dq/dt = 1/2 w q
-
-  double qhalf[4];
-  qhalf[0] = q[0] + 0.5*dtq * wq[0];
-  qhalf[1] = q[1] + 0.5*dtq * wq[1];
-  qhalf[2] = q[2] + 0.5*dtq * wq[2];
-  qhalf[3] = q[3] + 0.5*dtq * wq[3];
-
-  MathExtra::qnormalize(qhalf);
-
-  // udpate ex,ey,ez from qhalf
-  // re-compute omega at 1/2 step from m at 1/2 step and q at 1/2 step
-  // recompute wq
-
-  MathExtra::q_to_exyz(qhalf,ex,ey,ez);
-  MathExtra::angmom_to_omega(m,ex,ey,ez,moments,w);
-  MathExtra::vecquat(w,qhalf,wq);
-
-  // 2nd half update from dq/dt = 1/2 w q
-
-  qhalf[0] += 0.5*dtq * wq[0];
-  qhalf[1] += 0.5*dtq * wq[1];
-  qhalf[2] += 0.5*dtq * wq[2];
-  qhalf[3] += 0.5*dtq * wq[3];
-
-  MathExtra::qnormalize(qhalf);
-
-  // corrected Richardson update
-
-  q[0] = 2.0*qhalf[0] - qfull[0];
-  q[1] = 2.0*qhalf[1] - qfull[1];
-  q[2] = 2.0*qhalf[2] - qfull[2];
-  q[3] = 2.0*qhalf[3] - qfull[3];
-
-  MathExtra::qnormalize(q);
-  MathExtra::q_to_exyz(q,ex,ey,ez);
-}
-
 /* ----------------------------------------------------------------------
    apply evolution operators to quat, quat momentum
    see Miller paper cited in fix rigid/nvt and fix rigid/npt

src/fix_rigid.h

@@ -102,8 +102,6 @@ class FixRigid : public Fix {
   int ORIENT_DIPOLE,ORIENT_QUAT;
   int OMEGA,ANGMOM,TORQUE;
 
-  void richardson(double *, double *, double *, double *,
-		  double *, double *, double *);
   void no_squish_rotate(int, double *, double *, double *, double);
   void set_xv();
   void set_v();

src/math_extra.cpp

@@ -174,6 +174,58 @@ void rotate(double matrix[3][3], int i, int j, int k, int l,
   matrix[k][l] = h+s*(g-h*tau);
 }
 
+/* ----------------------------------------------------------------------
+   Richardson iteration to update quaternion from angular momentum
+   return new normalized quaternion q
+   also returns
+------------------------------------------------------------------------- */
+
+void richardson(double *q, double *m, double *w, double *moments, double dtq)
+{
+  // full update from dq/dt = 1/2 w q
+
+  double wq[4];
+  MathExtra::vecquat(w,q,wq);
+
+  double qfull[4];
+  qfull[0] = q[0] + dtq * wq[0];
+  qfull[1] = q[1] + dtq * wq[1];
+  qfull[2] = q[2] + dtq * wq[2];
+  qfull[3] = q[3] + dtq * wq[3];
+  MathExtra::qnormalize(qfull);
+
+  // 1st half update from dq/dt = 1/2 w q
+
+  double qhalf[4];
+  qhalf[0] = q[0] + 0.5*dtq * wq[0];
+  qhalf[1] = q[1] + 0.5*dtq * wq[1];
+  qhalf[2] = q[2] + 0.5*dtq * wq[2];
+  qhalf[3] = q[3] + 0.5*dtq * wq[3];
+  MathExtra::qnormalize(qhalf);
+
+  // re-compute omega at 1/2 step from m at 1/2 step and q at 1/2 step
+  // recompute wq
+
+  MathExtra::mq_to_omega(m,qhalf,moments,w);
+  MathExtra::vecquat(w,qhalf,wq);
+
+  // 2nd half update from dq/dt = 1/2 w q
+
+  qhalf[0] += 0.5*dtq * wq[0];
+  qhalf[1] += 0.5*dtq * wq[1];
+  qhalf[2] += 0.5*dtq * wq[2];
+  qhalf[3] += 0.5*dtq * wq[3];
+  MathExtra::qnormalize(qhalf);
+
+  // corrected Richardson update
+
+  q[0] = 2.0*qhalf[0] - qfull[0];
+  q[1] = 2.0*qhalf[1] - qfull[1];
+  q[2] = 2.0*qhalf[2] - qfull[2];
+  q[3] = 2.0*qhalf[3] - qfull[3];
+  MathExtra::qnormalize(q);
+}
+
 /* ----------------------------------------------------------------------
    compute rotation matrix from quaternion
    quat = [w i j k]
@@ -264,6 +316,30 @@ void angmom_to_omega(double *m, double *ex, double *ey, double *ez,
   w[2] = wbody[0]*ex[2] + wbody[1]*ey[2] + wbody[2]*ez[2];
 }
 
+/* ----------------------------------------------------------------------
+   compute omega from angular momentum
+   w = omega = angular velocity in space frame
+   wbody = angular velocity in body frame
+   project space-frame angular momentum onto body axes
+     and divide by principal moments
+------------------------------------------------------------------------- */
+
+void mq_to_omega(double *m, double *q, double *moments, double *w)
+{
+  double wbody[3];
+  double rot[3][3];
+
+  MathExtra::quat_to_mat(q,rot);
+  MathExtra::transpose_times_column3(rot,m,wbody);
+  if (moments[0] == 0.0) wbody[0] = 0.0;
+  else wbody[0] /= moments[0];
+  if (moments[1] == 0.0) wbody[1] = 0.0;
+  else wbody[1] /= moments[1];
+  if (moments[2] == 0.0) wbody[2] = 0.0;
+  else wbody[2] /= moments[2];
+  MathExtra::times_column3(rot,wbody,w);
+}
+
 /* ----------------------------------------------------------------------
    compute angular momentum from omega, both in space frame
    only know Idiag so need to do M = Iw in body frame

src/math_extra.h

@@ -68,6 +68,7 @@ namespace MathExtra {
   int jacobi(double matrix[3][3], double *evalues, double evectors[3][3]);
   void rotate(double matrix[3][3], int i, int j, int k, int l,
 	      double s, double tau);
+  void richardson(double *q, double *m, double *w, double *moments, double dtq);
 
   // shape matrix operations
   // upper-triangular 3x3 matrix stored in Voigt notation as 6-vector
@@ -94,6 +95,7 @@ namespace MathExtra {
 		       double *idiag, double *w);
   void omega_to_angmom(double *w, double *ex, double *ey, double *ez,
 		       double *idiag, double *m);
+  void mq_to_omega(double *m, double *q, double *moments, double *w);
   void exyz_to_q(double *ex, double *ey, double *ez, double *q);
   void q_to_exyz(double *q, double *ex, double *ey, double *ez);
   void quat_to_mat(const double *quat, double mat[3][3]);