finalized dump command support for general triclinic

doc/src/create_atoms.rst

@@ -112,6 +112,21 @@ added to the specified atom *type* (e.g., if *type* = 2 and the file
 specifies atom types 1, 2, and 3, then each created molecule will have
 atom types 3, 4, and 5).
 
+.. note::
+
+   You cannot use this command to create atoms that are outside the
+   simulation box; they will just be ignored by LAMMPS.  This is true
+   even if you are using shrink-wrapped box boundaries, as specified
+   by the :doc:`boundary <boundary>` command.  However, you can first
+   use the :doc:`change_box <change_box>` command to temporarily
+   expand the box, then add atoms via create_atoms, then finally use
+   change_box command again if needed to re-shrink-wrap the new atoms.
+   See the :doc:`change_box <change_box>` doc page for an example of
+   how to do this, using the create_atoms *single* style to insert a
+   new atom outside the current simulation box.
+
+----------
+
 For the *box* style, the create_atoms command fills the entire
 simulation box with particles on the lattice.  If your simulation box
 is periodic, you should ensure its size is a multiple of the lattice
@@ -182,7 +197,64 @@ styles should be made in that context.
 For the *single* style, a single particle is added to the system at
 the specified coordinates.  This can be useful for debugging purposes
 or to create a tiny system with a handful of particles at specified
-positions.
+positions.  For a 2d simulation the specified z coordinate must be
+0.0.
+
+.. versionchanged:: 2Jun2022
+
+The *porosity* style has been renamed to *random* with added functionality.
+
+For the *random* style, *N* particles are added to the system at
+randomly generated coordinates, which can be useful for generating an
+amorphous system.  For 2d simulations, the z coordinates of all added
+atoms will be 0.0.
+
+The particles are created one by one using the specified random number
+*seed*, resulting in the same set of particle coordinates, independent
+of how many processors are being used in the simulation.  Unless the
+*overlap* keyword is specified, particles created by the *random*
+style will typically be highly overlapped.  Various additional
+criteria can be used to accept or reject a random particle insertion;
+see the keyword discussion below.  Multiple attempts per particle are
+made (see the *maxtry* keyword) until the insertion is either
+successful or fails.  If this command fails to add all requested *N*
+particles, a warning will be output.
+
+If the *region-ID* argument is specified as NULL, then the randomly
+created particles will be anywhere in the simulation box.  If a
+*region-ID* is specified, a geometric volume is filled that is both
+inside the simulation box and is also consistent with the region
+volume.  See the :doc:`region <region>` command for details.  Note
+that a region can be specified so that its "volume" is either inside
+or outside its geometric boundary.
+
+Note that the create_atoms command adds particles to those that
+already exist.  This means it can be used to add particles to a system
+previously read in from a data or restart file.  Or the create_atoms
+command can be used multiple times, to add multiple sets of particles
+to the simulation.  For example, grain boundaries can be created, by
+interleaving the create_atoms command with :doc:`lattice <lattice>`
+commands specifying different orientations.
+
+When this command is used, care should be taken to ensure the
+resulting system does not contain particles that are highly
+overlapped.  Such overlaps will cause many interatomic potentials to
+compute huge energies and forces, leading to bad dynamics.  There are
+several strategies to avoid this problem:
+
+* Use the :doc:`delete_atoms overlap <delete_atoms>` command after
+  create_atoms.  For example, this strategy can be used to overlay and
+  surround a large protein molecule with a volume of water molecules,
+  then delete water molecules that overlap with the protein atoms.
+
+* For the *random* style, use the optional *overlap* keyword to avoid
+  overlaps when each new particle is created.
+
+* Before running dynamics on an overlapped system, perform an
+  :doc:`energy minimization <minimize>`.  Or run initial dynamics with
+  :doc:`pair_style soft <pair_soft>` or with :doc:`fix nve/limit
+  <fix_nve_limit>` to un-overlap the particles, before running normal
+  dynamics.
 
 .. figure:: img/marble_race.jpg
             :figwidth: 33%
@@ -242,73 +314,6 @@ to the area of that triangle.
    beneficial to exclude computing interactions between the created
    particles using :doc:`neigh_modify exclude <neigh_modify>`.
 
