more updates to doc page and read_data functionality

doc/src/Howto_body.rst

@@ -102,8 +102,8 @@ particles of different styles
 | :doc:`dump image <dump_image>`                 | output body particle attributes as an image         |
 +------------------------------------------------+-----------------------------------------------------+
 
-The pair styles defined for use with specific body styles are listed
+The pair styles currently defined for use with specific body styles
-in the sections below.
+are listed in the sections below.
 
 ----------
 

doc/src/read_data.rst

@@ -914,7 +914,7 @@ zero.
    If the data file defines a general triclinic box, then the
    following per-atom values in the list above are per-atom vectors
    which imply an orientation: (mux,muy,muz), (normx,normy,normz),
-   (spx,spy,spz).  This menas they should be specified consistent with
+   (spx,spy,spz).  This means they should be specified consistent with
    the general triclinic box and its orientation relative to the
    standard x,y,z coordinate axes.  For example a dipole moment vector
    which will be in the +x direction once LAMMPS converts from a
@@ -1352,8 +1352,8 @@ specified as a unit vector.
 If the data file defines a general triclinic box, then the quaternion
 for each ellipsoid should be specified for its orientation relative to
 the standard x,y,z coordinate axes.  When the system is converted to a
-restricted triclinic box, the ellipsoid quaternions will be altered
+restricted triclinic box, the ellipsoid quaternions will be altered to
-appropriately.
+reflect the new orientation of the ellipsoid.
 
 The *Ellipsoids* section must appear after the *Atoms* section.
 
@@ -1471,24 +1471,23 @@ is used and any atoms are listed in the *Atoms* section with a
 lineflag = 1.  The number of lines should be specified in the header
 section via the "lines" keyword.
 
-The 2 end points are the end points of the line segment.  The ordering
+The 2 end points are the end points of the line segment.  They should
-of the 2 points should be such that using a right-hand rule to cross
+be values close to the center point of the line segment specified in
-the line segment with a unit vector in the +z direction, gives an
+the Atoms section of the data file, even if individual end points are
-"outward" normal vector perpendicular to the line segment.
+outside the simulation box.
+
+The ordering of the 2 points should be such that using a right-hand
+rule to cross the line segment with a unit vector in the +z direction,
+gives an "outward" normal vector perpendicular to the line segment.
 I.e. normal = (c2-c1) x (0,0,1).  This orientation may be important
 for defining some interactions.
 
-If the data file defines a general triclinic box, then the quaternion
+If the data file defines a general triclinic box, then the x1,y1 and
-for each ellipsoid should be specified for its orientation relative to
+x2,y2 values for each line segment should be specified for its
-the standard x,y,z coordinate axes.  When the system is converted to a
+orientation relative to the standard x,y,z coordinate axes.  When the
-restricted triclinic box, the ellipsoid quaternions will be altered
+system is converted to a restricted triclinic box, the x1,y1,x2,y2
-appropriately.
+values will be altered to reflect the new orientation of the line
-
+segment.
-If the data file defines a general triclinic box, the (x1,y1) and
-(x2,y2) values should be within (or near) its parallelogram area,
-i.e. near the x,y coordinates of the line segment as defined in the
-Atoms section.  See the :doc:`Howto triclinic <Howto_triclininc>` doc
-page for more details.
 
 The *Lines* section must appear after the *Atoms* section.
 
