add general triclinic options to dump custom

src/atom.cpp

@@ -1181,7 +1181,7 @@ void Atom::data_atoms(int n, char *buf, tagint id_offset, tagint mol_offset,
 
       // convert atom coords from general triclinic to restricted triclinic
       
-      if (triclinic_general) domain->general_to_restricted(xdata);
+      if (triclinic_general) domain->general_to_restricted_coords(xdata);
 
       // apply shift if requested by read_data command
       

src/create_atoms.cpp

@@ -1236,31 +1236,31 @@ void CreateAtoms::add_lattice()
     double point[3];
 
     point[0] = bboxlo[0]; point[1] = bboxlo[1]; point[2] = bboxlo[2];
-    domain->restricted_to_general(point);
+    domain->restricted_to_general_coords(point);
     domain->lattice->bbox(1, point[0], point[1], point[2], xmin, ymin, zmin, xmax, ymax, zmax);
     point[0] = bboxhi[0]; point[1] = bboxlo[1]; point[2] = bboxlo[2];
-    domain->restricted_to_general(point);
+    domain->restricted_to_general_coords(point);
     domain->lattice->bbox(1, point[0], point[1], point[2], xmin, ymin, zmin, xmax, ymax, zmax);
 
     point[0] = bboxlo[0]; point[1] = bboxhi[1]; point[2] = bboxlo[2];
-    domain->restricted_to_general(point);
+    domain->restricted_to_general_coords(point);
     domain->lattice->bbox(1, point[0], point[1], point[2], xmin, ymin, zmin, xmax, ymax, zmax);
     point[0] = bboxhi[0]; point[1] = bboxhi[1]; point[2] = bboxlo[2];
-    domain->restricted_to_general(point);
+    domain->restricted_to_general_coords(point);
     domain->lattice->bbox(1, point[0], point[1], point[2], xmin, ymin, zmin, xmax, ymax, zmax);
 
     point[0] = bboxlo[0]; point[1] = bboxlo[1]; point[2] = bboxhi[2];
-    domain->restricted_to_general(point);
+    domain->restricted_to_general_coords(point);
     domain->lattice->bbox(1, point[0], point[1], point[2], xmin, ymin, zmin, xmax, ymax, zmax);
     point[0] = bboxhi[0]; point[1] = bboxlo[1]; point[2] = bboxhi[2];
-    domain->restricted_to_general(point);
+    domain->restricted_to_general_coords(point);
     domain->lattice->bbox(1, point[0], point[1], point[2], xmin, ymin, zmin, xmax, ymax, zmax);
 
     point[0] = bboxlo[0]; point[1] = bboxhi[1]; point[2] = bboxhi[2];
-    domain->restricted_to_general(point);
+    domain->restricted_to_general_coords(point);
     domain->lattice->bbox(1, point[0], point[1], point[2], xmin, ymin, zmin, xmax, ymax, zmax);
     point[0] = bboxhi[0]; point[1] = bboxhi[1]; point[2] = bboxhi[2];
-    domain->restricted_to_general(point);
+    domain->restricted_to_general_coords(point);
     domain->lattice->bbox(1, point[0], point[1], point[2], xmin, ymin, zmin, xmax, ymax, zmax);
   }
   
@@ -1373,7 +1373,7 @@ void CreateAtoms::loop_lattice(int action)
 
           // convert from general to restricted triclinic coords
 
-          if (triclinic_general) domain->general_to_restricted(x);
+          if (triclinic_general) domain->general_to_restricted_coords(x);
           
           // if a region was specified, test if atom is in it
 

src/create_box.cpp

@@ -96,7 +96,7 @@ void CreateBox::command(int narg, char **arg)
 
   // setup general triclinic box (with no region)
   // read next box extent arguments to create ABC edge vectors + origin
-  // setup_general_triclinic() converts
+  // define_general_triclinic() converts
   //   ABC edge vectors + origin to restricted triclinic
 
   } else if (triclinic_general) {
@@ -146,9 +146,9 @@ void CreateBox::command(int narg, char **arg)
     cvec[1] = py - origin[1];
     cvec[2] = pz - origin[2];
 
-    // setup general triclinic box within Domain
+    // define general triclinic box within Domain
     
-    domain->setup_general_triclinic(avec,bvec,cvec,origin);
+    domain->define_general_triclinic(avec,bvec,cvec,origin);
   }
 
   // if molecular, zero out topology info

src/domain.cpp

@@ -290,6 +290,29 @@ void Domain::set_global_box()
     boxhi_bound[1] = MAX(boxhi[1],boxhi[1]+yz);
     boxhi_bound[2] = boxhi[2];
   }
+
+  // update general triclinic box if defined
+  // reset ABC edge vectors from restricted triclinic box
+  // boxlo = lower left corner of general triclinic box
+  
+  if (triclinic_general) {
+    double aprime[3],bprime[3],cprime[3];
+
+    aprime[0] = boxhi[0] - boxlo[0];
+    aprime[1] = aprime[2] = 0.0;
+    bprime[0] = xy;
+    bprime[1] = boxhi[1] - boxlo[1];
+    bprime[2] = 0.0;
+    cprime[0] = xz;
+    cprime[1] = yz;
+    cprime[2] = boxhi[2] - boxlo[2];
+
+    // transform restricted A'B'C' to general triclinic A,B,C
+    
+    MathExtra::matvec(rotate_r2g,aprime,avec);
+    MathExtra::matvec(rotate_r2g,bprime,bvec);
+    MathExtra::matvec(rotate_r2g,cprime,cvec);
+  }
 }
 
 /* ----------------------------------------------------------------------
@@ -527,7 +550,7 @@ void Domain::reset_box()
    create rotation matrices for general <--> restricted transformations
 ------------------------------------------------------------------------- */
 
