From dfafdff2093fa347c44ff23fea381b620e670df8 Mon Sep 17 00:00:00 2001 From: Steve Plimpton Date: Mon, 20 Nov 2023 12:06:58 -0700 Subject: [PATCH] finalized dump command support for general triclinic --- doc/src/create_atoms.rst | 141 ++++++++++++++++++++------------------- doc/src/dump.rst | 108 +++++++++++++++++++++++++++--- doc/src/dump_modify.rst | 45 +++++++++++-- src/create_atoms.cpp | 9 ++- src/dump_atom.cpp | 20 +++--- src/dump_atom.h | 2 +- src/dump_custom.cpp | 34 ++++++---- src/dump_custom.h | 3 +- src/lattice.cpp | 61 ++++++++--------- src/lattice.h | 3 +- tools/binary2txt.cpp | 46 ++++++++++--- 11 files changed, 319 insertions(+), 153 deletions(-) diff --git a/doc/src/create_atoms.rst b/doc/src/create_atoms.rst index 2668146852..71303a7d4f 100644 --- a/doc/src/create_atoms.rst +++ b/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 ` command. However, you can first + use the :doc:`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 ` 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 ` 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 ` +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 ` 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 `. Or run initial dynamics with + :doc:`pair_style soft ` or with :doc:`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 `. -.. 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 ` 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 ` -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 ` 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 `. Or run initial dynamics with - :doc:`pair_style soft ` or with :doc:`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 ` command. - However, you can first use the :doc:`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 - ` 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* diff --git a/doc/src/dump.rst b/doc/src/dump.rst index e5885dc25d..93bffa3fdc 100644 --- a/doc/src/dump.rst +++ b/doc/src/dump.rst @@ -278,16 +278,20 @@ format ` 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 -viewing a dump file. Many post-processing tools either included with -LAMMPS or third-party tools can read this format, as does the +simple LAMMPS-specific text format that is mostly self-explanatory +when viewing a dump file. Many post-processing tools either included +with LAMMPS or third-party tools can read this format, as does the :doc:`rerun ` command. See tools described on the :doc:`Tools ` doc page for examples, including `Pizza.py `_. For all these styles, the dimensions of the simulation box are -included in each snapshot. For an orthogonal simulation box this -information is formatted as: +included in each snapshot. The simulation box in LAMMPS can be +defined in one of 3 ways: orthogonal, restricted triclinic, and +general triclinic. See the :doc:`Howto triclinic ` +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 ` command for details. -For triclinic simulation boxes (non-orthogonal), an orthogonal -bounding box which encloses the triclinic simulation box is output, -along with the three tilt factors (*xy*, *xz*, *yz*) of the triclinic box, -formatted as follows: +For a restricted triclinic simulation box, an orthogonal bounding box +which encloses the restricted triclinic simulation box is output, +along with the three tilt factors (*xy*, *xz*, *yz*) of the triclinic +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 +` 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 ` 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 +` 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 +` 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 ` + 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 diff --git a/doc/src/dump_modify.rst b/doc/src/dump_modify.rst index 2d84f28836..31f1cd214e 100644 --- a/doc/src/dump_modify.rst +++ b/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 -reactive force field). +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 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 ` 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 +` 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 ` 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 ` 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 diff --git a/src/create_atoms.cpp b/src/create_atoms.cpp index d510014191..59b0ad4ab6 100644 --- a/src/create_atoms.cpp +++ b/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 - if (mesh_style == BISECT) { // by lattice spacings if units = lattice, similar to xone,overlap + } else if (style == MESH) { // NOTE to Axel - here is the rescaling of both params + 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); diff --git a/src/dump_atom.cpp b/src/dump_atom.cpp index 38c8431269..a029a68135 100644 --- a/src/dump_atom.cpp +++ b/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"); - return 2; + error->all(FLERR,"Dump_modify triclinic/general cannot be used " + "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" - "ITEM: NUMBER OF ATOMS\n{}\n", - update->ntimestep, ndump); + fmt::print(fp,"ITEM: TIMESTEP\n{}\nITEM: 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" diff --git a/src/dump_atom.h b/src/dump_atom.h index df21c7788c..0b0d3a4f05 100644 --- a/src/dump_atom.h +++ b/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 diff --git a/src/dump_custom.cpp b/src/dump_custom.cpp index 9554247aff..6708a7c410 100644 --- a/src/dump_custom.cpp +++ b/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; - fwrite(&general_tri,sizeof(int),1,fp); + int triclinic_general_flag = 2; + 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" - "ITEM: NUMBER OF ATOMS\n{}\n", - update->ntimestep, ndump); + fmt::print(fp,"ITEM: TIMESTEP\n{}\nITEM: 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); - buf[n] = torquetri[0]; + domain->restricted_to_general_vector(torque[clist[i]],tqtri); + 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); - buf[n] = torquetri[1]; + domain->restricted_to_general_vector(torque[clist[i]],tqtri); + 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); - buf[n] = torquetri[2]; + domain->restricted_to_general_vector(torque[clist[i]],tqtri); + buf[n] = tqtri[2]; n += size_one; } } diff --git a/src/dump_custom.h b/src/dump_custom.h index 60070ddf62..b22d03f9b5 100644 --- a/src/dump_custom.h +++ b/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,x,v,f for general triclinic + int triclinic_general; // 1 if output box & per-atom info for general triclinic int nthresh; // # of defined thresholds int nthreshlast; // # of defined thresholds with value = LAST diff --git a/src/lattice.cpp b/src/lattice.cpp index d29a378589..6992828e3a 100644 --- a/src/lattice.cpp +++ b/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 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 - // right-handed requirement is checked by domain->general_to_restricted_rotation() + // additional requirements for a general triclinic lattice // 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); - double volume = fabs(dot(a1,vec)); + MathExtra::cross3(a2,a3,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); - if (dot(vec,vec) == 0.0) return 1; - cross(a2,a3,vec); - if (dot(vec,vec) == 0.0) return 1; - cross(a1,a3,vec); - if (dot(vec,vec) == 0.0) return 1; + MathExtra::cross3(a1,a2,vec); + if (MathExtra::len3(vec) == 0.0) return 1; + MathExtra::cross3(a2,a3,vec); + if (MathExtra::len3(vec) == 0.0) return 1; + MathExtra::cross3(a1,a3,vec); + 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) diff --git a/src/lattice.h b/src/lattice.h index 74f213bdb4..8f545de2ab 100644 --- a/src/lattice.h +++ b/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); diff --git a/tools/binary2txt.cpp b/tools/binary2txt.cpp index b77614a8fc..2369057324 100644 --- a/tools/binary2txt.cpp +++ b/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); - fread(&ylo, sizeof(double), 1, fp); - fread(&yhi, sizeof(double), 1, fp); - fread(&zlo, sizeof(double), 1, fp); - fread(&zhi, sizeof(double), 1, fp); - if (triclinic) { + + if (triclinic == 0) { + 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); + } 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)