upgrades to read/write data commands

src/atom_vec_body.cpp

@@ -586,11 +586,11 @@ void AtomVecBody::write_data_restricted_to_general()
 {
   AtomVec::write_data_restricted_to_general();
 
-  double quat[4],quat_r2g[4];
+  double quat[4];
   double *bquat;
-  
+  double *quat_r2g = domain->quat_r2g;
+    
   memory->create(quat_hold,nlocal_bonus,4,"atomvec:quat_hold");
-  MathExtra::qconjugate(domain->quat_g2r,quat_r2g);
   
   for (int i = 0; i < nlocal_bonus; i++) {
     bquat = bonus[i].quat;

src/atom_vec_ellipsoid.cpp

@@ -495,11 +495,11 @@ void AtomVecEllipsoid::write_data_restricted_to_general()
 {
   AtomVec::write_data_restricted_to_general();
 
-  double quat[4],quat_r2g[4];
+  double quat[4];
   double *bquat;
-  
+  double *quat_r2g = domain->quat_r2g;
+
   memory->create(quat_hold,nlocal_bonus,4,"atomvec:quat_hold");
-  MathExtra::qconjugate(domain->quat_g2r,quat_r2g);
   
   for (int i = 0; i < nlocal_bonus; i++) {
     bquat = bonus[i].quat;

src/atom_vec_tri.cpp

@@ -694,11 +694,11 @@ void AtomVecTri::write_data_restricted_to_general()
 {
   AtomVec::write_data_restricted_to_general();
 
-  double quat[4],quat_r2g[4];
+  double quat[4];
   double *bquat;
+  double *quat_r2g = domain->quat_r2g;
   
   memory->create(quat_hold,nlocal_bonus,4,"atomvec:quat_hold");
-  MathExtra::qconjugate(domain->quat_g2r,quat_r2g);
   
   for (int i = 0; i < nlocal_bonus; i++) {
     bquat = bonus[i].quat;

src/domain.cpp

@@ -580,13 +580,18 @@ void Domain::define_general_triclinic(double *avec_caller, double *bvec_caller,
     error->all(FLERR,"General triclinic box edge vector C invalid for 2d system");
 
   // error check for co-planar A,B,C
+  // if A x B is in C direction, only a rotation to restricted tri is needed
+  //   since restricted tri box will obey right-hand rule
+  // if not, an inversion (flip) of C vector is also needed
 
   double abcross[3];
   MathExtra::cross3(avec,bvec,abcross);
   double dot = MathExtra::dot3(abcross,cvec);
   if (dot == 0.0)
     error->all(FLERR,"General triclinic box edge vectors are co-planar");
-
+  if (dot > 0.0) triclinic_general_flip = 0;
+  else triclinic_general_flip = 1;
+  
   // quat1 = convert A into A' along +x-axis
   // rot1 = unit vector to rotate A around
   // theta1 = angle of rotation calculated from
@@ -627,18 +632,19 @@ void Domain::define_general_triclinic(double *avec_caller, double *bvec_caller,
   if (bvec1[2] > 0.0) theta2 = -theta2;
   MathExtra::axisangle_to_quat(xaxis,theta2,quat2);
 
-  // quat_g2r = transformation via single quat = quat2 * quat1
+  // quat_g2r = rotation via single quat = quat2 * quat1
-  // rotate_g2r = general to restricted transformation matrix
+  // quat_r2g = rotation from restricted to general
-  // if dot < 0.0 (A x B not in C direction)
+  // rotate_g2r = general to restricted rotation matrix
-  //   flip sign of z component of transform,
+  //   include flip of C vector if needed to obey right-hand rule
-  //   by flipping sign of 3rd row of rotate_g2r matrix
+  // rotate_r2g = restricted to general rotation matrix
-  // rotate_r2g = restricted to general transformation matrix
   //   simply a transpose of rotate_g2r since orthonormal
 
   MathExtra::quatquat(quat2,quat1,quat_g2r);
+  MathExtra::qconjugate(quat_g2r,quat_r2g);
+
   MathExtra::quat_to_mat(quat_g2r,rotate_g2r);
 
-  if (dot < 0.0) {
+  if (triclinic_general_flip) {
     rotate_g2r[2][0] = -rotate_g2r[2][0];
     rotate_g2r[2][1] = -rotate_g2r[2][1];
     rotate_g2r[2][2] = -rotate_g2r[2][2];
@@ -646,7 +652,7 @@ void Domain::define_general_triclinic(double *avec_caller, double *bvec_caller,
   
   MathExtra::transpose3(rotate_g2r,rotate_r2g);
 
-  // transform general ABC to restricted triclinic A'B'C'
+  // rotate general ABC to restricted triclinic A'B'C'
 
   double aprime[3],bprime[3],cprime[3];
   MathExtra::matvec(rotate_g2r,avec,aprime);

src/domain.h

@@ -41,6 +41,7 @@ class Domain : protected Pointers {
 
   int triclinic;          // 0 = orthog box, 1 = triclinic (restricted or general)
   int triclinic_general;  // 1 if mapping to/from general triclinic is stored, 0 if not
+  int triclinic_general_flip;  // 1 if general tri rotation needs to invert C edge vector
   
   // orthogonal box
 
@@ -90,9 +91,9 @@ class Domain : protected Pointers {
 
   // general triclinic box
   // boxlo = lower left corner
-  
+
   double avec[3], bvec[3], cvec[3];  // ABC edge vectors of general triclinic box
-  double quat_g2r[4];                // quaternion for general --> restricted rotation
+  double quat_g2r[4], quat_r2g[4];   // quaternions for general <--> restricted rotations
   double rotate_g2r[3][3];           // rotation matrix from general --> restricted tri
   double rotate_r2g[3][3];           // rotation matrix from restricted --> general tri
 

src/dump_atom.cpp

@@ -334,13 +334,14 @@ void DumpAtom::header_binary_triclinic_general(bigint ndump)
 
   fwrite(&update->ntimestep,sizeof(bigint),1,fp);
   fwrite(&ndump,sizeof(bigint),1,fp);
-  int general_tri = 2;
+  int triclinic_general_flag = 2;
-  fwrite(&general_tri,sizeof(int),1,fp);
+  fwrite(&triclinic_general_flag,sizeof(int),1,fp);
   fwrite(&domain->boundary[0][0],6*sizeof(int),1,fp);
   fwrite(domain->avec,3*sizeof(double),1,fp);
   fwrite(domain->bvec,3*sizeof(double),1,fp);
   fwrite(domain->cvec,3*sizeof(double),1,fp);
   fwrite(domain->boxlo,3*sizeof(double),1,fp);
+  fwrite(&size_one,sizeof(int),1,fp);
 
   header_unit_style_binary();
   header_time_binary();

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,ROTATE_G2R,ROTATE_R2G};
+     TRICLINIC_GENERAL,TRICLINIC_GENERAL_FLIP,QUAT_G2R};
 
 #define LB_FACTOR 1.1
 

src/read_restart.cpp

@@ -28,6 +28,7 @@
 #include "improper.h"
 #include "irregular.h"
 #include "label_map.h"
+#include "math_extra.h"
 #include "memory.h"
 #include "modify.h"
 #include "pair.h"
@@ -136,6 +137,21 @@ 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();
@@ -785,10 +801,10 @@ void ReadRestart::header()
 
     } else if (flag == TRICLINIC_GENERAL) {
       domain->triclinic_general = read_int();
-    } else if (flag == ROTATE_G2R) {
+    } else if (flag == TRICLINIC_GENERAL_FLIP) {
-      read_double_vec(9,&domain->rotate_g2r[0][0]);
+      domain->triclinic_general_flip = read_int();
-    } else if (flag == ROTATE_R2G) {
+    } else if (flag == QUAT_G2R) {
-      read_double_vec(9,&domain->rotate_r2g[0][0]);
+      read_double_vec(4,domain->quat_g2r);
 
     } else if (flag == SPECIAL_LJ) {
       read_int();

src/write_restart.cpp

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