more work on restart support

src/atom.cpp

@@ -1190,7 +1190,7 @@ void Atom::data_atoms(int n, char *buf, tagint id_offset, tagint mol_offset,
         xdata[2] = 0.0;
       }
       
-      // convert atom coords from general triclinic to restricted triclinic
+      // convert atom coords from general to restricted triclinic
       // so can decide which proc owns the atom
       
       if (triclinic_general) domain->general_to_restricted_coords(xdata);

src/atom_vec_line.cpp

@@ -352,6 +352,22 @@ void AtomVecLine::data_atom_bonus(int m, const std::vector<std::string> &values)
   double y1 = utils::numeric(FLERR, values[ivalue++], true, lmp);
   double x2 = utils::numeric(FLERR, values[ivalue++], true, lmp);
   double y2 = utils::numeric(FLERR, values[ivalue++], true, lmp);
+
+  // convert x1/y1 and x2/y2 from general to restricted triclniic
+  // x is already restricted triclinic
+  
+  if (domain->triclinic_general) {
+    double coords[3];
+    coords[0] = x1; coords[1] = y1; coords[2] = 0.0;
+    domain->general_to_restricted_coords(coords);
+    x1 = coords[0]; y1 = coords[1];
+    coords[0] = x2; coords[1] = y2; coords[2] = 0.0;
+    domain->general_to_restricted_coords(coords);
+    x2 = coords[0]; y2 = coords[1];
+  }
+
+  // calculate length and theta
+  
   double dx = x2 - x1;
   double dy = y2 - y1;
   double length = sqrt(dx * dx + dy * dy);
@@ -437,29 +453,6 @@ void AtomVecLine::data_atom_post(int ilocal)
   omega[ilocal][2] = 0.0;
 }
 
-/* ----------------------------------------------------------------------
-   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 bonus theta
-------------------------------------------------------------------------- */
-
-void AtomVecLine::read_data_general_to_restricted(int nlocal_previous, int nlocal)
-{
-  AtomVec::read_data_general_to_restricted(nlocal_previous, nlocal);
-
-  double btheta;
-  double theta_g2r = 2.0*acos(domain->quat_g2r[0]);
-  
-  for (int i = nlocal_previous; i < nlocal; i++) {
-    if (line[i] < 0) continue;
-    btheta = bonus[line[i]].theta;
-    btheta += theta_g2r;
-    if (btheta > MathConst::MY_PI) btheta -= MathConst::MY_2PI;
-    else if (btheta <= -MathConst::MY_PI) btheta += MathConst::MY_2PI;
-    bonus[line[i]].theta = btheta;
-  }
-}
-
 /* ----------------------------------------------------------------------
    modify values for AtomVec::pack_data() to pack
 ------------------------------------------------------------------------- */
@@ -501,8 +494,14 @@ int AtomVecLine::pack_data_bonus(double *buf, int /*flag*/)
   int i, j;
   double length, theta;
   double xc, yc, x1, x2, y1, y2;
+  double coords[3];
+
+  int triclinic_general = domain->triclinic_general;
+
+  double **x_bonus;
+  if (triclinic_general) x_bonus = x_hold;
+  else x_bonus = x;
 
-  double **x = atom->x;
   tagint *tag = atom->tag;
   int nlocal = atom->nlocal;
 
@@ -514,8 +513,9 @@ int AtomVecLine::pack_data_bonus(double *buf, int /*flag*/)
       j = line[i];
       length = bonus[j].length;
       theta = bonus[j].theta;
-      xc = x[i][0];
+      
-      yc = x[i][1];
+      xc = x_bonus[i][0];
+      yc = x_bonus[i][1];
       x1 = xc - 0.5 * cos(theta) * length;
       y1 = yc - 0.5 * sin(theta) * length;
       x2 = xc + 0.5 * cos(theta) * length;
@@ -524,6 +524,20 @@ int AtomVecLine::pack_data_bonus(double *buf, int /*flag*/)
       buf[m++] = y1;
       buf[m++] = x2;
       buf[m++] = y2;
+
+      // if triclinic_general:
+      // rotate 4 buf values from restricted to general triclinic
+      // output by write_data_bonus() as x1/y1 and x2/y2
+      
+      if (triclinic_general) {
+        coords[0] = buf[m-4]; coords[1] = buf[m-3]; coords[2] = 0.0;
+        domain->restricted_to_general_coords(coords);
+        buf[m-4] = coords[0]; buf[m-3] = coords[1];
+        coords[0] = buf[m-2]; coords[1] = buf[m-1]; coords[2] = 0.0;
+        domain->restricted_to_general_coords(coords);
+        buf[m-2] = coords[0]; buf[m-1] = coords[1];
+      }
+
     } else
       m += size_data_bonus;
   }