@@ -1616,17 +1615,22 @@ The *Triangles* section must appear if :doc:`atom_style tri
 with a triangleflag = 1.  The number of lines should be specified in
 the header section via the "triangles" keyword.
 
-The 3 corner points are the corner points of the triangle.  The
+The 3 corner points are the corner points of the triangle.  They
-ordering of the 3 points should be such that using a right-hand rule
+should be values close to the center point of the triangle specified
-to go from point1 to point2 to point3 gives an "outward" normal vector
+in the Atoms section of the data file, even if individual corner
-to the face of the triangle.  I.e. normal = (c2-c1) x (c3-c1).  This
+points are outside the simulation box.
-orientation may be important for defining some interactions.
 
-If the data file defines a general triclinic box, the (x1,y1,z1),
+The ordering of the 3 points should be such that using a right-hand
-(x2,y2,z2), (x3,y3,z3) corner points should be within (or near) its
+rule to go from point1 to point2 to point3 gives an "outward" normal
-parallelepiped volume, i.e. near the x,y,z coordinates of the triangle
+vector to the face of the triangle.  I.e. normal = (c2-c1) x (c3-c1).
-as defined in the Atoms section.  See the :doc:`Howto triclinic
+This orientation may be important for defining some interactions.
-<Howto_triclininc>` doc page for more details.
+
+If the data file defines a general triclinic box, then the x1,y1,z1
+and x2,y2,z2 and x3,y3,z3 values for each triangle should be specified
+for its orientation relative to the standard x,y,z coordinate axes.
+When the system is converted to a restricted triclinic box, the
+x1,y1,z1,x2,y2,z2,x3,y3,z3 values will be altered to reflect the new
+orientation of the triangle.
 
 The *Triangles* section must appear after the *Atoms* section.
 

src/DIELECTRIC/atom_vec_dielectric.cpp

@@ -208,7 +208,7 @@ void AtomVecDielectric::read_data_general_to_restricted(int nlocal_previous, int
 /* ----------------------------------------------------------------------
    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 dipole momemt mu
+   child class operates on dipole momemt mu which has 4 values per atom
 ------------------------------------------------------------------------- */
 
 void AtomVecDielectric::write_data_restricted_to_general()
@@ -217,16 +217,17 @@ void AtomVecDielectric::write_data_restricted_to_general()
 
   int nlocal = atom->nlocal;
   memory->create(mu_hold,nlocal,3,"atomvec:mu_hold");
-  if (nlocal) memcpy(&mu_hold[0][0],&mu[0][0],3*nlocal*sizeof(double));
+    for (int i = 0; i < nlocal; i++) {
-  for (int i = 0; i < nlocal; i++)
+    memcpy(&mu_hold[i],&mu[i],3*sizeof(double));
     domain->restricted_to_general_vector(mu[i]);
+  }
 }
 
 /* ----------------------------------------------------------------------
    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 dipole moment mu
+   child class operates on dipole moment mu which has 4 values per atom
 ------------------------------------------------------------------------- */
 
 void AtomVecDielectric::write_data_restore_restricted()
@@ -236,7 +237,8 @@ void AtomVecDielectric::write_data_restore_restricted()
   if (!mu_hold) return;
   
   int nlocal = atom->nlocal;
-  memcpy(&mu[0][0],&mu_hold[0][0],3*nlocal*sizeof(double));
+  for (int i = 0; i < nlocal; i++)
+    memcpy(&mu[i],&mu_hold[i],3*sizeof(double));
   memory->destroy(mu_hold);
   mu_hold = nullptr;
 }

src/DIPOLE/atom_vec_dipole.cpp

@@ -90,7 +90,7 @@ void AtomVecDipole::read_data_general_to_restricted(int nlocal_previous, int nlo
 /* ----------------------------------------------------------------------
    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 dipole momemt mu
+   child class operates on dipole momemt mu which has 4 values per atom
 ------------------------------------------------------------------------- */
 
 void AtomVecDipole::write_data_restricted_to_general()
@@ -99,16 +99,17 @@ void AtomVecDipole::write_data_restricted_to_general()
 
   int nlocal = atom->nlocal;
   memory->create(mu_hold,nlocal,3,"atomvec:mu_hold");
-  if (nlocal) memcpy(&mu_hold[0][0],&mu[0][0],3*nlocal*sizeof(double));
+  for (int i = 0; i < nlocal; i++) {
-  for (int i = 0; i < nlocal; i++)
+    memcpy(&mu_hold[i],&mu[i],3*sizeof(double));
     domain->restricted_to_general_vector(mu[i]);
+  }
 }
 
 /* ----------------------------------------------------------------------
    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 dipole moment mu
+   child class operates on dipole moment mu which has 4 values per atom
 ------------------------------------------------------------------------- */
 
 void AtomVecDipole::write_data_restore_restricted()
@@ -118,7 +119,8 @@ void AtomVecDipole::write_data_restore_restricted()
   if (!mu_hold) return;
   
   int nlocal = atom->nlocal;
-  memcpy(&mu[0][0],&mu_hold[0][0],3*nlocal*sizeof(double));
+  for (int i = 0; i < nlocal; i++)
+    memcpy(&mu[i],&mu_hold[i],3*sizeof(double));
   memory->destroy(mu_hold);
   mu_hold = nullptr;
 }

src/MACHDYN/atom_vec_smd.cpp

@@ -164,7 +164,7 @@ void AtomVecSMD::data_atom_post(int ilocal)
   // x and x0 are in Atoms section of data file
   // reset x0 b/c x may have been modified in Atom::data_atoms()
   //   for PBC, shift, etc
-  // this also means no need for read_data_general_to_restricted() method
+  // this means no need for read_data_general_to_restricted() method
   //   to rotate x0 for general triclinic
   
   x0[ilocal][0] = x[ilocal][0];

src/SPIN/atom_vec_spin.cpp

@@ -120,7 +120,7 @@ void AtomVecSpin::read_data_general_to_restricted(int nlocal_previous, int nloca
 /* ----------------------------------------------------------------------
    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 spin vector sp
+   child class operates on spin vector sp which has 4 values per atom
 ------------------------------------------------------------------------- */
 
 void AtomVecSpin::write_data_restricted_to_general()
@@ -129,16 +129,17 @@ void AtomVecSpin::write_data_restricted_to_general()
 
   int nlocal = atom->nlocal;
   memory->create(sp_hold,nlocal,3,"atomvec:sp_hold");
-  if (nlocal) memcpy(&sp_hold[0][0],&sp[0][0],3*nlocal*sizeof(double));
+  for (int i = 0; i < nlocal; i++) {
-  for (int i = 0; i < nlocal; i++)
+    memcpy(&sp_hold[i],&sp[i],3*sizeof(double));
     domain->restricted_to_general_vector(sp[i]);
+  }
 }
 
 /* ----------------------------------------------------------------------
    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 spin vector sp
+   child class operates on spin vector sp which has 4 values per atom
 ------------------------------------------------------------------------- */
 
 void AtomVecSpin::write_data_restore_restricted()
@@ -148,7 +149,8 @@ void AtomVecSpin::write_data_restore_restricted()
   if (!sp_hold) return;
   
   int nlocal = atom->nlocal;
-  memcpy(&sp[0][0],&sp_hold[0][0],3*nlocal*sizeof(double));
+  for (int i = 0; i < nlocal; i++)
+    memcpy(&sp[i],&sp_hold[i],3*sizeof(double));
   memory->destroy(sp_hold);
   sp_hold = nullptr;
 }

src/atom_vec_line.cpp

@@ -355,9 +355,10 @@ void AtomVecLine::data_atom_bonus(int m, const std::vector<std::string> &values)
 
   // convert x1/y1 and x2/y2 from general to restricted triclniic
   // x is already restricted triclinic
-  
+
+  double coords[3];
+
   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];
@@ -366,8 +367,20 @@ void AtomVecLine::data_atom_bonus(int m, const std::vector<std::string> &values)
     x2 = coords[0]; y2 = coords[1];
   }
 
+  // remap end points to be near x
+  // necessary if atom x was remapped into periodic box
+
+  coords[0] = x1; coords[1] = y1; coords[2] = 0.0;
+  domain->remap_near(coords,x[m]);
+  x1 = coords[0]; y1 = coords[1];
+  coords[0] = x2; coords[1] = y2; coords[2] = 0.0;
+  domain->remap_near(c2,x[m]);
+  x2 = coords[0]; y2 = coords[1];
+
   // calculate length and theta
-  
+  // error if segment center is not within EPSILON of atom x
+  // reset atom x to center point
+
   double dx = x2 - x1;
   double dy = y2 - y1;
   double length = sqrt(dx * dx + dy * dy);

src/atom_vec_line.h

@@ -76,7 +76,6 @@ class AtomVecLine : public AtomVec {
 
   void grow_bonus();
   void copy_bonus_all(int, int);
-  // void consistency_check(int, char *);
 };
 
 }    // namespace LAMMPS_NS

src/atom_vec_tri.cpp

@@ -52,8 +52,6 @@ AtomVecTri::AtomVecTri(LAMMPS *lmp) : AtomVec(lmp)
   nlocal_bonus = nghost_bonus = nmax_bonus = 0;
   bonus = nullptr;
 
-  double **quat_hold = nullptr;
-  
   // strings with peratom variables to include in each AtomVec method
   // strings cannot contain fields in corresponding AtomVec default strings
   // order of fields in a string does not matter
@@ -518,9 +516,17 @@ void AtomVecTri::data_atom_bonus(int m, const std::vector<std::string> &values)
     domain->general_to_restricted_coords(c2);
     domain->general_to_restricted_coords(c3);
   }
+
+  // remap corner points to be near x
+  // necessary if atom x was remapped into periodic box
+
+  domain->remap_near(c1,x[m]);
+  domain->remap_near(c2,x[m]);
+  domain->remap_near(c3,x[m]);
   
   // centroid = 1/3 of sum of vertices
-  // error if centroid is not within EPSILON of Atoms section coord
+  // error if centroid is not within EPSILON of atom x
+  // reset atom x to centroid
   
   double centroid[3];
   centroid[0] = (c1[0] + c2[0] + c3[0]) / 3.0;

src/atom_vec_tri.h

@@ -72,8 +72,6 @@ class AtomVecTri : public AtomVec {
   double *radius, *rmass;
   double **omega, **angmom;
 
-  double **quat_hold;
-
   int nghost_bonus, nmax_bonus;
   int tri_flag;
   double rmass_one;

src/domain.cpp

@@ -1361,7 +1361,7 @@ int Domain::closest_image(const double * const pos, int j)
 /* ----------------------------------------------------------------------
    find and return Xj image = periodic image of Xj that is closest to Xi
    for triclinic, add/subtract tilt factors in other dims as needed
-   called by ServerMD class and LammpsInterface in lib/atc.
+   called by ServerMD class and LammpsInterface in lib/atc
 ------------------------------------------------------------------------- */
 
 void Domain::closest_image(const double * const xi, const double * const xj, double * const xjimage)

src/write_data.cpp

@@ -98,6 +98,9 @@ void WriteData::command(int narg, char **arg)
       fixflag = 0;
       iarg++;
     } else if (strcmp(arg[iarg],"triclinic/general") == 0) {
+      if (!domain->triclinic_general)
+        error->all(FLERR,"Write_data triclinic/general for system "
+                   "that is not general triclinic");
       triclinic_general = 1;
       iarg++;
     } else if (strcmp(arg[iarg],"nolabelmap") == 0) {
@@ -213,10 +216,10 @@ void WriteData::write(const std::string &file)
     if (coeffflag) force_fields();
   }
 
-  // if general triclinic:
+  // if general triclinic requested:
   // reset internal per-atom data that needs rotation
 
-  atom->avec->write_data_restricted_to_general();
+  if (triclinic_general) atom->avec->write_data_restricted_to_general();
   
   // per atom info in Atoms and Velocities sections
   
@@ -247,10 +250,10 @@ void WriteData::write(const std::string &file)
       if (ifix->wd_section)
         for (int m = 0; m < ifix->wd_section; m++) fix(ifix,m);
 
-  // if general triclinic:
+  // if general triclinic requested:
   // restore internal per-atom data that was rotated
 
-  atom->avec->write_data_restore_restricted();
+  if (triclinic_general) atom->avec->write_data_restore_restricted();
   
   // close data file