-void Domain::setup_general_triclinic(double *avec_caller, double *bvec_caller,
+void Domain::define_general_triclinic(double *avec_caller, double *bvec_caller,
                                       double *cvec_caller, double *origin_caller)
 {
   if (triclinic || triclinic_general)
@@ -547,9 +570,9 @@ void Domain::setup_general_triclinic(double *avec_caller, double *bvec_caller,
   cvec[1] = cvec_caller[1];
   cvec[2] = cvec_caller[2];
 
-  gtri_origin[0] = origin_caller[0];
+  boxlo[0] = origin_caller[0];
-  gtri_origin[1] = origin_caller[1];
+  boxlo[1] = origin_caller[1];
-  gtri_origin[2] = origin_caller[2];
+  boxlo[2] = origin_caller[2];
 
   // error check on cvec for 2d systems
 
@@ -624,7 +647,7 @@ void Domain::setup_general_triclinic(double *avec_caller, double *bvec_caller,
   
   MathExtra::transpose3(rotate_g2r,rotate_r2g);
 
-  // A',B',C' = transformation of A,B,C to restricted triclinic
+  // transform general ABC to restricted triclinic A'B'C'
 
   double aprime[3],bprime[3],cprime[3];
   MathExtra::matvec(rotate_g2r,avec,aprime);
@@ -633,10 +656,6 @@ void Domain::setup_general_triclinic(double *avec_caller, double *bvec_caller,
 
   // set restricted triclinic boxlo, boxhi, and tilt factors
 
-  boxlo[0] = gtri_origin[0];
-  boxlo[1] = gtri_origin[1];
-  boxlo[2] = gtri_origin[2];
-
   boxhi[0] = boxlo[0] + aprime[0];
   boxhi[1] = boxlo[1] + bprime[1];
   boxhi[2] = boxlo[2] + cprime[2];
@@ -667,28 +686,72 @@ void Domain::setup_general_triclinic(double *avec_caller, double *bvec_caller,
    transform atom coords from general triclinic to restricted triclinic
 ------------------------------------------------------------------------- */
 
-void Domain::general_to_restricted(double *x)
+void Domain::general_to_restricted_coords(double *x)
 {
-  double xnew[3];
+  double xshift[3],xnew[3];
   
-  MathExtra::matvec(rotate_g2r,x,xnew);
+  xshift[0] = x[0] - boxlo[0];
-  x[0] = xnew[0] + gtri_origin[0];
+  xshift[1] = x[1] - boxlo[1];
-  x[1] = xnew[1] + gtri_origin[1];
+  xshift[2] = x[2] - boxlo[2];
-  x[2] = xnew[2] + gtri_origin[2];
+  MathExtra::matvec(rotate_g2r,xshift,xnew);
+  x[0] = xnew[0] + boxlo[0];
+  x[1] = xnew[1] + boxlo[1];
+  x[2] = xnew[2] + boxlo[2];
 }
 
 /* ----------------------------------------------------------------------
    transform atom coords from restricted triclinic to general triclinic
 ------------------------------------------------------------------------- */
 
-void Domain::restricted_to_general(double *x)
+void Domain::restricted_to_general_coords(double *x)
 {
   double xshift[3],xnew[3];
   
-  xshift[0] = x[0] - gtri_origin[0];
+  xshift[0] = x[0] - boxlo[0];
-  xshift[1] = x[1] - gtri_origin[1];
+  xshift[1] = x[1] - boxlo[1];
-  xshift[2] = x[2] - gtri_origin[2];
+  xshift[2] = x[2] - boxlo[2];
   MathExtra::matvec(rotate_r2g,xshift,xnew);
+  x[0] = xnew[0] + boxlo[0];
+  x[1] = xnew[1] + boxlo[1];
+  x[2] = xnew[2] + boxlo[2];
+}
+
+void Domain::restricted_to_general_coords(double *x, double *xnew)
+{
+  double xshift[3];
+  
+  xshift[0] = x[0] - boxlo[0];
+  xshift[1] = x[1] - boxlo[1];
+  xshift[2] = x[2] - boxlo[2];
+  MathExtra::matvec(rotate_r2g,xshift,xnew);
+  xnew[0] += boxlo[0];
+  xnew[1] += boxlo[1];
+  xnew[2] += boxlo[2];
+}
+
+/* ----------------------------------------------------------------------
+   transform atom vector from general triclinic to restricted triclinic
+------------------------------------------------------------------------- */
+
+void Domain::general_to_restricted_vector(double *x)
+{
+  double xnew[3];
+  
+  MathExtra::matvec(rotate_g2r,x,xnew);
+  x[0] = xnew[0];
+  x[1] = xnew[1];
+  x[2] = xnew[2];
+}
+
+/* ----------------------------------------------------------------------
+   transform atom vector from restricted triclinic to general triclinic
+------------------------------------------------------------------------- */
+
+void Domain::restricted_to_general_vector(double *x)
+{
+  double xnew[3];
+  
+  MathExtra::matvec(rotate_r2g,x,xnew);
   x[0] = xnew[0];
   x[1] = xnew[1];
   x[2] = xnew[2];

src/domain.h

@@ -89,9 +89,9 @@ class Domain : protected Pointers {
   double h_rate[6], h_ratelo[3];    // rate of box size/shape change
 
   // general triclinic box
+  // boxlo = lower left corner
   
   double avec[3], bvec[3], cvec[3];  // ABC edge vectors of general triclinic box
-  double gtri_origin[3];             // origin of general triclinic box
   double rotate_g2r[3][3];           // rotation matrix from general --> restricted tri
   double rotate_r2g[3][3];           // rotation matrix from restricted --> general tri
 
@@ -141,9 +141,12 @@ class Domain : protected Pointers {
   void image_flip(int, int, int);
   int ownatom(int, double *, imageint *, int);
 
-  void setup_general_triclinic(double *, double *, double *, double *);
+  void define_general_triclinic(double *, double *, double *, double *);
-  void general_to_restricted(double *);
+  void general_to_restricted_coords(double *); 
-  void restricted_to_general(double *);
+  void restricted_to_general_coords(double *);
+  void restricted_to_general_coords(double *, double *);
+  void general_to_restricted_vector(double *); 
+  void restricted_to_general_vector(double *);
  
   void set_lattice(int, char **);
   void add_region(int, char **);

src/dump_custom.cpp

@@ -34,7 +34,7 @@
 using namespace LAMMPS_NS;
 
 // customize by adding keyword
-// also customize compute_atom_property.cpp
+// also customize compute_property_atom.cpp
 
 enum{ID,MOL,PROC,PROCP1,TYPE,ELEMENT,MASS,
      X,Y,Z,XS,YS,ZS,XSTRI,YSTRI,ZSTRI,XU,YU,ZU,XUTRI,YUTRI,ZUTRI,
@@ -88,6 +88,7 @@ DumpCustom::DumpCustom(LAMMPS *lmp, int narg, char **arg) :
   buffer_allow = 1;
   buffer_flag = 1;
 
+  triclinic_general = 0;
   nthresh = 0;
   nthreshlast = 0;
 
@@ -294,10 +295,14 @@ void DumpCustom::init_style()
 
   if (binary && domain->triclinic == 0)
     header_choice = &DumpCustom::header_binary;
+  else if (binary && triclinic_general == 1)
+    header_choice = &DumpCustom::header_binary_triclinic_general;
   else if (binary && domain->triclinic == 1)
     header_choice = &DumpCustom::header_binary_triclinic;
   else if (!binary && domain->triclinic == 0)
     header_choice = &DumpCustom::header_item;
+  else if (!binary && triclinic_general == 1)
+    header_choice = &DumpCustom::header_item_triclinic_general;
   else if (!binary && domain->triclinic == 1)
     header_choice = &DumpCustom::header_item_triclinic;
 
@@ -489,6 +494,31 @@ void DumpCustom::header_binary_triclinic(bigint ndump)
 
 /* ---------------------------------------------------------------------- */
 
+void DumpCustom::header_binary_triclinic_general(bigint ndump)
+{
+  header_format_binary();
+
+  fwrite(&update->ntimestep,sizeof(bigint),1,fp);
+  fwrite(&ndump,sizeof(bigint),1,fp);
+  fwrite(&domain->triclinic,sizeof(int),1,fp);
+  fwrite(&domain->triclinic_general,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(&nfield,sizeof(int),1,fp);
+
+  header_unit_style_binary();
+  header_time_binary();
+  header_columns_binary();
+
+  if (multiproc) fwrite(&nclusterprocs,sizeof(int),1,fp);
+  else fwrite(&nprocs,sizeof(int),1,fp);
+}
+
+/* ---------------------------------------------------------------------- */
+
 void DumpCustom::header_item(bigint ndump)
 {
   if (unit_flag && !unit_count) {
@@ -535,6 +565,32 @@ void DumpCustom::header_item_triclinic(bigint ndump)
 
 /* ---------------------------------------------------------------------- */
 
+void DumpCustom::header_item_triclinic_general(bigint ndump)
+{
+  if (unit_flag && !unit_count) {
+    ++unit_count;
+    fmt::print(fp,"ITEM: UNITS\n{}\n",update->unit_style);
+  }
+  if (time_flag) fmt::print(fp,"ITEM: TIME\n{:.16}\n",compute_time());
+
+  fmt::print(fp,"ITEM: TIMESTEP\n{}\n"
+             "ITEM: NUMBER OF ATOMS\n{}\n",
+             update->ntimestep, ndump);
+
+  fmt::print(fp,"ITEM: BOX BOUNDS abc origin {}\n"
+             "{:>1.16e} {:>1.16e} {:>1.16e} {:>1.16e}\n"
+             "{:>1.16e} {:>1.16e} {:>1.16e} {:>1.16e}\n"
+             "{:>1.16e} {:>1.16e} {:>1.16e} {:>1.16e}\n",
+             boundstr,
+             domain->avec[0],domain->avec[1],domain->avec[2],domain->boxlo[0],
+             domain->bvec[0],domain->bvec[1],domain->bvec[2],domain->boxlo[1],
+             domain->cvec[0],domain->cvec[1],domain->cvec[2],domain->boxlo[2]);
+
+  fmt::print(fp,"ITEM: ATOMS {}\n",columns);
+}
+
+/* ---------------------------------------------------------------------- */
+
 int DumpCustom::count()
 {
   int i;
@@ -1286,13 +1342,16 @@ int DumpCustom::parse_fields(int narg, char **arg)
       vtype[iarg] = Dump::DOUBLE;
 
     } else if (strcmp(arg[iarg],"x") == 0) {
-      pack_choice[iarg] = &DumpCustom::pack_x;
+      if (triclinic_general) pack_choice[iarg] = &DumpCustom::pack_x_triclinic_general;
+      else pack_choice[iarg] = &DumpCustom::pack_x;
       vtype[iarg] = Dump::DOUBLE;
     } else if (strcmp(arg[iarg],"y") == 0) {
-      pack_choice[iarg] = &DumpCustom::pack_y;
+      if (triclinic_general) pack_choice[iarg] = &DumpCustom::pack_y_triclinic_general;
+      else pack_choice[iarg] = &DumpCustom::pack_y;
       vtype[iarg] = Dump::DOUBLE;
     } else if (strcmp(arg[iarg],"z") == 0) { 
-      pack_choice[iarg] = &DumpCustom::pack_z;
+      if (triclinic_general) pack_choice[iarg] = &DumpCustom::pack_z_triclinic_general;
+      else pack_choice[iarg] = &DumpCustom::pack_z;
       vtype[iarg] = Dump::DOUBLE;
     } else if (strcmp(arg[iarg],"xs") == 0) {
       if (domain->triclinic) pack_choice[iarg] = &DumpCustom::pack_xs_triclinic;
@@ -1307,15 +1366,18 @@ int DumpCustom::parse_fields(int narg, char **arg)
       else pack_choice[iarg] = &DumpCustom::pack_zs;
       vtype[iarg] = Dump::DOUBLE;
     } else if (strcmp(arg[iarg],"xu") == 0) {
-      if (domain->triclinic) pack_choice[iarg] = &DumpCustom::pack_xu_triclinic;
+      if (triclinic_general) pack_choice[iarg] = &DumpCustom::pack_xu_triclinic_general;
+      else if (domain->triclinic) pack_choice[iarg] = &DumpCustom::pack_xu_triclinic;
       else pack_choice[iarg] = &DumpCustom::pack_xu;
       vtype[iarg] = Dump::DOUBLE;
     } else if (strcmp(arg[iarg],"yu") == 0) {
-      if (domain->triclinic) pack_choice[iarg] = &DumpCustom::pack_yu_triclinic;
+      if (triclinic_general) pack_choice[iarg] = &DumpCustom::pack_yu_triclinic_general;
+      else if (domain->triclinic) pack_choice[iarg] = &DumpCustom::pack_yu_triclinic;
       else pack_choice[iarg] = &DumpCustom::pack_yu;
       vtype[iarg] = Dump::DOUBLE;
     } else if (strcmp(arg[iarg],"zu") == 0) {
-      if (domain->triclinic) pack_choice[iarg] = &DumpCustom::pack_zu_triclinic;
+      if (triclinic_general) pack_choice[iarg] = &DumpCustom::pack_zu_triclinic_general;
+      else if (domain->triclinic) pack_choice[iarg] = &DumpCustom::pack_zu_triclinic;
       else pack_choice[iarg] = &DumpCustom::pack_zu;
       vtype[iarg] = Dump::DOUBLE;
     } else if (strcmp(arg[iarg],"xsu") == 0) {
@@ -1341,22 +1403,28 @@ int DumpCustom::parse_fields(int narg, char **arg)
       vtype[iarg] = Dump::INT;
 
     } else if (strcmp(arg[iarg],"vx") == 0) {
-      pack_choice[iarg] = &DumpCustom::pack_vx;
+      if (triclinic_general) pack_choice[iarg] = &DumpCustom::pack_vx_triclinic_general;
+      else pack_choice[iarg] = &DumpCustom::pack_vx;
       vtype[iarg] = Dump::DOUBLE;
     } else if (strcmp(arg[iarg],"vy") == 0) {
-      pack_choice[iarg] = &DumpCustom::pack_vy;
+      if (triclinic_general) pack_choice[iarg] = &DumpCustom::pack_vy_triclinic_general;
+      else pack_choice[iarg] = &DumpCustom::pack_vy;
       vtype[iarg] = Dump::DOUBLE;
     } else if (strcmp(arg[iarg],"vz") == 0) {
-      pack_choice[iarg] = &DumpCustom::pack_vz;
+      if (triclinic_general) pack_choice[iarg] = &DumpCustom::pack_vz_triclinic_general;
+      else pack_choice[iarg] = &DumpCustom::pack_vz;
       vtype[iarg] = Dump::DOUBLE;
     } else if (strcmp(arg[iarg],"fx") == 0) {
-      pack_choice[iarg] = &DumpCustom::pack_fx;
+      if (triclinic_general) pack_choice[iarg] = &DumpCustom::pack_fx_triclinic_general;
+      else pack_choice[iarg] = &DumpCustom::pack_fx;
       vtype[iarg] = Dump::DOUBLE;
     } else if (strcmp(arg[iarg],"fy") == 0) {
-      pack_choice[iarg] = &DumpCustom::pack_fy;
+      if (triclinic_general) pack_choice[iarg] = &DumpCustom::pack_fy_triclinic_general;
+      else pack_choice[iarg] = &DumpCustom::pack_fy;
       vtype[iarg] = Dump::DOUBLE;
     } else if (strcmp(arg[iarg],"fz") == 0) {
-      pack_choice[iarg] = &DumpCustom::pack_fz;
+      if (triclinic_general) pack_choice[iarg] = &DumpCustom::pack_fz_triclinic_general;
+      else pack_choice[iarg] = &DumpCustom::pack_fz;
       vtype[iarg] = Dump::DOUBLE;
 
     } else if (strcmp(arg[iarg],"q") == 0) {
@@ -1395,6 +1463,7 @@ int DumpCustom::parse_fields(int narg, char **arg)
         error->all(FLERR,"Dumping an atom property that isn't allocated");
       pack_choice[iarg] = &DumpCustom::pack_diameter;
       vtype[iarg] = Dump::DOUBLE;
+      
     } else if (strcmp(arg[iarg],"heatflow") == 0) {
       if (!atom->heatflow_flag)
         error->all(FLERR,"Dumping an atom property that isn't allocated");
@@ -1405,6 +1474,7 @@ int DumpCustom::parse_fields(int narg, char **arg)
         error->all(FLERR,"Dumping an atom property that isn't allocated");
       pack_choice[iarg] = &DumpCustom::pack_temperature;
       vtype[iarg] = Dump::DOUBLE;
+      
     } else if (strcmp(arg[iarg],"omegax") == 0) {
       if (!atom->omega_flag)
         error->all(FLERR,"Dumping an atom property that isn't allocated");
@@ -1702,6 +1772,14 @@ int DumpCustom::modify_param(int narg, char **arg)
     return 2;
   }
 
+  if (strcmp(arg[0],"triclinic/general") == 0) {
+    if (narg < 2) error->all(FLERR,"Illegal dump_modify command");
+    triclinic_general = utils::logical(FLERR,arg[1],false,lmp);
+    if (triclinic_general && !domain->triclinic_general)
+      error->all(FLERR,"Dump_modify triclinic/general invalid b/c simulation box is not");
+    return 2;
+  }
+
   if (strcmp(arg[0],"format") == 0) {
     if (narg < 2) utils::missing_cmd_args(FLERR, "dump_modify format", error);
 
@@ -2283,6 +2361,48 @@ void DumpCustom::pack_z(int n)
 
 /* ---------------------------------------------------------------------- */
 
+void DumpCustom::pack_x_triclinic_general(int n)
+{
+  double **x = atom->x;
+  double xtri[3];
+  
+  for (int i = 0; i < nchoose; i++) {
+    domain->restricted_to_general_coords(x[clist[i]],xtri);
+    buf[n] = xtri[0];
+    n += size_one;
+  }
+}
+
+/* ---------------------------------------------------------------------- */
+
+void DumpCustom::pack_y_triclinic_general(int n)
+{
+  double **x = atom->x;
+  double xtri[3];
+
+  for (int i = 0; i < nchoose; i++) {
+    domain->restricted_to_general_coords(x[clist[i]],xtri);
+    buf[n] = xtri[1];
+    n += size_one;
+  }
+}
+
+/* ---------------------------------------------------------------------- */
+
+void DumpCustom::pack_z_triclinic_general(int n)
+{
+  double **x = atom->x;
+  double xtri[3];
+
+  for (int i = 0; i < nchoose; i++) {
+    domain->restricted_to_general_coords(x[clist[i]],xtri);
+    buf[n] = xtri[2];
+    n += size_one;
+  }
+}
+
+/* ---------------------------------------------------------------------- */
+
 void DumpCustom::pack_xs(int n)
 {
   double **x = atom->x;
@@ -2489,6 +2609,84 @@ void DumpCustom::pack_zu_triclinic(int n)
 
 /* ---------------------------------------------------------------------- */
 
+void DumpCustom::pack_xu_triclinic_general(int n)
+{
+  int j;
+  double **x = atom->x;
+  imageint *image = atom->image;
+
+  double *h = domain->h;
+  double xu[3];
+  int xbox,ybox,zbox;
+
+  for (int i = 0; i < nchoose; i++) {
+    j = clist[i];
+    xbox = (image[j] & IMGMASK) - IMGMAX;
+    ybox = (image[j] >> IMGBITS & IMGMASK) - IMGMAX;
+    zbox = (image[j] >> IMG2BITS) - IMGMAX;
+    xu[0] = x[j][0] + h[0]*xbox + h[5]*ybox + h[4]*zbox;
+    xu[1] = x[j][1] + h[1]*ybox + h[3]*zbox;
+    xu[2] = x[j][2] + h[2]*zbox;
+    domain->restricted_to_general_coords(xu);
+    buf[n] = xu[0];
+    n += size_one;
+  }
+}
+
+/* ---------------------------------------------------------------------- */
+
+void DumpCustom::pack_yu_triclinic_general(int n)
+{
+  int j;
+  double **x = atom->x;
+  imageint *image = atom->image;
+
+  double *h = domain->h;
+  double xu[3];
+  int xbox,ybox,zbox;
+
+  for (int i = 0; i < nchoose; i++) {
+    j = clist[i];
+    xbox = (image[j] & IMGMASK) - IMGMAX;
+    ybox = (image[j] >> IMGBITS & IMGMASK) - IMGMAX;
+    zbox = (image[j] >> IMG2BITS) - IMGMAX;
+    xu[0] = x[j][0] + h[0]*xbox + h[5]*ybox + h[4]*zbox;
+    xu[1] = x[j][1] + h[1]*ybox + h[3]*zbox;
+    xu[2] = x[j][2] + h[2]*zbox;
+    domain->restricted_to_general_coords(xu);
+    buf[n] = xu[1];
+    n += size_one;
+  }
+}
+
+/* ---------------------------------------------------------------------- */
+
+void DumpCustom::pack_zu_triclinic_general(int n)
+{
+  int j;
+  double **x = atom->x;
+  imageint *image = atom->image;
+
+  double *h = domain->h;
+  double xu[3];
+  int xbox,ybox,zbox;
+
+  for (int i = 0; i < nchoose; i++) {
+    j = clist[i];
+    xbox = (image[j] & IMGMASK) - IMGMAX;
+    ybox = (image[j] >> IMGBITS & IMGMASK) - IMGMAX;
+    zbox = (image[j] >> IMG2BITS) - IMGMAX;
+    xu[0] = x[j][0] + h[0]*xbox + h[5]*ybox + h[4]*zbox;
+    xu[1] = x[j][1] + h[1]*ybox + h[3]*zbox;
+    xu[2] = x[j][2] + h[2]*zbox;
+    domain->restricted_to_general_coords(xu);
+    buf[n] = xu[2];
+    n += size_one;
+  }
+}
+
+/* ---------------------------------------------------------------------- */
+
 void DumpCustom::pack_xsu(int n)
 {
   int j;
@@ -2671,6 +2869,48 @@ void DumpCustom::pack_vz(int n)
 
 /* ---------------------------------------------------------------------- */
 
+void DumpCustom::pack_vx_triclinic_general(int n)
+{
+  double **v = atom->v;
+  double vtri[3];
+
+  for (int i = 0; i < nchoose; i++) {
+    domain->restricted_to_general_coords(v[clist[i]],vtri);
+    buf[n] = vtri[0];
+    n += size_one;
+  }
+}
+
+/* ---------------------------------------------------------------------- */
+
+void DumpCustom::pack_vy_triclinic_general(int n)
+{
+  double **v = atom->v;
+  double vtri[3];
+
+  for (int i = 0; i < nchoose; i++) {
+    domain->restricted_to_general_coords(v[clist[i]],vtri);
+    buf[n] = vtri[1];
+    n += size_one;
+  }
+}
+
+/* ---------------------------------------------------------------------- */
+
+void DumpCustom::pack_vz_triclinic_general(int n)
+{
+  double **v = atom->v;
+  double vtri[3];
+
+  for (int i = 0; i < nchoose; i++) {
+    domain->restricted_to_general_coords(v[clist[i]],vtri);
+    buf[n] = vtri[2];
+    n += size_one;
+  }
+}
+
+/* ---------------------------------------------------------------------- */
+
 void DumpCustom::pack_fx(int n)
 {
   double **f = atom->f;
@@ -2707,6 +2947,48 @@ void DumpCustom::pack_fz(int n)
 
 /* ---------------------------------------------------------------------- */
 
+void DumpCustom::pack_fx_triclinic_general(int n)
+{
+  double **f = atom->f;
+  double ftri[3];
+  
+  for (int i = 0; i < nchoose; i++) {
+    domain->restricted_to_general_coords(f[clist[i]],ftri);
+    buf[n] = ftri[0];
+    n += size_one;
+  }
+}
+
+/* ---------------------------------------------------------------------- */
+
+void DumpCustom::pack_fy_triclinic_general(int n)
+{
+  double **f = atom->f;
+  double ftri[3];
+
+  for (int i = 0; i < nchoose; i++) {
+    domain->restricted_to_general_coords(f[clist[i]],ftri);
+    buf[n] = ftri[1];
+    n += size_one;
+  }
+}
+
+/* ---------------------------------------------------------------------- */
+
+void DumpCustom::pack_fz_triclinic_general(int n)
+{
+  double **f = atom->f;
+  double ftri[3];
+
+  for (int i = 0; i < nchoose; i++) {
+    domain->restricted_to_general_coords(f[clist[i]],ftri);
+    buf[n] = ftri[2];
+    n += size_one;
+  }
+}
+
+/* ---------------------------------------------------------------------- */
+
 void DumpCustom::pack_q(int n)
 {
   double *q = atom->q;

src/dump_custom.h

@@ -37,6 +37,8 @@ class DumpCustom : public Dump {
   int nevery;        // dump frequency for output
   char *idregion;    // region ID, nullptr if no region
 
+  int triclinic_general;              // set by dump_modify
+  
   int nthresh;                        // # of defined thresholds
   int nthreshlast;                    // # of defined thresholds with value = LAST
                                       //
@@ -124,8 +126,10 @@ class DumpCustom : public Dump {
   FnPtrHeader header_choice;    // ptr to write header functions
   void header_binary(bigint);
   void header_binary_triclinic(bigint);
+  void header_binary_triclinic_general(bigint);
   void header_item(bigint);
   void header_item_triclinic(bigint);
+  void header_item_triclinic_general(bigint);
 
   typedef void (DumpCustom::*FnPtrWrite)(int, double *);
   FnPtrWrite write_choice;    // ptr to write data functions
@@ -153,24 +157,34 @@ class DumpCustom : public Dump {
   void pack_x(int);
   void pack_y(int);
   void pack_z(int);
+  void pack_x_triclinic_general(int);
+  void pack_y_triclinic_general(int);
+  void pack_z_triclinic_general(int);
+  
   void pack_xs(int);
   void pack_ys(int);
   void pack_zs(int);
   void pack_xs_triclinic(int);
   void pack_ys_triclinic(int);
   void pack_zs_triclinic(int);
+
   void pack_xu(int);
   void pack_yu(int);
   void pack_zu(int);
   void pack_xu_triclinic(int);
   void pack_yu_triclinic(int);
   void pack_zu_triclinic(int);
+  void pack_xu_triclinic_general(int);
+  void pack_yu_triclinic_general(int);
+  void pack_zu_triclinic_general(int);
+
   void pack_xsu(int);
   void pack_ysu(int);
   void pack_zsu(int);
   void pack_xsu_triclinic(int);
   void pack_ysu_triclinic(int);
   void pack_zsu_triclinic(int);
+
   void pack_ix(int);
   void pack_iy(int);
   void pack_iz(int);
@@ -178,9 +192,17 @@ class DumpCustom : public Dump {
   void pack_vx(int);
   void pack_vy(int);
   void pack_vz(int); 
+  void pack_vx_triclinic_general(int);
+  void pack_vy_triclinic_general(int);
+  void pack_vz_triclinic_general(int);
+   
   void pack_fx(int);
   void pack_fy(int);
   void pack_fz(int);
+  void pack_fx_triclinic_general(int);
+  void pack_fy_triclinic_general(int);
+  void pack_fz_triclinic_general(int);
+
   void pack_q(int);
   void pack_mux(int);
   void pack_muy(int);
@@ -188,6 +210,7 @@ class DumpCustom : public Dump {
   void pack_mu(int);
   void pack_radius(int);
   void pack_diameter(int);
+  
   void pack_heatflow(int);
   void pack_temperature(int);
 

src/lmprestart.h

@@ -29,6 +29,7 @@ enum{VERSION,SMALLINT,TAGINT,BIGINT,
      NDIHEDRALS,NDIHEDRALTYPES,DIHEDRAL_PER_ATOM,
      NIMPROPERS,NIMPROPERTYPES,IMPROPER_PER_ATOM,
      TRICLINIC,BOXLO,BOXHI,XY,XZ,YZ,
+     TRICLINIC_GENERAL,ROTATE_G2R,ROTATE_R2G,
      SPECIAL_LJ,SPECIAL_COUL,
      MASS,PAIR,BOND,ANGLE,DIHEDRAL,IMPROPER,
      MULTIPROC,MPIIO,PROCSPERFILE,PERPROC,

src/pair_lj_cut.cpp

@@ -111,8 +111,6 @@ void PairLJCut::compute(int eflag, int vflag)
       jtype = type[j];
 
       if (rsq < cutsq[itype][jtype]) {
-        printf("AAA tags %d %d ij %d %d nlocal %d dist %g\n",
-               atom->tag[i],atom->tag[j],i,j,atom->nlocal,sqrt(rsq));
         r2inv = 1.0 / rsq;
         r6inv = r2inv * r2inv * r2inv;
         forcelj = r6inv * (lj1[itype][jtype] * r6inv - lj2[itype][jtype]);

src/read_data.cpp

@@ -480,7 +480,7 @@ void ReadData::command(int narg, char **arg)
   ellipsoidflag = lineflag = triflag = bodyflag = 0;
 
   xloxhi_flag = yloyhi_flag = zlozhi_flag = tilt_flag = 0;
-  avec_flag = bvec_flag = cvec_flag = gtri_origin_flag = 0;
+  avec_flag = bvec_flag = cvec_flag = abc_origin_flag = 0;
 
   // values in this data file
 
@@ -495,7 +495,7 @@ void ReadData::command(int narg, char **arg)
   bvec[0] = bvec[1] = bvec[2] = 0.0;
   cvec[0] = cvec[1] = cvec[2] = 0.0;
   avec[0] = bvec[1] = cvec[2] = 1.0;
-  gtri_origin[0] = gtri_origin[1] = gtri_origin[2] = 0.0;
+  abc_origin[0] = abc_origin[1] = abc_origin[2] = 0.0;
   
   keyword[0] = '\0';
 
@@ -518,7 +518,7 @@ void ReadData::command(int narg, char **arg)
 
     // check if simulation box specified consistently
 
-    if (!avec_flag && !bvec_flag && !cvec_flag && !gtri_origin_flag) {
+    if (!avec_flag && !bvec_flag && !cvec_flag && !abc_origin_flag) {
       triclinic = triclinic_general = 0;
       if (tilt_flag) triclinic = 1;
     } else {
@@ -579,11 +579,11 @@ void ReadData::command(int narg, char **arg)
         }
 
       // general triclinic box
-      // setup_general_triclinic() converts
+      // define_general_triclinic() converts
-      //   ABC edge vectors + gtri_origin to restricted triclinic
+      //   ABC edge vectors + abc_origin to restricted triclinic
         
       } else if (triclinic_general) {
-        domain->setup_general_triclinic(avec,bvec,cvec,gtri_origin);
+        domain->define_general_triclinic(avec,bvec,cvec,abc_origin);
       }
     }
 
@@ -596,7 +596,7 @@ void ReadData::command(int narg, char **arg)
 
       // general triclinic
       // first data file must also be general triclinic
-      // avec,bvec,vec must match first data file
+      // avec,bvec,vec and origin must match first data file
       // shift not allowed
       
       if (triclinic_general) {
@@ -609,8 +609,8 @@ void ReadData::command(int narg, char **arg)
           errflag = 1;
         if (cvec[0] != domain->cvec[0] || cvec[1] != domain->cvec[1] || cvec[2] != domain->cvec[2])
           errflag = 1;
-        if (gtri_origin[0] != domain->gtri_origin[0] || gtri_origin[1] != domain->gtri_origin[1] ||
+        if (abc_origin[0] != domain->boxlo[0] || abc_origin[1] != domain->boxlo[1] ||
-            gtri_origin[2] != domain->gtri_origin[2])
+            abc_origin[2] != domain->boxlo[2])
           errflag = 1;
         if (errflag)
           error->all(FLERR,"Read_data subsequent file ABC vectors must be same as first file");
@@ -1432,11 +1432,11 @@ void ReadData::header(int firstpass)
       cvec[1] = utils::numeric(FLERR, words[1], false, lmp);
       cvec[2] = utils::numeric(FLERR, words[2], false, lmp);
 
-    } else if (utils::strmatch(line, "^\\s*\\f+\\s+\\f+\\s+\\f+\\s+\\gtri\\s+origin\\s")) {
+    } else if (utils::strmatch(line, "^\\s*\\f+\\s+\\f+\\s+\\f+\\s+\\abc\\s+origin\\s")) {
-      gtri_origin_flag = 1;
+      abc_origin_flag = 1;
-      gtri_origin[0] = utils::numeric(FLERR, words[0], false, lmp);
+      abc_origin[0] = utils::numeric(FLERR, words[0], false, lmp);
-      gtri_origin[1] = utils::numeric(FLERR, words[1], false, lmp);
+      abc_origin[1] = utils::numeric(FLERR, words[1], false, lmp);
-      gtri_origin[2] = utils::numeric(FLERR, words[2], false, lmp);
+      abc_origin[2] = utils::numeric(FLERR, words[2], false, lmp);
 
     } else
       break;

src/read_data.h

@@ -63,10 +63,11 @@ class ReadData : public Command {
 
   double boxlo[3], boxhi[3];
   double xy, xz, yz;
-  double avec[3], bvec[3], cvec[3], gtri_origin[3];
+  double avec[3], bvec[3], cvec[3];
+  double abc_origin[3];
   int triclinic, triclinic_general;
   int xloxhi_flag, yloyhi_flag, zlozhi_flag, tilt_flag;
-  int avec_flag, bvec_flag, cvec_flag, gtri_origin_flag;
+  int avec_flag, bvec_flag, cvec_flag, abc_origin_flag;
   
   // optional args
 

src/read_restart.cpp

@@ -783,6 +783,13 @@ void ReadRestart::header()
     } else if (flag == YZ) {
       domain->yz = read_double();
 
+    } else if (flag == TRICLINIC_GENERAL) {
+      domain->triclinic_general = read_int();
+    } else if (flag == ROTATE_G2R) {
+      read_double_vec(9,&domain->rotate_g2r[0][0]);
+    } else if (flag == ROTATE_R2G) {
+      read_double_vec(9,&domain->rotate_r2g[0][0]);
+
     } else if (flag == SPECIAL_LJ) {
       read_int();
       read_double_vec(3,&force->special_lj[1]);

src/write_data.cpp

@@ -227,7 +227,7 @@ void WriteData::write(const std::string &file)
     memory->create(xstore,nlocal,3,"write_data:xstore");
     if (nlocal) memcpy(&xstore[0][0],&x[0][0],3*nlocal*sizeof(double));
     for (int i = 0; i < nlocal; i++)
-      domain->restricted_to_general(x[i]);
+      domain->restricted_to_general_coords(x[i]);
   }
   
   if (natoms) atoms();
@@ -335,8 +335,8 @@ void WriteData::header()
                domain->avec[0],domain->avec[1],domain->avec[2],
                domain->bvec[0],domain->bvec[1],domain->bvec[2],
                domain->cvec[0],domain->cvec[1],domain->cvec[2]);
-    fmt::print(fp,"{} {} {} gtri origin\n",
+    fmt::print(fp,"{} {} {} abc origin\n",
-               domain->gtri_origin[0],domain->gtri_origin[1],domain->gtri_origin[2]);
+               domain->boxlo[0],domain->boxlo[1],domain->boxhi[2]);
   }
 }
 

src/write_restart.cpp

@@ -448,6 +448,12 @@ void WriteRestart::header()
   write_double(XZ,domain->xz);
   write_double(YZ,domain->yz);
 
+  write_int(TRICLINIC_GENERAL,domain->triclinic_general);
+  if (domain->triclinic_general) {
+    write_double_vec(ROTATE_G2R,9,&domain->rotate_g2r[0][0]);
+    write_double_vec(ROTATE_R2G,9,&domain->rotate_r2g[0][0]);
+  }
+  
   write_double_vec(SPECIAL_LJ,3,&force->special_lj[1]);
   write_double_vec(SPECIAL_COUL,3,&force->special_coul[1]);