src/atom_vec_line.h

@@ -53,7 +53,6 @@ class AtomVecLine : public AtomVec {
 
   void create_atom_post(int) override;
   void data_atom_post(int) override;
-  void read_data_general_to_restricted(int, int);
   void pack_data_pre(int) override;
   void pack_data_post(int) override;
 

src/atom_vec_tri.cpp

@@ -510,8 +510,18 @@ void AtomVecTri::data_atom_bonus(int m, const std::vector<std::string> &values)
   MathExtra::sub3(c3, c1, c3mc1);
   double size = MAX(MathExtra::len3(c2mc1), MathExtra::len3(c3mc1));
 
+  // convert c1,c2,c3 from general to restricted triclniic
+  // x is already restricted triclinic
+  
+  if (domain->triclinic_general) {
+    domain->general_to_restricted_coords(c1);
+    domain->general_to_restricted_coords(c2);
+    domain->general_to_restricted_coords(c3);
+  }
+  
   // centroid = 1/3 of sum of vertices
-
+  // error if centroid is not within EPSILON of Atoms section coord
+  
   double centroid[3];
   centroid[0] = (c1[0] + c2[0] + c3[0]) / 3.0;
   centroid[1] = (c1[1] + c2[1] + c3[1]) / 3.0;
@@ -659,90 +669,6 @@ void AtomVecTri::data_atom_post(int ilocal)
   angmom[ilocal][2] = 0.0;
 }
 
-/* ----------------------------------------------------------------------
-   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 bonus quat
-------------------------------------------------------------------------- */
-
-void AtomVecTri::read_data_general_to_restricted(int nlocal_previous, int nlocal)
-{
-  AtomVec::read_data_general_to_restricted(nlocal_previous, nlocal);
-
-  double quat[4];
-  double *bquat;
-  
-  for (int i = nlocal_previous; i < nlocal; i++) {
-    if (tri[i] < 0) continue;
-    bquat = bonus[tri[i]].quat;
-    MathExtra::quatquat(domain->quat_g2r,bquat,quat);
-    bquat[0] = quat[0];
-    bquat[1] = quat[1];
-    bquat[2] = quat[2];
-    bquat[3] = quat[3];
-    MathExtra::qnormalize(bquat);
-  }
-}
-
-/* ----------------------------------------------------------------------
-   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 bonus quat
-------------------------------------------------------------------------- */
-
-void AtomVecTri::write_data_restricted_to_general()
-{
-  AtomVec::write_data_restricted_to_general();
-
-  double quat[4];
-  double *bquat;
-  double *quat_r2g = domain->quat_r2g;
-  
-  memory->create(quat_hold,nlocal_bonus,4,"atomvec:quat_hold");
-  
-  for (int i = 0; i < nlocal_bonus; i++) {
-    bquat = bonus[i].quat;
-    quat_hold[i][0] = bquat[0];
-    quat_hold[i][1] = bquat[1];
-    quat_hold[i][2] = bquat[2];
-    quat_hold[i][3] = bquat[3];
-
-    MathExtra::quatquat(quat_r2g,bquat,quat);
-    bquat[0] = quat[0];
-    bquat[1] = quat[1];
-    bquat[2] = quat[2];
-    bquat[3] = quat[3];
-    MathExtra::qnormalize(bquat);
-  }
-}
-
-/* ----------------------------------------------------------------------
-   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 bonus quat
-------------------------------------------------------------------------- */
-
-void AtomVecTri::write_data_restore_restricted()
-{
-  AtomVec::write_data_restore_restricted();
-
-  if (!quat_hold) return;
-
-  double *bquat;
-  
-  for (int i = 0; i < nlocal_bonus; i++) {
-    bquat = bonus[i].quat;
-    bquat[0] = quat_hold[i][0];
-    bquat[1] = quat_hold[i][1];
-    bquat[2] = quat_hold[i][2];
-    bquat[3] = quat_hold[i][3];
-  }
-  
-  memory->destroy(quat_hold);
-  quat_hold = nullptr;
-}
-
 /* ----------------------------------------------------------------------
    modify values for AtomVec::pack_data() to pack
 ------------------------------------------------------------------------- */
@@ -792,7 +718,12 @@ int AtomVecTri::pack_data_bonus(double *buf, int /*flag*/)
   double dc1[3], dc2[3], dc3[3];
   double p[3][3];
 
-  double **x = atom->x;
+  int triclinic_general = domain->triclinic_general;
+
+  double **x_bonus;
+  if (triclinic_general) x_bonus = x_hold;
+  else x_bonus = x;
+    
   tagint *tag = atom->tag;
   int nlocal = atom->nlocal;
 
