reset quats for ellipsoids for general <-> restriced triclinic

src/atom_vec_ellipsoid.cpp

@@ -536,6 +536,83 @@ void AtomVecEllipsoid::write_data_bonus(FILE *fp, int n, double *buf, int /*flag
   }
 }
 
+/* ----------------------------------------------------------------------
+   convert read_data file info from general to restricted triclinic
+   parent class operates on data from Velocities section of data file
+   child class operates on ellipsoid quaternion
+------------------------------------------------------------------------- */
+
+void AtomVecEllipsoid::read_data_general_to_restricted(int nlocal_previous, int nlocal)
+{
+  int j;
+  
+  AtomVec::read_data_general_to_restricted(nlocal_previous, nlocal);
+
+  // quat_g2r = quat that rotates from general to restricted triclinic
+  // quat_new = ellipsoid quat converted to restricted triclinic
+  
+  double quat_g2r[4],quat_new[4];
+  MathExtra::mat_to_quat(domain->rotate_g2r,quat_g2r);
+  
+  for (int i = nlocal_previous; i < nlocal; i++) {
+    if (ellipsoid[i] < 0) continue;
+    j = ellipsoid[i];
+    MathExtra::quatquat(quat_g2r,bonus[j].quat,quat_new);
+    bonus[j].quat[0] = quat_new[0];
+    bonus[j].quat[1] = quat_new[1];
+    bonus[j].quat[2] = quat_new[2];
+    bonus[j].quat[3] = quat_new[3];
+  }
+}
+
+/* ----------------------------------------------------------------------
+   convert info output by write_data from restricted to general triclinic
+   parent class operates on x and data from Velocities section of data file
+   child class operates on ellipsoid quaternion
+------------------------------------------------------------------------- */
+
+void AtomVecEllipsoid::write_data_restricted_to_general()
+{
+  AtomVec::write_data_restricted_to_general();
+
+  memory->create(quat_hold,nlocal_bonus,4,"atomvec:quat_hold");
+
+  for (int i = 0; i < nlocal_bonus; i++)
+    memcpy(quat_hold[i],bonus[i].quat,4*sizeof(double));
+
+  // quat_r2g = quat that rotates from restricted to general triclinic
+  // quat_new = ellipsoid quat converted to general triclinic
+  
+  double quat_r2g[4],quat_new[4];
+  MathExtra::mat_to_quat(domain->rotate_r2g,quat_r2g);
+
+  for (int i = 0; i < nlocal_bonus; i++) {
+    MathExtra::quatquat(quat_r2g,bonus[i].quat,quat_new);
+    bonus[i].quat[0] = quat_new[0];
+    bonus[i].quat[1] = quat_new[1];
+    bonus[i].quat[2] = quat_new[2];
+    bonus[i].quat[3] = quat_new[3];
+  }
+}
+
+/* ----------------------------------------------------------------------
+   restore info output by write_data to restricted triclinic
+   original data is in "hold" arrays
+   parent class operates on x and data from Velocities section of data file
+   child class operates on ellipsoid quaternion
+------------------------------------------------------------------------- */
+
+void AtomVecEllipsoid::write_data_restore_restricted()
+{
+  AtomVec::write_data_restore_restricted();
+
+  for (int i = 0; i < nlocal_bonus; i++)
+    memcpy(bonus[i].quat,quat_hold[i],4*sizeof(double));
+
+  memory->destroy(quat_hold);
+  quat_hold = nullptr;
+}
+
 /* ----------------------------------------------------------------------
    set shape values in bonus data for particle I
    oriented aligned with xyz axes

src/atom_vec_ellipsoid.h

@@ -59,6 +59,10 @@ class AtomVecEllipsoid : public AtomVec {
   int pack_data_bonus(double *, int) override;
   void write_data_bonus(FILE *, int, double *, int) override;
 
+  void read_data_general_to_restricted(int, int);
+  void write_data_restricted_to_general();
+  void write_data_restore_restricted();
+
   // unique to AtomVecEllipsoid
 
   void set_shape(int, double, double, double);
@@ -69,6 +73,7 @@ class AtomVecEllipsoid : public AtomVec {
   int *ellipsoid;
   double *rmass;
   double **angmom;
+  double **quat_hold;
   
   int nghost_bonus, nmax_bonus;
   int ellipsoid_flag;

src/domain.cpp

@@ -307,7 +307,7 @@ void Domain::set_global_box()
     cprime[1] = yz;
     cprime[2] = boxhi[2] - boxlo[2];
 
-    // transform restricted A'B'C' to general triclinic A,B,C
+    // transform restricted A'B'C' to general triclinic ABC
     
     MathExtra::matvec(rotate_r2g,aprime,avec);
     MathExtra::matvec(rotate_r2g,bprime,bvec);

src/math_extra.cpp

@@ -362,6 +362,29 @@ void exyz_to_q(double *ex, double *ey, double *ez, double *q)
   qnormalize(q);
 }
 
+/* ----------------------------------------------------------------------
+   create unit quaternion from a rotation matrix
+   just a wrapper on exyz_to_q()
+   ex,ey,ez are columns of a rotation matrix
+------------------------------------------------------------------------- */
+
+void mat_to_quat(double mat[3][3], double *q)
+{
+  double ex[3],ey[3],ez[3];
+
+  ex[0] = mat[0][0];
+  ex[1] = mat[1][0];
+  ex[2] = mat[2][0];
+  ey[0] = mat[0][1];
+  ey[1] = mat[1][1];
+  ey[2] = mat[2][1];
+  ez[0] = mat[0][2];
+  ez[1] = mat[1][2];
+  ez[2] = mat[2][2];
+  
+  MathExtra::exyz_to_q(ex,ey,ez,q);
+}
+
 /* ----------------------------------------------------------------------
    compute space-frame ex,ey,ez from current quaternion q
    ex,ey,ez = space-frame coords of 1st,2nd,3rd principal axis
@@ -417,6 +440,7 @@ void quat_to_mat(const double *quat, double mat[3][3])
 /* ----------------------------------------------------------------------
    compute rotation matrix from quaternion conjugate
    quat = [w i j k]
+   similar logic to quat_to_mat()
 ------------------------------------------------------------------------- */
 
 void quat_to_mat_trans(const double *quat, double mat[3][3])
@@ -647,5 +671,4 @@ void tribbox(double *h, double radius, double *dist)
 
 /* ---------------------------------------------------------------------- */
 
-
 }

src/math_extra.h

@@ -100,6 +100,7 @@ void angmom_to_omega(double *m, double *ex, double *ey, double *ez, double *idia
 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 mat_to_quat(double mat[3][3], double *quat);
 void q_to_exyz(double *q, double *ex, double *ey, double *ez);
 void quat_to_mat(const double *quat, double mat[3][3]);
 void quat_to_mat_trans(const double *quat, double mat[3][3]);