diff --git a/src/atom.cpp b/src/atom.cpp
index c4521a244e..0249547253 100644
--- a/src/atom.cpp
+++ b/src/atom.cpp
@@ -1041,7 +1041,8 @@ void Atom::deallocate_topology()
 ------------------------------------------------------------------------- */
 
 void Atom::data_atoms(int n, char *buf, tagint id_offset, tagint mol_offset,
-                      int type_offset, int shiftflag, double *shift,
+                      int type_offset, int triclinic_general,
+                      int shiftflag, double *shift,
                       int labelflag, int *ilabel)
 {
   int xptr,iptr;
@@ -1053,6 +1054,7 @@ void Atom::data_atoms(int n, char *buf, tagint id_offset, tagint mol_offset,
   auto location = "Atoms section of data file";
 
   // use the first line to detect and validate the number of words/tokens per line
+  
   next = strchr(buf,'\n');
   if (!next) error->all(FLERR, "Missing data in {}", location);
   *next = '\0';
@@ -1069,6 +1071,7 @@ void Atom::data_atoms(int n, char *buf, tagint id_offset, tagint mol_offset,
     error->all(FLERR,"Incorrect format in {}: {}", location, utils::trim(buf));
 
   *next = '\n';
+  
   // set bounds for my proc
   // if periodic and I am lo/hi proc, adjust bounds by EPSILON
   // ensures all data atoms will be owned even with round-off
@@ -1143,11 +1146,19 @@ void Atom::data_atoms(int n, char *buf, tagint id_offset, tagint mol_offset,
     *next = '\0';
     auto values = Tokenizer(buf).as_vector();
     int nvalues = values.size();
-    if ((nvalues == 0) || (utils::strmatch(values[0],"^#.*")))  {
-      // skip over empty or comment lines
+
+    // skip comment lines
+    
+    if ((nvalues == 0) || (utils::strmatch(values[0],"^#.*"))) {
+
+    // check that line has correct # of words
+        
     } else if ((nvalues < nwords) ||
                ((nvalues > nwords) && (!utils::strmatch(values[nwords],"^#")))) {
       error->all(FLERR, "Incorrect format in {}: {}", location, utils::trim(buf));
+
+    // extract the atom coords and image flags (if they exist)
+      
     } else {
       int imx = 0, imy = 0, imz = 0;
       if (imageflag) {
@@ -1167,13 +1178,25 @@ void Atom::data_atoms(int n, char *buf, tagint id_offset, tagint mol_offset,
       xdata[0] = utils::numeric(FLERR,values[xptr],false,lmp);
       xdata[1] = utils::numeric(FLERR,values[xptr+1],false,lmp);
       xdata[2] = utils::numeric(FLERR,values[xptr+2],false,lmp);
+
+      // convert atom coords from general triclinic to restricted triclinic
+      
+      if (triclinic_general) domain->general_to_restricted(xdata);
+
+      // apply shift if requested by read_data command
+      
       if (shiftflag) {
         xdata[0] += shift[0];
         xdata[1] += shift[1];
         xdata[2] += shift[2];
       }
 
+      // map atom into simulation box for periodic dimensions
+      
       domain->remap(xdata,imagedata);
+
+      // determine if this proc owns the atom
+      
       if (triclinic) {
         domain->x2lamda(xdata,lamda);
         coord = lamda;
diff --git a/src/atom.h b/src/atom.h
index 548168ac59..cf5fa10814 100644
--- a/src/atom.h
+++ b/src/atom.h
@@ -328,7 +328,8 @@ class Atom : protected Pointers {
 
   void deallocate_topology();
 
-  void data_atoms(int, char *, tagint, tagint, int, int, double *, int, int *);
+  void data_atoms(int, char *, tagint, tagint, int, int,
+                  int, double *, int, int *);
   void data_vels(int, char *, tagint);
   void data_bonds(int, char *, int *, tagint, int, int, int *);
   void data_angles(int, char *, int *, tagint, int, int, int *);
diff --git a/src/domain.cpp b/src/domain.cpp
index 3627af26cf..6523bfc8c7 100644
--- a/src/domain.cpp
+++ b/src/domain.cpp
@@ -28,6 +28,7 @@
 #include "force.h"
 #include "kspace.h"
 #include "lattice.h"
+#include "math_extra.h"
 #include "memory.h"
 #include "modify.h"
 #include "molecule.h"
@@ -41,6 +42,7 @@
 #include <cmath>
 
 using namespace LAMMPS_NS;
+using namespace MathExtra;
 
 #define BIG   1.0e20
 #define SMALL 1.0e-4
@@ -81,7 +83,7 @@ Domain::Domain(LAMMPS *lmp) : Pointers(lmp)
   minylo = minyhi = 0.0;
   minzlo = minzhi = 0.0;
 
-  triclinic = 0;
+  triclinic = triclinic_general = 0;
 
   boxlo[0] = boxlo[1] = boxlo[2] = -0.5;
   boxhi[0] = boxhi[1] = boxhi[2] = 0.5;
@@ -516,6 +518,180 @@ void Domain::reset_box()
   }
 }
 
+/* ----------------------------------------------------------------------
+   define and store a general triclinic simulation box
+   3 edge vectors of box = avec/bvec/cvec caller
+   origin of edge vectors = origin_caller = lower left corner of box
+   create mapping to restricted triclinic box
+   set boxlo[3], boxhi[3] and 3 tilt factors
+   create rotation matrices for general <--> restricted transformations
+------------------------------------------------------------------------- */
+
+void Domain::set_general_triclinic(double *avec_caller, double *bvec_caller,
+                                   double *cvec_caller, double *origin_caller)
+{
+  if (triclinic || triclinic_general)
+    error->all(FLERR,"General triclinic box edge vectors are already set");
+  
+  triclinic = triclinic_general = 1;
+
+  avec[0] = avec_caller[0];
+  avec[1] = avec_caller[1];
+  avec[2] = avec_caller[2];
+
+  bvec[0] = bvec_caller[0];
+  bvec[1] = bvec_caller[1];
+  bvec[2] = bvec_caller[2];
+
+  cvec[0] = cvec_caller[0];
+  cvec[1] = cvec_caller[1];
+  cvec[2] = cvec_caller[2];
+
+  tri_origin[0] = origin_caller[0];
+  tri_origin[1] = origin_caller[1];
+  tri_origin[2] = origin_caller[2];
+
+  // error check for co-planar A,B,C
+
+  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");
+    
+  // quat1 = convert A into A' along x-axis
+  // rot1 = unit vector to rotate A around
+  // theta1 = angle of rotation calculated from
+  //   A dot xunit = Ax = |A| cos(theta1)
+
+  double rot1[3],quat1[4];
+  double xaxis[3] = {1.0, 0.0, 0.0};
+
+  double avec_len = MathExtra::len3(avec);
+  MathExtra::cross3(avec,xaxis,rot1);
+  MathExtra::norm3(rot1);
+  double theta1 = acos(avec[0]/avec_len);
+  MathExtra::axisangle_to_quat(rot1,theta1,quat1);
+
+  // rotmat1 = rotation matrix associated with quat1
+
+  double rotmat1[3][3];
+  MathExtra::quat_to_mat(quat1,rotmat1);
+
+  // B1 = rotation of B by quat1 rotation matrix
+
+  double bvec1[3];
+
+  MathExtra::matvec(rotmat1,bvec,bvec1);
+
+  // quat2 = rotation to convert B1 into B' in xy plane
+  // Byz1 = projection of B1 into yz plane
+  // rot2 = unit vector to rotate B1 around = -x axis
+  // theta2 = angle of rotation calculated from
+  //   Byz1 dot yunit = B1y = |Byz1} cos(theta2)
+
+  double byzvec1[3],quat2[4];
+  MathExtra::copy3(bvec1,byzvec1);
+  byzvec1[0] = 0.0;
+  double byzvec1_len = MathExtra::len3(byzvec1);
+  double rot2[3] = {-1.0, 0.0, 0.0};
+  double theta2 = acos(bvec1[1]/byzvec1_len);
+  MathExtra::axisangle_to_quat(rot2,theta2,quat2);
+
+  // quat = product of quat2 * quat1 = transformation via single quat
+  // rotate_g2r = general to restricted transformation matrix
+  // rotate_r2g = restricted to general transformation matrix
+  // if A x B not in direction of C, flip sign of z component of transform
+  //   done by flipping sign of 3rd row of rotate_g2r matrix
+
+  double quat[4];
+  MathExtra::quatquat(quat2,quat1,quat);
+  MathExtra::quat_to_mat(quat,rotate_g2r);
+
+  if (dot < 0.0) {
+    rotate_g2r[2][0] = -rotate_g2r[2][0];
+    rotate_g2r[2][1] = -rotate_g2r[2][1];
+    rotate_g2r[2][2] = -rotate_g2r[2][2];
+  }
+  
+  MathExtra::transpose3(rotate_g2r,rotate_r2g);
+  
+  // A',B',C' = transformation of A,B,C to restricted triclinic
+  
+  double aprime[3],bprime[3],cprime[3];
+  MathExtra::matvec(rotate_g2r,avec,aprime);
+  MathExtra::matvec(rotate_g2r,bvec,bprime);
+  MathExtra::matvec(rotate_g2r,cvec,cprime);
+
+  // set restricted triclinic boxlo, boxhi, and tilt factors
+
+  boxlo[0] = tri_origin[0];
+  boxlo[1] = tri_origin[1];
+  boxlo[2] = tri_origin[2];
+
+  boxhi[0] = boxlo[0] + aprime[0];
+  boxhi[1] = boxlo[1] + bprime[1];
+  boxhi[2] = boxlo[2] + cprime[2];
+
+  xy = bprime[1];
+  xz = cprime[0];
+  yz = cprime[1];
+  
+  // debug
+
+  /*
+  printf("Quat: %g %g %g %g\n",quat[0],quat[1],quat[2],quat[3]);
+  double angle = 2.0*acos(quat[0]);
+  printf("Theta: %g\n",angle);
+  printf("Rotvec: %g %g %g\n",quat[1]/sin(0.5*angle),quat[2]/sin(0.5*angle),quat[3]/sin(0.5*angle));
+  printf("Aprime: %g %g %g\n",aprime[0],aprime[1],aprime[2]);
+  printf("Bprime: %g %g %g\n",bprime[0],bprime[1],bprime[2]);
+  printf("Cprime: %g %g %g\n",cprime[0],cprime[1],cprime[2]);
+  printf("Length A: %g %g\n",MathExtra::len3(avec),MathExtra::len3(aprime));
+  printf("Length B: %g %g\n",MathExtra::len3(bvec),MathExtra::len3(bprime));
+  printf("Length C: %g %g\n",MathExtra::len3(cvec),MathExtra::len3(cprime));
+
+  double coord1[3] = {0.5,0.0,0.0};
+  double coord2[3] = {0.5,0.0,0.3};
+  double newcoord[3];
+  MathExtra::matvec(rotate_g2r,coord1,newcoord);
+  printf("Atom1: %g %g %g\n",newcoord[0],newcoord[1],newcoord[2]);
+  MathExtra::matvec(rotate_g2r,coord2,newcoord);
+  printf("Atom2: %g %g %g\n",newcoord[0],newcoord[1],newcoord[2]);
+  */
+}
+
+/* ----------------------------------------------------------------------
+   transform one atom's coords from general triclinic to restricted triclinic
+------------------------------------------------------------------------- */
+
+void Domain::general_to_restricted(double *x)
+{
+  double xnew[3];
+  
+  MathExtra::matvec(rotate_g2r,x,xnew);
+  x[0] = xnew[0] + tri_origin[0];
+  x[1] = xnew[1] + tri_origin[1];
+  x[2] = xnew[2] + tri_origin[2];
+}
+
+/* ----------------------------------------------------------------------
+   transform one atom's coords from restricted triclinic to general triclinic
+------------------------------------------------------------------------- */
+
+void Domain::restricted_to_general(double *x)
+{
+  double xshift[3],xnew[3];
+  
+  xshift[0] = x[0] - tri_origin[0];
+  xshift[1] = x[1] - tri_origin[1];
+  xshift[2] = x[2] - tri_origin[2];
+  MathExtra::matvec(rotate_r2g,xshift,xnew);
+  x[0] = xnew[0];
+  x[1] = xnew[1];
+  x[2] = xnew[2];
+}
+
 /* ----------------------------------------------------------------------
    enforce PBC and modify box image flags for each atom
    called every reneighboring and by other commands that change atoms
@@ -676,9 +852,7 @@ int Domain::inside(double* x)
         lamda[1] < lo[1] || lamda[1] >= hi[1] ||
         lamda[2] < lo[2] || lamda[2] >= hi[2]) return 0;
     else return 1;
-
   }
-
 }
 
 /* ----------------------------------------------------------------------
@@ -713,7 +887,6 @@ int Domain::inside_nonperiodic(double* x)
     if (!zperiodic && (lamda[2] < lo[2] || lamda[2] >= hi[2])) return 0;
     return 1;
   }
-
 }
 
 /* ----------------------------------------------------------------------
diff --git a/src/domain.h b/src/domain.h
index ab054f1b50..e77af5d968 100644
--- a/src/domain.h
+++ b/src/domain.h
@@ -39,8 +39,9 @@ class Domain : protected Pointers {
                          // 2 = shrink-wrap non-periodic
                          // 3 = shrink-wrap non-per w/ min
 
-  int triclinic;    // 0 = orthog box, 1 = triclinic
-
+  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
+  
   // orthogonal box
 
   double xprd, yprd, zprd;                   // global box dimensions
@@ -48,8 +49,8 @@ class Domain : protected Pointers {
   double prd[3];                             // array form of dimensions
   double prd_half[3];                        // array form of half dimensions
 
-  // triclinic box
-  // xyzprd,xyzprd_half and prd,prd_half = same as if untilted
+  // restricted triclinic box
+  // xyzprd,xyzprd_half and prd,prd_half = same as if not tilted
 
   double prd_lamda[3];         // lamda box = (1,1,1)
   double prd_half_lamda[3];    // lamda half box = (0.5,0.5,0.5)
@@ -58,14 +59,14 @@ class Domain : protected Pointers {
 
   double boxlo[3], boxhi[3];
 
-  // triclinic box
-  // boxlo/hi = same as if untilted
+  // restricted triclinic box
+  // boxlo/hi = same as if not tilted
 
   double boxlo_lamda[3], boxhi_lamda[3];    // lamda box = (0,1)
   double boxlo_bound[3], boxhi_bound[3];    // bounding box of tilted domain
   double corners[8][3];                     // 8 corner points
 
-  // orthogonal box & triclinic box
+  // orthogonal box & restricted triclinic box
 
   double minxlo, minxhi;    // minimum size of global box
   double minylo, minyhi;    // when shrink-wrapping
@@ -75,18 +76,27 @@ class Domain : protected Pointers {
 
   double sublo[3], subhi[3];    // sub-box bounds on this proc
 
-  // triclinic box
+  // restricted triclinic box
   // sublo/hi = undefined
 
   double sublo_lamda[3], subhi_lamda[3];    // bounds of subbox in lamda
 
-  // triclinic box
+  // restricted triclinic box
 
   double xy, xz, yz;                // 3 tilt factors
   double h[6], h_inv[6];            // shape matrix in Voigt ordering
                                     // Voigt = xx,yy,zz,yz,xz,xy
   double h_rate[6], h_ratelo[3];    // rate of box size/shape change
 
+  // general triclinic box
+
+  double avec[3], bvec[3], cvec[3];  // ABC edge vectors of general triclinic box
+  double tri_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
+
+  // box flags
+  
   int box_change;           // 1 if any of next 3 flags are set, else 0
   int box_change_size;      // 1 if box size changes, 0 if not
   int box_change_shape;     // 1 if box shape changes, 0 if not
@@ -131,6 +141,10 @@ class Domain : protected Pointers {
   void image_flip(int, int, int);
   int ownatom(int, double *, imageint *, int);
 
+  void set_general_triclinic(double *, double *, double *, double *);
+  void general_to_restricted(double *);
+  void restricted_to_general(double *);
+  
   void set_lattice(int, char **);
   void add_region(int, char **);
   void delete_region(Region *);
diff --git a/src/math_extra.h b/src/math_extra.h
index 1efacea463..49e128c3df 100644
--- a/src/math_extra.h
+++ b/src/math_extra.h
@@ -47,6 +47,7 @@ inline void cross3(const double *v1, const double *v2, double *ans);
 
 inline void zeromat3(double m[3][3]);
 inline void zeromat3(double **m);
+inline void transpose3(const double m[3][3], double ans[3][3]);
 
 inline void col2mat(const double *ex, const double *ey, const double *ez, double m[3][3]);
 inline double det3(const double mat[3][3]);
@@ -763,6 +764,21 @@ inline void MathExtra::zeromat3(double **m)
   m[2][0] = m[2][1] = m[2][2] = 0.0;
 }
 
+// transpose a matrix
+
+inline void MathExtra::transpose3(const double m[3][3], double ans[3][3])
+{
+  ans[0][0] = m[0][0];
+  ans[0][1] = m[1][0];
+  ans[0][2] = m[2][0];
+  ans[1][0] = m[0][1];
+  ans[1][1] = m[1][1];
+  ans[1][2] = m[2][1];
+  ans[2][0] = m[0][2];
+  ans[2][1] = m[1][2];
+  ans[2][2] = m[2][2];
+}
+
 // add two matrices
 
 inline void MathExtra::plus3(const double m[3][3], double **m2, double **ans)
diff --git a/src/read_data.cpp b/src/read_data.cpp
index 4b4602f1fc..8df901062c 100644
--- a/src/read_data.cpp
+++ b/src/read_data.cpp
@@ -474,11 +474,14 @@ void ReadData::command(int narg, char **arg)
   int atomflag, topoflag;
   int bondflag, angleflag, dihedralflag, improperflag;
   int ellipsoidflag, lineflag, triflag, bodyflag;
-
+  
   atomflag = topoflag = 0;
   bondflag = angleflag = dihedralflag = improperflag = 0;
   ellipsoidflag = lineflag = triflag = bodyflag = 0;
 
+  xloxhi_flag = yloyhi_flag = zlozhi_flag = tilt_flag = 0;
+  avec_flag = bvec_flag = cvec_flag = 0;
+
   // values in this data file
 
   natoms = 0;
@@ -488,7 +491,12 @@ void ReadData::command(int narg, char **arg)
 
   boxlo[0] = boxlo[1] = boxlo[2] = -0.5;
   boxhi[0] = boxhi[1] = boxhi[2] = 0.5;
-  triclinic = 0;
+  avec[0] = avec[1] = avec[2] = 0.0;
+  bvec[0] = bvec[1] = bvec[2] = 0.0;
+  cvec[0] = cvec[1] = cvec[2] = 0.0;
+  avec[0] = bvec[1] = cvec[2] = 1.0;
+  tri_origin[0] = tri_origin[1] = tri_origin[2] = 0.0;
+  
   keyword[0] = '\0';
 
   nlocal_previous = atom->nlocal;
@@ -508,6 +516,17 @@ void ReadData::command(int narg, char **arg)
 
     header(firstpass);
 
+    // check if simulation box specified consistently
+
+    if (!avec_flag && !bvec_flag && !cvec_flag) {
+      triclinic = triclinic_general = 0;
+      if (tilt_flag) triclinic = 1;
+    } else {
+      if (xloxhi_flag || yloyhi_flag || zlozhi_flag || tilt_flag)
+        error->all(FLERR,"Read_data header cannot specify simulation box lo/hi/tilt and ABC vectors");
+      triclinic = triclinic_general = 1;
+    }
+    
     // problem setup using info from header
     // only done once, if firstpass and first data file
     // apply extra settings before grow(), even if no topology in file
@@ -536,31 +555,69 @@ void ReadData::command(int narg, char **arg)
       n = static_cast<int>(nbig);
       atom->avec->grow(n);
 
-      domain->boxlo[0] = boxlo[0];
-      domain->boxhi[0] = boxhi[0];
-      domain->boxlo[1] = boxlo[1];
-      domain->boxhi[1] = boxhi[1];
-      domain->boxlo[2] = boxlo[2];
-      domain->boxhi[2] = boxhi[2];
+      // setup simulation box
+      // 3 options: orthogonal, restricted triclinic, general triclinic
+      // for general triclinic: convect general ABC edge vectors to LAMMPS restricted triclinic
 
-      if (triclinic) {
-        domain->triclinic = 1;
-        domain->xy = xy;
-        domain->xz = xz;
-        domain->yz = yz;
+      if (!triclinic_general) {
+        domain->boxlo[0] = boxlo[0];
+        domain->boxhi[0] = boxhi[0];
+        domain->boxlo[1] = boxlo[1];
+        domain->boxhi[1] = boxhi[1];
+        domain->boxlo[2] = boxlo[2];
+        domain->boxhi[2] = boxhi[2];
+
+        if (triclinic) {
+          domain->triclinic = 1;
+          domain->xy = xy;
+          domain->xz = xz;
+          domain->yz = yz;
+        }
+
+      } else if (triclinic_general) {
+        domain->set_general_triclinic(avec,bvec,cvec,tri_origin);
       }
-
-      domain->print_box("  ");
-      domain->set_initial_box();
-      domain->set_global_box();
-      comm->set_proc_grid();
-      domain->set_local_box();
     }
 
     // change simulation box to be union of existing box and new box + shift
     // only done if firstpass and not first data file
+    // shift not allowed for general triclinic
 
     if (firstpass && addflag != NONE) {
+
+      // general triclinic
+      // first data file must also be general triclinic
+      // avec,bvec,vec must match first data file
+      // shift not allowed
+      
+      if (triclinic_general) {
+        if (!domain->triclinic_general)
+          error->all(FLERR,"Read_data subsequent file cannot switch to general triclinic");
+        int errflag = 0;
+        if (avec[0] != domain->avec[0] || avec[1] != domain->avec[1] || avec[2] != domain->avec[2])
+          errflag = 1;
+        if (bvec[0] != domain->bvec[0] || bvec[1] != domain->bvec[1] || bvec[2] != domain->bvec[2])
+          errflag = 1;
+        if (cvec[0] != domain->cvec[0] || cvec[1] != domain->cvec[1] || cvec[2] != domain->cvec[2])
+          errflag = 1;
+        if (tri_origin[0] != domain->tri_origin[0] || tri_origin[1] != domain->tri_origin[1] ||
+            tri_origin[2] != domain->tri_origin[2])
+          errflag = 1;
+        if (errflag)
+          error->all(FLERR,"Read_data subsequent file ABC vectors must be same as first file");
+        if (shift[0] != 0.0 || shift[1] != 0.0 || shift[2] != 0.0)
+          error->all(FLERR,"Read_data subsequent file with ABC vectors cannot define shift");
+
+      // restricted triclinic
+      // tilt factors must match first data file
+        
+      } else if (triclinic) {
+        if (!domain->triclinic || domain->triclinic_general)
+          error->all(FLERR,"Read_data subsequent file cannot switch to restricted triclinic");
+        if (xy != domain->xy || xz != domain->xz || yz != domain->yz)
+          error->all(FLERR,"Read_data subsequent file tilt factors must be same as first file");
+      }
+      
       double oldboxlo[3] = {domain->boxlo[0], domain->boxlo[1], domain->boxlo[2]};
       double oldboxhi[3] = {domain->boxhi[0], domain->boxhi[1], domain->boxhi[2]};
       domain->boxlo[0] = MIN(domain->boxlo[0], boxlo[0] + shift[0]);
@@ -570,7 +627,9 @@ void ReadData::command(int narg, char **arg)
       domain->boxlo[2] = MIN(domain->boxlo[2], boxlo[2] + shift[2]);
       domain->boxhi[2] = MAX(domain->boxhi[2], boxhi[2] + shift[2]);
 
-      // check of box has changed. If yes, warn about non-zero image flags
+      // check if box has changed
+      // if yes, warn about non-zero image flags
+      
       if ((oldboxlo[0] != domain->boxlo[0]) || (oldboxlo[1] != domain->boxlo[1]) ||
           (oldboxlo[2] != domain->boxlo[2]) || (oldboxhi[0] != domain->boxhi[0]) ||
           (oldboxhi[1] != domain->boxhi[1]) || (oldboxhi[2] != domain->boxhi[2])) {
@@ -588,19 +647,16 @@ void ReadData::command(int narg, char **arg)
         if ((flag_all > 0) && (comm->me == 0))
           error->warning(FLERR, "Non-zero image flags with growing box leads to bad coordinates");
       }
-
-      // NOTE: not sure what to do about tilt value in subsequent data files
-      //if (triclinic) {
-      //  domain->xy = xy; domain->xz = xz; domain->yz = yz;
-      // }
-
-      domain->print_box("  ");
-      domain->set_initial_box();
-      domain->set_global_box();
-      comm->set_proc_grid();
-      domain->set_local_box();
     }
 
+    // setup simulation box and paritioning in Domain and Comm classes
+    
+    domain->print_box("  ");
+    domain->set_initial_box();
+    domain->set_global_box();
+    comm->set_proc_grid();
+    domain->set_local_box();
+
     // allocate space for type label map
 
     if (firstpass) {
@@ -608,8 +664,9 @@ void ReadData::command(int narg, char **arg)
       lmap = new LabelMap(lmp, ntypes, nbondtypes, nangletypes, ndihedraltypes, nimpropertypes);
     }
 
+    // rest of data file is Sections
+    // read in any order, except where error checks
     // customize for new sections
-    // read rest of file in free format
 
     while (strlen(keyword)) {
 
@@ -1140,7 +1197,8 @@ void ReadData::header(int firstpass)
 
   // initialize type counts by the "extra" numbers so they get counted
   // in case the corresponding "types" line is missing and thus the extra
-  // value will not be processed.
+  // value will not be processed
+  
   if (addflag == NONE) {
     atom->ntypes = extra_atom_types;
     atom->nbondtypes = extra_bond_types;
@@ -1156,6 +1214,7 @@ void ReadData::header(int firstpass)
     if (eof == nullptr) error->one(FLERR, "Unexpected end of data file");
 
     // check for units keyword in first line and print warning on mismatch
+    
     auto units = Tokenizer(utils::strfind(line, "units = \\w+")).as_vector();
     if (units.size() > 2) {
       if (units[2] != update->unit_style)
@@ -1278,7 +1337,7 @@ void ReadData::header(int firstpass)
       else if (firstpass)
         atom->nimpropers += nimpropers;
 
-      // Atom class type settings are only set by first data file
+    // Atom class type settings are only set by first data file
 
     } else if (utils::strmatch(line, "^\\s*\\d+\\s+atom\\s+types\\s")) {
       ntypes = utils::inumeric(FLERR, words[0], false, lmp);
@@ -1300,11 +1359,11 @@ void ReadData::header(int firstpass)
       nimpropertypes = utils::inumeric(FLERR, words[0], false, lmp);
       if (addflag == NONE) atom->nimpropertypes = nimpropertypes + extra_improper_types;
 
-      // these settings only used by first data file
-      // also, these are obsolescent. we parse them to maintain backward
-      // compatibility, but the recommended way is to set them via keywords
-      // in the LAMMPS input file. In case these flags are set in both,
-      // the input and the data file, we use the larger of the two.
+    // these settings only used by first data file
+    // also, these are obsolescent. we parse them to maintain backward
+    // compatibility, but the recommended way is to set them via keywords
+    // in the LAMMPS input file. In case these flags are set in both,
+    // the input and the data file, we use the larger of the two.
 
     } else if (strstr(line, "extra bond per atom")) {
       if (addflag == NONE) extra_flag_value = utils::inumeric(FLERR, words[0], false, lmp);
@@ -1322,26 +1381,50 @@ void ReadData::header(int firstpass)
       if (addflag == NONE) extra_flag_value = utils::inumeric(FLERR, words[0], false, lmp);
       force->special_extra = MAX(force->special_extra, extra_flag_value);
 
-      // local copy of box info
-      // so can treat differently for first vs subsequent data files
+    // local copy of box info
+    // so can treat differently for first vs subsequent data files
 
     } else if (utils::strmatch(line, "^\\s*\\f+\\s+\\f+\\s+xlo\\s+xhi\\s")) {
+      xloxhi_flag = 1;
       boxlo[0] = utils::numeric(FLERR, words[0], false, lmp);
       boxhi[0] = utils::numeric(FLERR, words[1], false, lmp);
 
     } else if (utils::strmatch(line, "^\\s*\\f+\\s+\\f+\\s+ylo\\s+yhi\\s")) {
+      yloyhi_flag = 1;
       boxlo[1] = utils::numeric(FLERR, words[0], false, lmp);
       boxhi[1] = utils::numeric(FLERR, words[1], false, lmp);
 
     } else if (utils::strmatch(line, "^\\s*\\f+\\s+\\f+\\s+zlo\\s+zhi\\s")) {
+      zlozhi_flag = 1;
       boxlo[2] = utils::numeric(FLERR, words[0], false, lmp);
       boxhi[2] = utils::numeric(FLERR, words[1], false, lmp);
 
     } else if (utils::strmatch(line, "^\\s*\\f+\\s+\\f+\\s+\\f+\\s+xy\\s+xz\\s+yz\\s")) {
-      triclinic = 1;
+      tilt_flag = 1;
       xy = utils::numeric(FLERR, words[0], false, lmp);
       xz = utils::numeric(FLERR, words[1], false, lmp);
       yz = utils::numeric(FLERR, words[2], false, lmp);
+      
+    } else if (utils::strmatch(line, "^\\s*\\f+\\s+\\f+\\s+\\f+\\s+\\f+\\s+\\avec\\s")) {
+      avec_flag = 1;
+      avec[0] = utils::numeric(FLERR, words[0], false, lmp);
+      avec[1] = utils::numeric(FLERR, words[1], false, lmp);
+      avec[2] = utils::numeric(FLERR, words[2], false, lmp);
+      tri_origin[0] = utils::numeric(FLERR, words[3], false, lmp);
+
+    } else if (utils::strmatch(line, "^\\s*\\f+\\s+\\f+\\s+\\f+\\s+\\f+\\s+\\bvec\\s")) {
+      bvec_flag = 1;
+      bvec[0] = utils::numeric(FLERR, words[0], false, lmp);
+      bvec[1] = utils::numeric(FLERR, words[1], false, lmp);
+      bvec[2] = utils::numeric(FLERR, words[2], false, lmp);
+      tri_origin[1] = utils::numeric(FLERR, words[3], false, lmp);
+
+    } else if (utils::strmatch(line, "^\\s*\\f+\\s+\\f+\\s+\\f+\\s+\\f+\\s+\\cvec\\s")) {
+      cvec_flag = 1;
+      cvec[0] = utils::numeric(FLERR, words[0], false, lmp);
+      cvec[1] = utils::numeric(FLERR, words[1], false, lmp);
+      cvec[2] = utils::numeric(FLERR, words[2], false, lmp);
+      tri_origin[2] = utils::numeric(FLERR, words[3], false, lmp);
 
     } else
       break;
@@ -1407,8 +1490,8 @@ void ReadData::atoms()
     if (eof) error->all(FLERR, "Unexpected end of data file");
     if (tlabelflag && !lmap->is_complete(Atom::ATOM))
       error->all(FLERR, "Label map is incomplete: all types must be assigned a unique type label");
-    atom->data_atoms(nchunk, buffer, id_offset, mol_offset, toffset, shiftflag, shift, tlabelflag,
-                     lmap->lmap2lmap.atom);
+    atom->data_atoms(nchunk, buffer, id_offset, mol_offset, toffset, triclinic_general,
+                     shiftflag, shift, tlabelflag, lmap->lmap2lmap.atom);
     nread += nchunk;
   }
 
diff --git a/src/read_data.h b/src/read_data.h
index e97f225997..6439fbd096 100644
--- a/src/read_data.h
+++ b/src/read_data.h
@@ -63,8 +63,11 @@ class ReadData : public Command {
 
   double boxlo[3], boxhi[3];
   double xy, xz, yz;
-  int triclinic;
-
+  double avec[3], bvec[3], cvec[3], tri_origin[3];
+  int triclinic, triclinic_general;
+  int xloxhi_flag, yloyhi_flag, zlozhi_flag, tilt_flag;
+  int avec_flag, bvec_flag, cvec_flag;
+  
   // optional args
 
   int addflag, offsetflag, shiftflag, coeffflag, settypeflag;
diff --git a/src/write_data.cpp b/src/write_data.cpp
index dd5b056ae8..0ab72f857a 100644
--- a/src/write_data.cpp
+++ b/src/write_data.cpp
@@ -72,8 +72,11 @@ void WriteData::command(int narg, char **arg)
   pairflag = II;
   coeffflag = 1;
   fixflag = 1;
+  triclinic_general = 0;
   lmapflag = 1;
-  // store current (default) setting since we may change it.
+  
+  // store current (default) setting since we may change it
+  
   int types_style = atom->types_style;
   int noinit = 0;
 
@@ -94,6 +97,9 @@ void WriteData::command(int narg, char **arg)
     } else if (strcmp(arg[iarg],"nofix") == 0) {
       fixflag = 0;
       iarg++;
+    } else if (strcmp(arg[iarg],"triclinic") == 0) {
+      triclinic_general = 1;
+      iarg++;
     } else if (strcmp(arg[iarg],"nolabelmap") == 0) {
       lmapflag = 0;
       iarg++;
@@ -135,7 +141,9 @@ void WriteData::command(int narg, char **arg)
   }
 
   write(file);
+
   // restore saved setting
+  
   atom->types_style = types_style;
 }
 
@@ -206,8 +214,32 @@ void WriteData::write(const std::string &file)
   }
 
   // per atom info in Atoms and Velocities sections
-
+  // if general triclinic:
+  //   save restricted triclinic atom coords
+  //   transform atom coords from restricted to general
+  //   restore save atom coords after output
+  // NOTE: do same for velocities as well ?
+  
+  double **xstore = nullptr;
+  
+  if (triclinic_general) {
+    double **x = atom->x;
+    int nlocal = atom->nlocal;
+    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]);
+  }
+  
   if (natoms) atoms();
+
+  if (triclinic_general) {
+    double **x = atom->x;
+    int nlocal = atom->nlocal;
+    if (nlocal) memcpy(&x[0][0],&xstore[0][0],3*nlocal*sizeof(double));
+    memory->destroy(xstore);
+  }
+
   if (natoms) velocities();
 
   // molecular topology info if defined
@@ -289,15 +321,22 @@ void WriteData::header()
         for (int m = 0; m < ifix->wd_header; m++)
           ifix->write_data_header(fp,m);
 
-  // box info
+  // box info: orthogonal, restricted triclinic, or general triclinic (if requested)
 
-  auto box = fmt::format("\n{} {} xlo xhi\n{} {} ylo yhi\n{} {} zlo zhi\n",
-                         domain->boxlo[0],domain->boxhi[0],
-                         domain->boxlo[1],domain->boxhi[1],
-                         domain->boxlo[2],domain->boxhi[2]);
-  if (domain->triclinic)
-    box += fmt::format("{} {} {} xy xz yz\n",domain->xy,domain->xz,domain->yz);
-  fputs(box.c_str(),fp);
+  if (!triclinic_general) {
+    fmt::print(fp,"\n{} {} xlo xhi\n{} {} ylo yhi\n{} {} zlo zhi\n",
+               domain->boxlo[0],domain->boxhi[0],
+               domain->boxlo[1],domain->boxhi[1],
+               domain->boxlo[2],domain->boxhi[2]);
+    if (domain->triclinic)
+      fmt::print(fp,"{} {} {} xy xz yz\n",domain->xy,domain->xz,domain->yz);
+    
+  } else if (triclinic_general) {
+    fmt::print(fp,"\n{} {} {} {} avec\n{} {} {} {} bvec\n{} {} {} {} cvec\n",
+               domain->avec[0],domain->avec[1],domain->avec[2],domain->tri_origin[0],
+               domain->bvec[0],domain->bvec[1],domain->bvec[2],domain->tri_origin[1],
+               domain->cvec[0],domain->cvec[1],domain->cvec[2],domain->tri_origin[2]);
+  }
 }
 
 /* ----------------------------------------------------------------------
diff --git a/src/write_data.h b/src/write_data.h
index f0df9b0c5f..ebaa97ffc0 100644
--- a/src/write_data.h
+++ b/src/write_data.h
@@ -35,6 +35,7 @@ class WriteData : public Command {
   int pairflag;
   int coeffflag;
   int fixflag;
+  int triclinic_general;
   int lmapflag;
   FILE *fp;
   bigint nbonds_local, nbonds;