-.. versionchanged:: 2Jun2022
-
-The *porosity* style has been renamed to *random* with added functionality.
-
-For the *random* style, *N* particles are added to the system at
-randomly generated coordinates, which can be useful for generating an
-amorphous system.  The particles are created one by one using the
-specified random number *seed*, resulting in the same set of particle
-coordinates, independent of how many processors are being used in the
-simulation.  Unless the *overlap* keyword is specified, particles
-created by the *random* style will typically be highly overlapped.
-Various additional criteria can be used to accept or reject a random
-particle insertion; see the keyword discussion below.  Multiple
-attempts per particle are made (see the *maxtry* keyword) until the
-insertion is either successful or fails.  If this command fails to add
-all requested *N* particles, a warning will be output.
-
-If the *region-ID* argument is specified as NULL, then the randomly
-created particles will be anywhere in the simulation box.  If a
-*region-ID* is specified, a geometric volume is filled that is both
-inside the simulation box and is also consistent with the region
-volume.  See the :doc:`region <region>` command for details.  Note
-that a region can be specified so that its "volume" is either inside
-or outside its geometric boundary.
-
-Note that the create_atoms command adds particles to those that
-already exist.  This means it can be used to add particles to a system
-previously read in from a data or restart file.  Or the create_atoms
-command can be used multiple times, to add multiple sets of particles
-to the simulation.  For example, grain boundaries can be created, by
-interleaving the create_atoms command with :doc:`lattice <lattice>`
-commands specifying different orientations.
-
-When this command is used, care should be taken to ensure the
-resulting system does not contain particles that are highly
-overlapped.  Such overlaps will cause many interatomic potentials to
-compute huge energies and forces, leading to bad dynamics.  There are
-several strategies to avoid this problem:
-
-* Use the :doc:`delete_atoms overlap <delete_atoms>` command after
-  create_atoms.  For example, this strategy can be used to overlay and
-  surround a large protein molecule with a volume of water molecules,
-  then delete water molecules that overlap with the protein atoms.
-
-* For the *random* style, use the optional *overlap* keyword to avoid
-  overlaps when each new particle is created.
-
-* Before running dynamics on an overlapped system, perform an
-  :doc:`energy minimization <minimize>`.  Or run initial dynamics with
-  :doc:`pair_style soft <pair_soft>` or with :doc:`fix nve/limit
-  <fix_nve_limit>` to un-overlap the particles, before running normal
-  dynamics.
-
-.. note::
-
-   You cannot use any of the styles explained above to create atoms
-   that are outside the simulation box; they will just be ignored by
-   LAMMPS.  This is true even if you are using shrink-wrapped box
-   boundaries, as specified by the :doc:`boundary <boundary>` command.
-   However, you can first use the :doc:`change_box <change_box>`
-   command to temporarily expand the box, then add atoms via
-   create_atoms, then finally use change_box command again if needed
-   to re-shrink-wrap the new atoms.  See the :doc:`change_box
-   <change_box>` doc page for an example of how to do this, using the
-   create_atoms *single* style to insert a new atom outside the
-   current simulation box.
-
 ----------
 
 Individual atoms are inserted by this command, unless the *mol*

doc/src/dump.rst

@@ -278,16 +278,20 @@ format <dump_modify>` command and its options.
 Format of native LAMMPS format dump files:
 
 The *atom*, *custom*, *grid*, and *local* styles create files in a
-simple LAMMPS-specific text format that is self-explanatory when
+simple LAMMPS-specific text format that is mostly self-explanatory
-viewing a dump file.  Many post-processing tools either included with
+when viewing a dump file.  Many post-processing tools either included
-LAMMPS or third-party tools can read this format, as does the
+with LAMMPS or third-party tools can read this format, as does the
 :doc:`rerun <rerun>` command.  See tools described on the :doc:`Tools
 <Tools>` doc page for examples, including `Pizza.py
 <https://lammps.github.io/pizza>`_.
 
 For all these styles, the dimensions of the simulation box are
-included in each snapshot.  For an orthogonal simulation box this
+included in each snapshot.  The simulation box in LAMMPS can be
-information is formatted as:
+defined in one of 3 ways: orthogonal, restricted triclinic, and
+general triclinic.  See the :doc:`Howto triclinic <Howto_triclininc>`
+doc page for a detailed description of all 3 options.
+
+For an orthogonal simulation box the box information is formatted as:
 
 .. parsed-literal::
 
@@ -304,10 +308,10 @@ the six characters is one of *p* (periodic), *f* (fixed), *s* (shrink wrap),
 or *m* (shrink wrapped with a minimum value).  See the
 :doc:`boundary <boundary>` command for details.
 