@@ -806,9 +737,10 @@ int AtomVecTri::pack_data_bonus(double *buf, int /*flag*/)
       MathExtra::matvec(p, bonus[j].c1, dc1);
       MathExtra::matvec(p, bonus[j].c2, dc2);
       MathExtra::matvec(p, bonus[j].c3, dc3);
-      xc = x[i][0];
+
-      yc = x[i][1];
+      xc = x_bonus[i][0];
-      zc = x[i][2];
+      yc = x_bonus[i][1];
+      zc = x_bonus[i][2];
       buf[m++] = xc + dc1[0];
       buf[m++] = yc + dc1[1];
       buf[m++] = zc + dc1[2];
@@ -818,6 +750,17 @@ int AtomVecTri::pack_data_bonus(double *buf, int /*flag*/)
       buf[m++] = xc + dc3[0];
       buf[m++] = yc + dc3[1];
       buf[m++] = zc + dc3[2];
+
+      // if triclinic_general:
+      // rotate 9 buf values from restricted to general triclinic
+      // output by write_data_bonus() as c1,c2,c3
+      
+      if (triclinic_general) {
+        domain->restricted_to_general_coords(&buf[m-9]);
+        domain->restricted_to_general_coords(&buf[m-6]);
+        domain->restricted_to_general_coords(&buf[m-3]);
+      }
+
     } else
       m += size_data_bonus;
   }

src/atom_vec_tri.h

@@ -55,9 +55,6 @@ class AtomVecTri : public AtomVec {
 
   void create_atom_post(int) override;
   void data_atom_post(int) override;
-  void read_data_general_to_restricted(int, int);
-  void write_data_restricted_to_general();
-  void write_data_restore_restricted();
   void pack_data_pre(int) override;
   void pack_data_post(int) override;
 

src/lmprestart.h

@@ -38,7 +38,7 @@ enum{VERSION,SMALLINT,TAGINT,BIGINT,
      EXTRA_BOND_PER_ATOM,EXTRA_ANGLE_PER_ATOM,EXTRA_DIHEDRAL_PER_ATOM,
      EXTRA_IMPROPER_PER_ATOM,EXTRA_SPECIAL_PER_ATOM,ATOM_MAXSPECIAL,
      NELLIPSOIDS,NLINES,NTRIS,NBODIES,ATIME,ATIMESTEP,LABELMAP,
-     TRICLINIC_GENERAL,TRICLINIC_GENERAL_FLIP,QUAT_G2R};
+     TRICLINIC_GENERAL,TRICLINIC_GENERAL_FLIP,ROTATE_G2R};
 
 #define LB_FACTOR 1.1
 

src/read_restart.cpp

@@ -137,21 +137,6 @@ void ReadRestart::command(int narg, char **arg)
   atom->avec->grow(n);
   n = atom->nmax;
 
-  // setup simulation box
-  // for general triclinic, need to generate additional info which
-  //   Domain::define_general_triclinic() would have created
-
-  if (domain->triclinic_general) {
-    MathExtra::qconjugate(domain->quat_g2r,domain->quat_r2g);
-    MathExtra::quat_to_mat(domain->quat_g2r,domain->rotate_g2r);
-    if (domain->triclinic_general_flip) {
-      domain->rotate_g2r[2][0] = -domain->rotate_g2r[2][0];
-      domain->rotate_g2r[2][1] = -domain->rotate_g2r[2][1];
-      domain->rotate_g2r[2][2] = -domain->rotate_g2r[2][2];
-    }
-    MathExtra::transpose3(domain->rotate_g2r,domain->rotate_r2g);
-  }
-  
   domain->print_box("  ");
   domain->set_initial_box(0);
   domain->set_global_box();
@@ -803,8 +788,9 @@ void ReadRestart::header()
       domain->triclinic_general = read_int();
     } else if (flag == TRICLINIC_GENERAL_FLIP) {
       domain->triclinic_general_flip = read_int();
-    } else if (flag == QUAT_G2R) {
+    } else if (flag == ROTATE_G2R) {
-      read_double_vec(4,domain->quat_g2r);
+      read_double_vec(9,&domain->rotate_g2r[0][0]);
+      MathExtra::transpose3(domain->rotate_g2r,domain->rotate_r2g);
 
     } else if (flag == SPECIAL_LJ) {
       read_int();

src/write_restart.cpp

@@ -451,7 +451,7 @@ void WriteRestart::header()
   write_int(TRICLINIC_GENERAL,domain->triclinic_general);
   if (domain->triclinic_general) {
     write_int(TRICLINIC_GENERAL_FLIP,domain->triclinic_general_flip);
-    write_double_vec(QUAT_G2R,4,domain->quat_g2r);
+    write_double_vec(ROTATE_G2R,9,&domain->rotate_g2r[0][0]);
   }
   
   write_double_vec(SPECIAL_LJ,3,&force->special_lj[1]);