-For triclinic simulation boxes (non-orthogonal), an orthogonal
+For a restricted triclinic simulation box, an orthogonal bounding box
-bounding box which encloses the triclinic simulation box is output,
+which encloses the restricted triclinic simulation box is output,
-along with the three tilt factors (*xy*, *xz*, *yz*) of the triclinic box,
+along with the three tilt factors (*xy*, *xz*, *yz*) of the triclinic
-formatted as follows:
+box, formatted as follows:
 
 .. parsed-literal::
 
@@ -329,6 +333,10 @@ bounding box extents (xlo_bound, xhi_bound, etc.) are calculated from the
 triclinic parameters, and how to transform those parameters to and
 from other commonly used triclinic representations.
 
+For a general triclinic simulation box, see the "General triclinic"
+section below for a description of the ITEM: BOX BOUNDS format as well
+as how per-atom coordinates and per-atom vector quantities are output.
+
 The *atom* and *custom* styles output a "ITEM: NUMBER OF ATOMS" line
 with the count of atoms in the snapshot.  Likewise they output an
 "ITEM: ATOMS" line which includes column descriptors for the per-atom
@@ -400,7 +408,6 @@ command.
 
 Dump files in other popular formats:
 
-
 .. note::
 
    This section only discusses file formats relevant to this doc page.
@@ -656,6 +663,87 @@ how to control the compression level in both variants.
 
 ----------
 
+General triclinic simulation box output for the *atom* and *custom* styles:
+
+As mentioned above, the simulation box can be defined as a general
+triclinic box, which means that 3 arbitrary box edge vectors **A**,
+**B**, **C** can be specified.  See the :doc:`Howto triclinic
+<Howto_triclininc>` doc page for a detailed description of general
+triclinic boxes.
+
+This option is provided as a convenience for users who may be
+converting data from solid-state crystallograhic representations or
+from DFT codes for input to LAMMPS.  However, as explained on the
+:doc:`Howto_triclinic <Howto_triclinic>` doc page, internally, LAMMPS
+only uses restricted triclinic simulation boxes.  This means the box
+and per-atom information (e.g. coordinates, velocities) LAMMPS stores
+are converted (rotated) from general to restricted triclinic form when
+the system is created.
+
+For dump output, if the :doc:`dump_modify triclinic/general
+<dump_modify>` command is used, the box description and per-atom
+coordinates and other per-atom vectors will be converted (rotated)
+from restricted to general form when each dump file snapshots is
+output.  This option can only be used if the simulation box was
+initially created as general triclinic.  If the option is not used,
+and the simulation box is general triclinic, then the dump file
+snapshots will reflect the internal restricted triclinic geometry.
+
+The dump_modify triclinic/general option affects 3 aspects of the dump
+file output.
+
+First, the format for the BOX BOUNDS is as follows
+
+.. parsed-literal::
+
+   ITEM: BOX BOUNDS abc origin
+   ax ay az originx
+   bx by bz originy
+   cx cy cz originz
+
+where the **A** edge vector of the box is (ax,ay,az) and similarly
+for **B** and **C**.  The origin of all 3 edge vectors is (originx,
+originy, originz).
+
+Second, the coordinates of each atom are converted (rotated) so that
+the atom is inside (or near) the general triclinic box defined by the
+**A**, **B**, **C** edge vectors.  For style *atom*, this only alters
+output for unscaled atom coords, via the :doc:`dump_modify scaled no
+<dump_modify>` setting. For style *custom*, this alters output for
+either unscaled or unwrapped output of atom coords, via the *x,y,z* or
+*xu,yu,zu* attributes.  For output of scaled atom coords by both
+styles, there is no difference bewteen restricted and general
+triclinic values.
+
+Third, the output for any attribute of the *custom* style which
+represents a per-atom vector quantity will be converted (rotated) to
+be oriented consistent with the general tricinic box and its
+orientation relative to the standard xyz coordinate axes.
+ 
+This applies to the following *custom* style attributes:
+ 
+* vx,vy,vz = atom velocities
+* fx,fy,fz = forces on atoms
+* mux,muy,muz = orientation of dipole moment of atom
+* omegax,omegay,omegaz = angular velocity of spherical particle
+* angmomx,angmomy,angmomz = angular momentum of aspherical particle
+* tqx,tqy,tqz = torque on finite-size particles
+
+For example, if the velocity of an atom in a restricted triclinic box
+is along the x-axis, then it will be output for a general triclinic
+box as a vector along the **A** edge vector of the box.
+
+.. note::
+
+   For style *custom*, the :doc:`dump_modify thresh <dump_modify>`
+   command may access per-atom attributes either directly or
+   indirectly through a compute or variable.  If the attribute is an
+   atom coordinate or one of the vectors mentioned above, its value
+   will *NOT* be a general triclinic (rotated) value.  Rather it will
+   be a restricted triclinic value.
+
+----------
+
 Arguments for different styles:
 
 The sections below describe per-atom, local, and per grid cell

doc/src/dump_modify.rst

@@ -17,7 +17,7 @@ Syntax
 * one or more keyword/value pairs may be appended
 
 * these keywords apply to various dump styles
-* keyword = *append* or *at* or *balance* or *buffer* or *colname* or *delay* or *element* or *every* or *every/time* or *fileper* or *first* or *flush* or *format* or *header* or *image* or *label* or *maxfiles* or *nfile* or *pad* or *pbc* or *precision* or *region* or *refresh* or *scale* or *sfactor* or *skip* or *sort* or *tfactor* or *thermo* or *thresh* or *time* or *units* or *unwrap*
+* keyword = *append* or *at* or *balance* or *buffer* or *colname* or *delay* or *element* or *every* or *every/time* or *fileper* or *first* or *flush* or *format* or *header* or *image* or *label* or *maxfiles* or *nfile* or *pad* or *pbc* or *precision* or *region* or *refresh* or *scale* or *sfactor* or *skip* or *sort* or *tfactor* or *thermo* or *thresh* or *time* or *triclinic/general* or *units* or *unwrap*
 
   .. parsed-literal::
 
@@ -74,12 +74,13 @@ Syntax
           -N = sort per-atom lines in descending order by the Nth column
        *tfactor* arg = time scaling factor (> 0.0)
        *thermo* arg = *yes* or *no*
-       *time* arg = *yes* or *no*
        *thresh* args = attribute operator value
          attribute = same attributes (x,fy,etotal,sxx,etc) used by dump custom style
          operator = "<" or "<=" or ">" or ">=" or "==" or "!=" or "\|\^"
          value = numeric value to compare to, or LAST
          these 3 args can be replaced by the word "none" to turn off thresholding
+       *time* arg = *yes* or *no*
+       *triclinic/general* arg = none
        *units* arg = *yes* or *no*
        *unwrap* arg = *yes* or *no*
 
@@ -802,8 +803,9 @@ region since the last dump.
    dump_modify ... thresh v_charge |^ LAST
 
 This will dump atoms whose charge has changed from an absolute value
-less than :math:`\frac12` to greater than :math:`\frac12` (or vice versa) since the last dump (e.g., due to reactions and subsequent charge equilibration in a
+less than :math:`\frac12` to greater than :math:`\frac12` (or vice
-reactive force field).
+versa) since the last dump (e.g., due to reactions and subsequent
+charge equilibration in a reactive force field).
 
 The choice of operators listed above are the usual comparison
 operators.  The XOR operation (exclusive or) is also included as "\|\^".
@@ -811,6 +813,18 @@ In this context, XOR means that if either the attribute or value is
 0.0 and the other is non-zero, then the result is "true" and the
 threshold criterion is met.  Otherwise it is not met.
 
+.. note::
+
+   For style *custom*, the *triclinic/general* keyword alters dump
+   output for general triclinic simulation boxes and their atoms.  See
+   the :doc:`dump <dump>` command for details of how this changes the
+   format of dump file snapstots.  The thresh keyword may access
+   per-atom attributes either directly or indirectly through a compute
+   or variable.  If the attribute is an atom coordinate or a per-atom
+   vector (such as velocity, force, or dipole moment), its value will
+   *NOT* be a general triclinic (rotated) value.  Rather it will be a
+   restricted triclinic value.
+
 ----------
 
 The *time* keyword only applies to the dump *atom*, *custom*, *local*,
@@ -835,6 +849,27 @@ The default setting is *no*\ .
 
 ----------
 
+The *triclinic/general* keyword only applies to the dump *atom* and
+*custom* styles.  It can only be used if the simulation box was
+created as a general triclinic box.  See the :doc:`Howto_triclinic
+<Howto_triclinic>` doc page for a detailed explanation of orthogonal,
+restricted triclinic, and general triclinic simulation boxes.
+
+If this keyword is used, the box information at the beginning of each
+snapshot will include information about the 3 arbitrary edge vectors
+**A**, **B**, **C** that define the general triclinic box as well as
+their origin.  The format is described on the :doc:`dump <dump>` doc
+page.
+
+The coordinates of each atom will be output as values in (or near) the
+general triclinic box.  Likewise, per-atom vector quantities such as
+velocity, omega, dipole moment, etc will have orientations consistent
+with the general triclinic box, meaning they will be rotated relative
+to the standard xyz coordinate axes.  See the :doc:`dump <dump>` doc
+page for a full list of which dump attributes this affects.
+
+----------
+
 The *units* keyword only applies to the dump *atom*, *custom*, and
 *local* styles (and their COMPRESS package versions *atom/gz*,
 *custom/gz* and *local/gz*\ ). If set to *yes*, each individual dump
@@ -922,6 +957,8 @@ The option defaults are
 * sort = off for dump styles *atom*, *custom*, *cfg*, and *local*
 * sort = id for dump styles *dcd*, *xtc*, and *xyz*
 * thresh = none
+* time = no
+* triclinic/general not specified
 * units = no
 * unwrap = no
 

src/create_atoms.cpp

@@ -113,6 +113,8 @@ void CreateAtoms::command(int narg, char **arg)
     xone[0] = utils::numeric(FLERR, arg[2], false, lmp);
     xone[1] = utils::numeric(FLERR, arg[3], false, lmp);
     xone[2] = utils::numeric(FLERR, arg[4], false, lmp);
+    if (domain->dimension == 2 && xone[2] != 0.0)
+      error->all(FLERR,"Create_atoms single for 2d simulation requires z coord = 0.0");
     iarg = 5;
   } else if (strcmp(arg[1], "random") == 0) {
     style = RANDOM;
@@ -359,7 +361,8 @@ void CreateAtoms::command(int narg, char **arg)
   
   if ((style == BOX) || (style == REGION)) {
     if (nbasis == 0) error->all(FLERR, "Cannot create atoms with undefined lattice");
-
+  }
+  
   // apply scaling factor for styles that use distance-dependent factors
 
   if (scaleflag) {
@@ -369,8 +372,8 @@ void CreateAtoms::command(int narg, char **arg)
     xone[2] *= domain->lattice->zlattice;
     } else if (style == RANDOM) {
       if (overlapflag) overlap *= domain->lattice->xlattice;
-    } else if (style == MESH) {     // NOTE to Axel - here is the rescaling of both params
+    } else if (style == MESH) {         // NOTE to Axel - here is the rescaling of both params
-      if (mesh_style == BISECT) {   //   by lattice spacings if units = lattice, similar to xone,overlap
+      if (mesh_style == BISECTION) {    //   by lattice spacings if units = lattice, similar to xone,overlap
         radthresh *= domain->lattice->xlattice;
       } else if (mesh_style = QUASIRANDOM) {
         mesh_density /= (domain->lattice->xlattice * domain->lattice->xlattice);

src/dump_atom.cpp

@@ -156,17 +156,21 @@ int DumpAtom::modify_param(int narg, char **arg)
     scale_flag = utils::logical(FLERR,arg[1],false,lmp);
     for (auto &item : keyword_user) item.clear();
     return 2;
-  } else if (strcmp(arg[0],"image") == 0) {
+  }
+  
+  if (strcmp(arg[0],"image") == 0) {
     if (narg < 2) error->all(FLERR,"Illegal dump_modify command");
     image_flag = utils::logical(FLERR,arg[1],false,lmp);
     for (auto &item : keyword_user) item.clear();
     return 2;
-  } else 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 (strcmp(arg[0],"triclinic/general") == 0) {
+    triclinic_general = 1;
     if (triclinic_general && !domain->triclinic_general)
-      error->all(FLERR,"Dump_modify triclinic/general invalid b/c simulation box is not");
+      error->all(FLERR,"Dump_modify triclinic/general cannot be used "
-    return 2;
+                 "if simulation box is not general triclinic");
+    return 1;
   }
   
   return 0;
@@ -410,9 +414,7 @@ void DumpAtom::header_item_triclinic_general(bigint ndump)
   }
   if (time_flag) fmt::print(fp,"ITEM: TIME\n{:.16}\n",compute_time());
 
-  fmt::print(fp,"ITEM: TIMESTEP\n{}\n"
+  fmt::print(fp,"ITEM: TIMESTEP\n{}\nITEM: NUMBER OF ATOMS\n{}\n", update->ntimestep, ndump);
-             "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"

src/dump_atom.h

@@ -35,7 +35,7 @@ class DumpAtom : public Dump {
  protected:
   int scale_flag;        // 1 if atom coords are scaled, 0 if no
   int image_flag;        // 1 if append box count to atom coords, 0 if no
-  int triclinic_general; // 1 if output box,coords for general triclinic
+  int triclinic_general; // 1 if output box & coords for general triclinic, 0 if no
 
   std::string columns;    // column labels
 

src/dump_custom.cpp

@@ -500,8 +500,8 @@ void DumpCustom::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);
@@ -573,9 +573,7 @@ void DumpCustom::header_item_triclinic_general(bigint ndump)
   }
   if (time_flag) fmt::print(fp,"ITEM: TIME\n{:.16}\n",compute_time());
 
-  fmt::print(fp,"ITEM: TIMESTEP\n{}\n"
+  fmt::print(fp,"ITEM: TIMESTEP\n{}\nITEM: NUMBER OF ATOMS\n{}\n", update->ntimestep, ndump);
-             "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"
@@ -1788,6 +1786,14 @@ int DumpCustom::modify_param(int narg, char **arg)
     return 2;
   }
 
+  if (strcmp(arg[0],"triclinic/general") == 0) {
+    triclinic_general = 1;
+    if (triclinic_general && !domain->triclinic_general)
+      error->all(FLERR,"Dump_modify triclinic/general cannot be used "
+                 "if simulation box is not general triclinic");
+    return 1;
+  }
+  
   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);
@@ -3350,11 +3356,11 @@ void DumpCustom::pack_tqz(int n)
 void DumpCustom::pack_tqx_triclinic_general(int n)
 {
   double **torque = atom->torque;
-  double torquetri[3];
+  double tqtri[3];
   
   for (int i = 0; i < nchoose; i++) {
-    domain->restricted_to_general_vector(torque[clist[i]],torquetri);
+    domain->restricted_to_general_vector(torque[clist[i]],tqtri);
-    buf[n] = torquetri[0];
+    buf[n] = tqtri[0];
     n += size_one;
   }
 }
@@ -3364,11 +3370,11 @@ void DumpCustom::pack_tqx_triclinic_general(int n)
 void DumpCustom::pack_tqy_triclinic_general(int n)
 {
   double **torque = atom->torque;
-  double torquetri[3];
+  double tqtri[3];
   
   for (int i = 0; i < nchoose; i++) {
-    domain->restricted_to_general_vector(torque[clist[i]],torquetri);
+    domain->restricted_to_general_vector(torque[clist[i]],tqtri);
-    buf[n] = torquetri[1];
+    buf[n] = tqtri[1];
     n += size_one;
   }
 }
@@ -3378,11 +3384,11 @@ void DumpCustom::pack_tqy_triclinic_general(int n)
 void DumpCustom::pack_tqz_triclinic_general(int n)
 {
   double **torque = atom->torque;
-  double torquetri[3];
+  double tqtri[3];
   
   for (int i = 0; i < nchoose; i++) {
-    domain->restricted_to_general_vector(torque[clist[i]],torquetri);
+    domain->restricted_to_general_vector(torque[clist[i]],tqtri);
-    buf[n] = torquetri[2];
+    buf[n] = tqtri[2];
     n += size_one;
   }
 }

src/dump_custom.h

@@ -36,8 +36,7 @@ class DumpCustom : public Dump {
  protected:
   int nevery;        // dump frequency for output
   char *idregion;    // region ID, nullptr if no region
-
+  int triclinic_general;  // 1 if output box & per-atom info for general triclinic
-  int triclinic_general;              // 1 if output box,x,v,f for general triclinic
   
   int nthresh;                        // # of defined thresholds
   int nthreshlast;                    // # of defined thresholds with value = LAST

src/lattice.cpp

@@ -17,6 +17,7 @@
 #include "comm.h"
 #include "domain.h"
 #include "error.h"
+#include "math_extra.h"
 #include "memory.h"
 #include "update.h"
 
@@ -24,6 +25,7 @@
 #include <cstring>
 
 using namespace LAMMPS_NS;
+using namespace MathExtra;
 
 #define BIG 1.0e30
 
@@ -237,7 +239,7 @@ Lattice::Lattice(LAMMPS *lmp, int narg, char **arg) : Pointers(lmp)
   if (nbasis == 0) error->all(FLERR,"No basis atoms in lattice");
   if (!orthogonal())
     error->all(FLERR,"Lattice orient vectors are not orthogonal");
-  if (!right_handed())
+  if (!right_handed_orientation())
     error->all(FLERR,"Lattice orient vectors are not right-handed");
   if (collinear())
     error->all(FLERR,"Lattice primitive vectors are collinear");
@@ -263,8 +265,7 @@ Lattice::Lattice(LAMMPS *lmp, int narg, char **arg) : Pointers(lmp)
       error->all(FLERR,"Lattice spacings are invalid");
   }
 
-  // requirements for a general triclinic lattice
+  // additional requirements for a general triclinic lattice
-  // right-handed requirement is checked by domain->general_to_restricted_rotation()
   // a123 prime are used to compute lattice spacings
   
   if (triclinic_general) {
@@ -278,6 +279,8 @@ Lattice::Lattice(LAMMPS *lmp, int narg, char **arg) : Pointers(lmp)
     if (orientz[0] != 0 || orientz[1] != 0 || orientz[2] != 1) oriented = 1;
     if (oriented)
       error->all(FLERR,"Lattice triclnic/general must have default orientation");
+    if (!right_handed_primitive())
+      error->all(FLERR,"Lattice triclnic/general a1,a2,a3 must be right-handed");
 
     double rotmat[3][3];
     domain->general_to_restricted_rotation(a1,a2,a3,rotmat,a1_prime,a2_prime,a3_prime);
@@ -289,8 +292,8 @@ Lattice::Lattice(LAMMPS *lmp, int narg, char **arg) : Pointers(lmp)
   
   if (strcmp(update->unit_style,"lj") == 0) {
     double vec[3];
-    cross(a2,a3,vec);
+    MathExtra::cross3(a2,a3,vec);
-    double volume = fabs(dot(a1,vec));
+    double volume = fabs(MathExtra::dot3(a1,vec));
     scale = pow(nbasis/volume/scale,1.0/dimension);
   }
 
@@ -389,7 +392,7 @@ int Lattice::orthogonal()
    x cross y must be in same direction as z
 ------------------------------------------------------------------------- */
 
-int Lattice::right_handed()
+int Lattice::right_handed_orientation()
 {
   int xy0 = orientx[1]*orienty[2] - orientx[2]*orienty[1];
   int xy1 = orientx[2]*orienty[0] - orientx[0]*orienty[2];
@@ -398,19 +401,33 @@ int Lattice::right_handed()
   return 1;
 }
 
+/* ----------------------------------------------------------------------
+   check righthandedness of a1,a2,a3 primitive vectors
+   x cross y must be in same direction as z
+------------------------------------------------------------------------- */
+
+int Lattice::right_handed_primitive()
+{
+  double vec[3];
+  MathExtra::cross3(a1,a2,vec);
+  if (MathExtra::dot3(vec,a3) <= 0.0) return 0;
+  return 1;
+}
+
 /* ----------------------------------------------------------------------
    check collinearity of each pair of primitive vectors
+   also checks if any primitive vector is zero-length
 ------------------------------------------------------------------------- */
 
 int Lattice::collinear()
 {
   double vec[3];
-  cross(a1,a2,vec);
+  MathExtra::cross3(a1,a2,vec);
-  if (dot(vec,vec) == 0.0) return 1;
+  if (MathExtra::len3(vec) == 0.0) return 1;
-  cross(a2,a3,vec);
+  MathExtra::cross3(a2,a3,vec);
-  if (dot(vec,vec) == 0.0) return 1;
+  if (MathExtra::len3(vec) == 0.0) return 1;
-  cross(a1,a3,vec);
+  MathExtra::cross3(a1,a3,vec);
-  if (dot(vec,vec) == 0.0) return 1;
+  if (MathExtra::len3(vec) == 0.0) return 1;
   return 0;
 }
 
@@ -589,26 +606,6 @@ void Lattice::add_basis(double x, double y, double z)
   nbasis++;
 }
 
-/* ----------------------------------------------------------------------
-   return x dot y
-------------------------------------------------------------------------- */
-
-double Lattice::dot(double *x, double *y)
-{
-  return x[0]*y[0] + x[1]*y[1] + x[2]*y[2];
-}
-
-/* ----------------------------------------------------------------------
-   z = x cross y
-------------------------------------------------------------------------- */
-
-void Lattice::cross(double *x, double *y, double *z)
-{
-  z[0] = x[1]*y[2] - x[2]*y[1];
-  z[1] = x[2]*y[0] - x[0]*y[2];
-  z[2] = x[0]*y[1] - x[1]*y[0];
-}
-
 /* ----------------------------------------------------------------------
    convert x,y,z from lattice coords to box coords (flag = 0)
    or from box coords to lattice coords (flag = 1)

src/lattice.h

@@ -56,7 +56,8 @@ class Lattice : protected Pointers {
   double a3_prime[3];
   
   int orthogonal();
-  int right_handed();
+  int right_handed_orientation();
+  int right_handed_primitive();
   int collinear();
   void setup_transform(double *, double *, double *);
   void add_basis(double, double, double);

tools/binary2txt.cpp

@@ -60,6 +60,7 @@ int main(int narg, char **arg)
   bigint ntimestep, natoms;
   int size_one, nchunk, triclinic;
   double xlo, xhi, ylo, yhi, zlo, zhi, xy, xz, yz;
+  double ax, ay, az, bx, by, bz, cx, cy, cz, abcx, abcy, abcz;
   int boundary[3][2];
   char boundstr[9];
 
@@ -133,17 +134,39 @@ int main(int narg, char **arg)
       fread(&natoms, sizeof(bigint), 1, fp);
       fread(&triclinic, sizeof(int), 1, fp);
       fread(&boundary[0][0], 6 * sizeof(int), 1, fp);
-      fread(&xlo, sizeof(double), 1, fp);
+
-      fread(&xhi, sizeof(double), 1, fp);
+      if (triclinic == 0) {
-      fread(&ylo, sizeof(double), 1, fp);
+        fread(&xlo, sizeof(double), 1, fp);
-      fread(&yhi, sizeof(double), 1, fp);
+        fread(&xhi, sizeof(double), 1, fp);
-      fread(&zlo, sizeof(double), 1, fp);
+        fread(&ylo, sizeof(double), 1, fp);
-      fread(&zhi, sizeof(double), 1, fp);
+        fread(&yhi, sizeof(double), 1, fp);
-      if (triclinic) {
+        fread(&zlo, sizeof(double), 1, fp);
+        fread(&zhi, sizeof(double), 1, fp);
+      } else if (triclinic == 1) {
+        fread(&xlo, sizeof(double), 1, fp);
+        fread(&xhi, sizeof(double), 1, fp);
+        fread(&ylo, sizeof(double), 1, fp);
+        fread(&yhi, sizeof(double), 1, fp);
+        fread(&zlo, sizeof(double), 1, fp);
+        fread(&zhi, sizeof(double), 1, fp);
         fread(&xy, sizeof(double), 1, fp);
         fread(&xz, sizeof(double), 1, fp);
         fread(&yz, sizeof(double), 1, fp);
+      } else if (triclinic == 2) {
+        fread(&ax, sizeof(double), 1, fp);
+        fread(&ay, sizeof(double), 1, fp);
+        fread(&az, sizeof(double), 1, fp);
+        fread(&bx, sizeof(double), 1, fp);
+        fread(&by, sizeof(double), 1, fp);
+        fread(&bz, sizeof(double), 1, fp);
+        fread(&cx, sizeof(double), 1, fp);
+        fread(&cy, sizeof(double), 1, fp);
+        fread(&cz, sizeof(double), 1, fp);
+        fread(&abcx, sizeof(double), 1, fp);
+        fread(&abcy, sizeof(double), 1, fp);
+        fread(&abcz, sizeof(double), 1, fp);
       }
+      
       fread(&size_one, sizeof(int), 1, fp);
 
       if (magic_string && revision > 0x0001) {
@@ -201,16 +224,21 @@ int main(int narg, char **arg)
       }
       boundstr[8] = '\0';
 
-      if (!triclinic) {
+      if (triclinic == 0) {
         fprintf(fptxt, "ITEM: BOX BOUNDS %s\n", boundstr);
         fprintf(fptxt, "%-1.16e %-1.16e\n", xlo, xhi);
         fprintf(fptxt, "%-1.16e %-1.16e\n", ylo, yhi);
         fprintf(fptxt, "%-1.16e %-1.16e\n", zlo, zhi);
-      } else {
+      } else if (triclinic == 1) {
         fprintf(fptxt, "ITEM: BOX BOUNDS xy xz yz %s\n", boundstr);
         fprintf(fptxt, "%-1.16e %-1.16e %-1.16e\n", xlo, xhi, xy);
         fprintf(fptxt, "%-1.16e %-1.16e %-1.16e\n", ylo, yhi, xz);
         fprintf(fptxt, "%-1.16e %-1.16e %-1.16e\n", zlo, zhi, yz);
+      } else if (triclinic == 2) {
+        fprintf(fptxt, "ITEM: BOX BOUNDS abc origin %s\n", boundstr);
+        fprintf(fptxt, "%-1.16e %-1.16e %-1.16e %-1.16e\n", ax, ay, az, abcx);
+        fprintf(fptxt, "%-1.16e %-1.16e %-1.16e %-1.16e\n", bx, by, bz, abcy);
+        fprintf(fptxt, "%-1.16e %-1.16e %-1.16e %-1.16e\n", cx, cy, cz, abcz);
       }
 
       if (columns)