From 56b2c7ed4651c3bbb3ba7457d7aa9f5cc0ba4c03 Mon Sep 17 00:00:00 2001 From: Steve Plimpton Date: Fri, 17 Nov 2023 17:10:45 -0700 Subject: [PATCH] alter how lattice interacts with create_box and create_atoms for general triclinic --- doc/src/Howto_triclinic.rst | 6 +- doc/src/create_atoms.rst | 112 ++++++++++++++++++++----- doc/src/create_box.rst | 56 ++++++------- doc/src/lattice.rst | 68 ++++++++++++++- src/atom_vec_line.cpp | 2 +- src/create_atoms.cpp | 54 ++++++------ src/create_box.cpp | 12 ++- src/domain.cpp | 160 ++++++++++++++++++++---------------- src/domain.h | 3 + src/lattice.cpp | 91 ++++++++++++-------- src/lattice.h | 25 +++--- 11 files changed, 383 insertions(+), 206 deletions(-) diff --git a/doc/src/Howto_triclinic.rst b/doc/src/Howto_triclinic.rst index 6cbc0644cd..525c3e0f1b 100644 --- a/doc/src/Howto_triclinic.rst +++ b/doc/src/Howto_triclinic.rst @@ -56,9 +56,9 @@ at (xlo,ylo,zhi) and 3 edge vectors **A** = (ax,ay,az), **B** = (bx,by,bz), **C** = (cx,cy,cz) which can be arbitrary vectors, so long as they are non-zero, distinct, and not co-planar. In addition, they must define a right-handed system, such that (**A** cross **B**) -points in the direction of **C**. A left-handed system can be -converted to a right-handed system by simply swapping the order of any -pair of the **A**, **B**, **C** vectors. +points in the direction of **C**. Note that a left-handed system can +be converted to a right-handed system by simply swapping the order of +any pair of the **A**, **B**, **C** vectors. The 4 commands listed above for defining orthogonal simulation boxes have triclinic options which allow for specification of the origin and diff --git a/doc/src/create_atoms.rst b/doc/src/create_atoms.rst index 5d1e7c872c..c9e8f3e840 100644 --- a/doc/src/create_atoms.rst +++ b/doc/src/create_atoms.rst @@ -86,25 +86,25 @@ Description """"""""""" This command creates atoms (or molecules) within the simulation box, -either on a lattice, or a single atom (or molecule), or on a surface -defined by a triangulated mesh, or a random collection of atoms (or -molecules). It is an alternative to reading in atom coordinates -explicitly via a :doc:`read_data ` or :doc:`read_restart -` command. A simulation box must already exist, which is -typically created via the :doc:`create_box ` command. -Before using this command, a lattice must also be defined using the -:doc:`lattice ` command, unless you specify the *single* style -with units = box or the *random* style. For the remainder of this doc -page, a created atom or molecule is referred to as a "particle". +either on a lattice, or at a single specified location, or randomly, +or on a surface defined by a triangulated mesh. It is an alternative +to reading in atom coordinates explicitly via a :doc:`read_data +` or :doc:`read_restart ` command. A +simulation box must already exist, which is typically created via the +:doc:`create_box ` command. Before using this command, a +lattice must typically also be defined using the :doc:`lattice +` command, unless you specify the *single* style with units = +box or the *random* style. For the remainder of this doc page, a +created atom or molecule is referred to as a "particle". If created particles are individual atoms, they are assigned the specified atom *type*, though this can be altered via the *basis* keyword as discussed below. If molecules are being created, the type -of each atom in the created molecule is specified in the file read by -the :doc:`molecule ` command, and those values are 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). +of each atom in the created molecule is specified in a specified file +read by the :doc:`molecule ` command, and those values are +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). For the *box* style, the create_atoms command fills the entire simulation box with particles on the lattice. If your simulation box @@ -126,6 +126,69 @@ periodic boundaries. If this is desired, you should either use the *box* style, or tweak the region size to get precisely the particles you want. + +------------------ + +WORK on this + + +If a general triclinic simulation box is defined ... + + + +As noted above, general triclinic boxes in LAMMPS allow for arbitrary +edge vectors **A**, **B**, **C**. The only restrictions are that the +three vectors be distinct, non-zero, and not co-planar. They must +also define a right-handed system such that (**A** x **B**) points in +the direction of **C**. Note that a left-handed system can be +converted to a right-handed system by simply swapping the order of any +pair of the **A**, **B**, **C** vectors. + +To create a general triclinic boxes, the region is specified as NULL +and the next 6 parameters (alo,ahi,blo,bhi,clo,chi) define the three +edge vectors **A**, **B**, **C** using additional information +previously defind by the :doc:`lattice ` command. + +The lattice must be of style *custom* and use its *triclinic/general* +option. This insures the lattice satisfies the restrictions listed +above. The *a1, *a2*, *a3* settings of the :doc:`lattice ` +command define the edge vectors of a unit cell of the general +triclinic lattice. This command uses them to define the three edge +vectors and origin of the general triclinic box as: + + +explain region is applied after conversion to restricted triclinic atom coords + +explain general tri for box and region styles +must use lattice triclinic/general +paragraph about DFT motivation +doc that single, random, mesh operate on restricted triclinic box + +------------------ + + +.. note:: + + LAMMPS allows specification of general triclinic simulation boxes + 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 defined by this command and per-atom information + (e.g. coordinates, velocities) defined by the :doc:`create_atoms + ` command are converted from general to restricted + triclinic form when the two commands are invoked. The + ` doc page also discusses other LAMMPS commands + which can input/output general triclinic representations of the + simulation box and per-atom data. + + + + + + + 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 @@ -463,12 +526,19 @@ on a single CPU core. ----- The *units* keyword determines the meaning of the distance units used -to specify the coordinates of the one particle created by the *single* -style, or the overlap distance *Doverlap* by the *overlap* keyword. A -*box* value selects standard distance units as defined by the -:doc:`units ` command (e.g., :math:`\AA` for -units = *real* or *metal*\ . A *lattice* value means the distance units are in -lattice spacings. +by parameters for various styles. A *box* value selects standard +distance units as defined by the :doc:`units ` command (e.g., +:math:`\AA` for units = *real* or *metal*\ . A *lattice* value means +the distance units are in lattice spacings. These are affected settings: + +* for *single* style: coordinates of the particle created +* for *random* style: overlap distance *Doverlap* by the *overlap* keyword +* for *mesh* style: *bisect* threshold valeu for *meshmode* = *bisect* +* for *mesh* style: *radthresh* value for *meshmode* = *bisect* +* for *mesh* style: *density* value for *meshmode* = *qrand* + +Since *density* represents an area (distance ^2), the lattice spacing +factor is also squared. ---------- diff --git a/doc/src/create_box.rst b/doc/src/create_box.rst index b5ac4ac907..11b2878fe4 100644 --- a/doc/src/create_box.rst +++ b/doc/src/create_box.rst @@ -91,14 +91,14 @@ parallelepiped has an "origin" at (xlo,ylo,zlo) and three edge vectors starting from the origin given by :math:`\vec a = (x_\text{hi}-x_\text{lo},0,0)`; :math:`\vec b = (xy,y_\text{hi}-y_\text{lo},0)`; and :math:`\vec c = -(xz,yz,z_\text{hi}-z_\text{lo})`. This is a restricted triclinic box -because the three edge vectors cannot be defined in arbitrary -(general) directions. The parameters *xy*\ , *xz*\ , and *yz* can be -0.0 or positive or negative values and are called "tilt factors" -because they are the amount of displacement applied to faces of an -originally orthogonal box to transform it into the parallelepiped. -For a 2d simulation, the zlo and zhi values of the simulation box must -straddle zero. +(xz,yz,z_\text{hi}-z_\text{lo})`. In LAMMPS lingo, this is a +restricted triclinic box because the three edge vectors cannot be +defined in arbitrary (general) directions. The parameters *xy*\ , +*xz*\ , and *yz* can be 0.0 or positive or negative values and are +called "tilt factors" because they are the amount of displacement +applied to faces of an originally orthogonal box to transform it into +the parallelepiped. For a 2d simulation, the zlo and zhi values of +the simulation box must straddle zero. Typically a *prism* region used with the create_box command should have tilt factors :math:`(xy,xz,yz)` that do not skew the box more @@ -164,37 +164,35 @@ using the :doc:`change box ` command with its *ortho* and ---------- -As noted above, general triclinic boxes allow for arbitrary edge -vectors **A**, **B**, **C*. The only restrictions are that the three -vectors be distinct, non-zero, and not co-planar. They must also -define a right-handed system such that (**A** x **B**) points in the -direction of **C**. Note that a left-handed system can be converted -to a right-handed system by simply swapping the order of any pair of -the **A**, **B**, **C** vectors. +As noted above, general triclinic boxes in LAMMPS allow for arbitrary +edge vectors **A**, **B**, **C**. The only restrictions are that the +three vectors be distinct, non-zero, and not co-planar. They must +also define a right-handed system such that (**A** x **B**) points in +the direction of **C**. Note that a left-handed system can be +converted to a right-handed system by simply swapping the order of any +pair of the **A**, **B**, **C** vectors. To create a general triclinic boxes, the region is specified as NULL and the next 6 parameters (alo,ahi,blo,bhi,clo,chi) define the three edge vectors **A**, **B**, **C** using additional information -previousl set by the :doc:`lattice ` command. +previously defind by the :doc:`lattice ` command. -The lattice must be of style *custom* with a default orientation (no -use of the *orient* keyword to change the default values). The *a1, -*a2*, *a3* settings of the :doc:`lattice ` command define the -edge vectors of a unit cell of the lattice. The three edge vectors of -the general triclinic box are defined as +The lattice must be of style *custom* and use its *triclinic/general* +option. This insures the lattice satisfies the restrictions listed +above. The *a1, *a2*, *a3* settings of the :doc:`lattice ` +command define the edge vectors of a unit cell of the general +triclinic lattice. This command uses them to define the three edge +vectors and origin of the general triclinic box as: * **A** = (ahi-alo) * *a1* * **B** = (bhi-blo) * *a2* * **C** = (chi-clo) * *a3* +* origin = (alo*a1 + blo*a2 + clo*a3) -The origin of the general triclinic box is at (alo*a1 + blo*a2 + -clo*a3) with an additional offset due to the *origin* setting of the -lattice command (zero vector by default). - -For 2d simulations, **C** = (0,0,1) is required, and the z-component -of the simulation box origin must be -0.5. The easy way to do this is -to specify clo = -0.5 and chi = 0.5 with the :doc:`lattice ` -command default of a3 = (0,0,1). +For 2d general triclinic boxes, **C** = (0,0,1) is required, and the +z-component of the simulation box origin must be -0.5. The easy way +to do this is to specify clo = -0.5 and chi = 0.5 and use the +:doc:`lattice ` command default for a3 = (0,0,1). .. note:: diff --git a/doc/src/lattice.rst b/doc/src/lattice.rst index 3306c58e8c..235b8b30da 100644 --- a/doc/src/lattice.rst +++ b/doc/src/lattice.rst @@ -19,7 +19,7 @@ Syntax scale = lattice constant in distance units (for all other units) * zero or more keyword/value pairs may be appended -* keyword = *origin* or *orient* or *spacing* or *a1* or *a2* or *a3* or *basis* +* keyword = *origin* or *orient* or *spacing* or *a1* or *a2* or *a3* or *basis* or *triclinic/general* .. parsed-literal:: @@ -34,6 +34,7 @@ Syntax x,y,z = primitive vector components that define unit cell *basis* values = x y z x,y,z = fractional coords of a basis atom (0 <= x,y,z < 1) + *triclinic/general* values = no values, assume lattice tiles Examples """""""" @@ -44,7 +45,7 @@ Examples lattice hex 0.85 lattice sq 0.8 origin 0.0 0.5 0.0 orient x 1 1 0 orient y -1 1 0 lattice custom 3.52 a1 1.0 0.0 0.0 a2 0.5 1.0 0.0 a3 0.0 0.0 0.5 & - basis 0.0 0.0 0.0 basis 0.5 0.5 0.5 + basis 0.0 0.0 0.0 basis 0.5 0.5 0.5 triclinic/general lattice none 2.0 Description @@ -196,6 +197,57 @@ the Z direction. ---------- +The *triclinic/general* option specifies that the defined lattice is +for use with a general triclinic simulation box, as opposed to an +orthogonal or restricted triclinic box. The :doc:`Howto triclinic +` doc page explains all 3 kinds of simluation boxes +LAMMPS supports. + +If this option is specified, a *custom* lattice style must be used. +The *a1*, *a2*, *a3* vectors should define the edge vectors of a +single unit cell of the lattice with one or more basis atoms. They +edge vectors can be arbitrary so long as they are non-zero, distinct, +and not co-planar. In addition, they must define a right-handed +system, such that (*a1* cross *a2*) points in the direction of *a3*. +Note that a left-handed system can be converted to a right-handed +system by simply swapping the order of any pair of the *a1*, *a2*, +*a3* vectors. + +If this option is used, the *origin* and *orient* settings must have +their default values. + +The :doc:`create_box ` command can be used to create a +general triclinic box that replicates the *a1*, *a2*, *a3* unit cell +vectors in each direction to create the 3 arbitrary edge vectors of +the overall simulation box. It requires a lattice with the +*triclinic/general* option. + +Likewise, the :doc:`create_atoms ` command can be used +to add atoms (or molecules) to a general triclinic box which lie on +the lattice points defined by *a1*, *a2*, *a3* and the unit cell basis +atoms. To do this, it also requires a lattice with the +*triclinic/general* option. + +.. note:: + + LAMMPS allows specification of general triclinic lattices and + simulation boxes 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) + defined by the :doc:`create_box ` and + :doc:`create_atoms ` commands are converted from + general to restricted triclinic form when the two commands are + invoked. It also means that any other commands which use lattice + spacings from this command (e.g. the region command), will be + operating on a restricted triclinic simulation box, even if the + *triclinic/general* option was used to define the lattice. See the + next section for details. + +---------- + Several LAMMPS commands have the option to use distance units that are inferred from "lattice spacings" in the x,y,z box directions. E.g. the :doc:`region ` command can create a block of size @@ -219,6 +271,18 @@ coordinates of the 8 corner points of the modified unit cell (4 in 2d). Similarly, the Y and Z lattice spacings are defined as the difference in the min/max of the y and z coordinates. +.. note:: + + If the *triclinic/general* option is specified, the unit cell + defined by *a1*, *a2*, *a3* edge vectors is first converted to a + restricted triclinic orientation, which is a rotation operation. + The min/max extent of the 8 corner points is then determined, as + described in the preceeding paragraph, to set the lattice + spacings. As explained for the *triclinic/general* option above, + this is because any use of the lattice spacings by other commands + will be for a restricted triclinic simulation box, not a general + triclinic box. + Note that if the unit cell is orthogonal with axis-aligned edges (no rotation via the *orient* keyword), then the lattice spacings in each dimension are simply the scale factor (described above) multiplied by diff --git a/src/atom_vec_line.cpp b/src/atom_vec_line.cpp index f5544c08ee..21a9cc880a 100644 --- a/src/atom_vec_line.cpp +++ b/src/atom_vec_line.cpp @@ -374,7 +374,7 @@ void AtomVecLine::data_atom_bonus(int m, const std::vector &values) 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]); + domain->remap_near(coords,x[m]); x2 = coords[0]; y2 = coords[1]; // calculate length and theta diff --git a/src/create_atoms.cpp b/src/create_atoms.cpp index ef19085ce5..a24fba4612 100644 --- a/src/create_atoms.cpp +++ b/src/create_atoms.cpp @@ -154,9 +154,10 @@ void CreateAtoms::command(int narg, char **arg) overlapflag = 0; maxtry = DEFAULT_MAXTRY; radscale = 1.0; - radthresh = domain->lattice->xlattice; mesh_style = BISECTION; - mesh_density = 1.0; + radthresh = domain->lattice->xlattice; // NOTE to Axel - I think this should be 1.0 by default + mesh_density = 1.0; // similar to how this setting is 1.0 + // see rescaling of both below if units = lattice nbasis = domain->lattice->nbasis; basistype = new int[nbasis]; for (int i = 0; i < nbasis; i++) basistype[i] = ntype; @@ -354,19 +355,27 @@ void CreateAtoms::command(int narg, char **arg) } } - // demand non-none lattice be defined for BOX and REGION - // else setup scaling for SINGLE and RANDOM - // could use domain->lattice->lattice2box() to do conversion of - // lattice to box, but not consistent with other uses of units=lattice - // triclinic remapping occurs in add_single() - + // require non-none lattice be defined for BOX or REGION styles + if ((style == BOX) || (style == REGION)) { if (nbasis == 0) error->all(FLERR, "Cannot create atoms with undefined lattice"); - } else if (scaleflag == 1) { + + // apply scaling factor for styles that use distance-dependent factors + + if (scaleflag) { + if (style == SINGLE) { xone[0] *= domain->lattice->xlattice; xone[1] *= domain->lattice->ylattice; xone[2] *= domain->lattice->zlattice; - overlap *= domain->lattice->xlattice; + } 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 + radthresh *= domain->lattice->xlattice; + } else if (mesh_style = QUASIRANDOM) { + mesh_density /= (domain->lattice->xlattice * domain->lattice->xlattice); + } + } } // set bounds for my proc in sublo[3] & subhi[3] @@ -710,7 +719,7 @@ void CreateAtoms::add_single() void CreateAtoms::add_random() { - double xlo, ylo, zlo, xhi, yhi, zhi, zmid; + double xlo, ylo, zlo, xhi, yhi, zhi; double delx, dely, delz, distsq, odistsq; double lamda[3], *coord; double *boxlo, *boxhi; @@ -729,7 +738,6 @@ void CreateAtoms::add_random() for (int ii = 0; ii < 30; ii++) random->uniform(); // bounding box for atom creation - // in real units, even if triclinic // only limit bbox by region if its bboxflag is set (interior region) if (triclinic == 0) { @@ -739,7 +747,6 @@ void CreateAtoms::add_random() yhi = domain->boxhi[1]; zlo = domain->boxlo[2]; zhi = domain->boxhi[2]; - zmid = zlo + 0.5 * (zhi - zlo); } else { xlo = domain->boxlo_bound[0]; xhi = domain->boxhi_bound[0]; @@ -747,7 +754,6 @@ void CreateAtoms::add_random() yhi = domain->boxhi_bound[1]; zlo = domain->boxlo_bound[2]; zhi = domain->boxhi_bound[2]; - zmid = zlo + 0.5 * (zhi - zlo); boxlo = domain->boxlo_lamda; boxhi = domain->boxhi_lamda; } @@ -783,7 +789,7 @@ void CreateAtoms::add_random() xone[0] = xlo + random->uniform() * (xhi - xlo); xone[1] = ylo + random->uniform() * (yhi - ylo); xone[2] = zlo + random->uniform() * (zhi - zlo); - if (domain->dimension == 2) xone[2] = zmid; + if (domain->dimension == 2) xone[2] = 0.0; if (region && (region->match(xone[0], xone[1], xone[2]) == 0)) continue; if (varflag && vartest(xone) == 0) continue; @@ -1168,16 +1174,12 @@ void CreateAtoms::add_mesh(const char *filename) void CreateAtoms::add_lattice() { // add atoms on general triclinic lattice if Domain has setting for it - // verify lattice is valid for general triclinic - - int triclinic_general = domain->triclinic_general; + // verify lattice was defined with triclinic/general option - if (triclinic_general) { - if (!domain->lattice->is_custom()) - error->all(FLERR,"Create_atoms for general triclinic requires custom lattice"); - if (domain->lattice->is_oriented()) - error->all(FLERR,"Create_atoms for general triclinic requires lattice with default orientation"); - } + if (!domain->triclinic_general && domain->lattice->is_general_triclinic()) + error->all(FLERR,"Create_atoms for non general triclinic box cannot use triclinic/general lattice"); + if (domain->triclinic_general && !domain->lattice->is_general_triclinic()) + error->all(FLERR,"Create_atoms for general triclinic box requires triclnic/general lattice"); // convert 8 corners of my subdomain from box coords to lattice coords // for orthogonal, use corner pts of my subbox @@ -1224,7 +1226,7 @@ void CreateAtoms::add_lattice() // compute new bounding box (xyz min/max) in lattice coords // for orthogonal or restricted triclinic, use 8 corner points of bbox lo/hi - if (!triclinic_general) { + if (!domain->triclinic_general) { domain->lattice->bbox(1, bboxlo[0], bboxlo[1], bboxlo[2], xmin, ymin, zmin, xmax, ymax, zmax); domain->lattice->bbox(1, bboxhi[0], bboxlo[1], bboxlo[2], xmin, ymin, zmin, xmax, ymax, zmax); domain->lattice->bbox(1, bboxlo[0], bboxhi[1], bboxlo[2], xmin, ymin, zmin, xmax, ymax, zmax); @@ -1238,7 +1240,7 @@ void CreateAtoms::add_lattice() // new set of 8 points is no longer an orthogonal bounding box // instead invoke lattice->bbox() on each of 8 points - } else if (triclinic_general) { + } else if (domain->triclinic_general) { double point[3]; point[0] = bboxlo[0]; point[1] = bboxlo[1]; point[2] = bboxlo[2]; diff --git a/src/create_box.cpp b/src/create_box.cpp index 4df23d286f..e4ebe4c0cb 100644 --- a/src/create_box.cpp +++ b/src/create_box.cpp @@ -106,10 +106,8 @@ void CreateBox::command(int narg, char **arg) // ABC edge vectors + origin to restricted triclinic } else if (triclinic_general) { - if (!domain->lattice->is_custom()) - error->all(FLERR,"Create_box with no region requires custom lattice"); - if (domain->lattice->is_oriented()) - error->all(FLERR,"Create_box with no region requires lattice with default orientation"); + if (!domain->lattice->is_general_triclinic()) + error->all(FLERR,"Create_box for general triclinic requires triclnic/general lattice"); if (iarg + 6 > narg) utils::missing_cmd_args(FLERR, "create_box general triclinic", error); @@ -154,11 +152,11 @@ void CreateBox::command(int narg, char **arg) if (domain->dimension == 2) { if (cvec[0] != 0.0 || cvec[1] != 0.0 || cvec[2] != 1.0 || origin[2] != -0.5) - error->all(FLERR,"Create_box C vector and/or origin is invalid " + error->all(FLERR,"Create_box C edge vector and/or origin is invalid " "for 2d simulation with general triclinic box"); } - // define general triclinic box within Domain + // define general triclinic box within Domain class domain->define_general_triclinic(avec,bvec,cvec,origin); } @@ -247,7 +245,7 @@ void CreateBox::command(int narg, char **arg) error->all(FLERR, "Unknown create_box keyword: {}", arg[iarg]); } - // problem setup using info from header + // setup the simulation box and initial system // deallocate/grow ensures any extra settings are used for topology arrays // necessary in case no create_atoms is performed diff --git a/src/domain.cpp b/src/domain.cpp index e9a18b897e..3e49255b57 100644 --- a/src/domain.cpp +++ b/src/domain.cpp @@ -571,88 +571,24 @@ void Domain::define_general_triclinic(double *avec_caller, double *bvec_caller, cvec[1] = cvec_caller[1]; cvec[2] = cvec_caller[2]; - boxlo[0] = origin_caller[0]; - boxlo[1] = origin_caller[1]; - boxlo[2] = origin_caller[2]; - // error check on cvec for 2d systems if (dimension == 2 && (cvec[0] != 0.0 || cvec[1] != 0.0)) error->all(FLERR,"General triclinic box edge vector C invalid for 2d system"); - // error checks for 3d systems - // A,B,C cannot be co-planar - // A x B must point in C direction (right-handed) - - 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"); - if (dot < 0.0) - error->all(FLERR,"General triclinic box edge vectors must be right-handed"); - - // 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 - // +xaxis = unit vector to rotate B1 around - // theta2 = angle of rotation calculated from - // Byz1 dot yunit = B1y = |Byz1| cos(theta2) - // theta2 via acos() is positive (0 to PI) - // positive is valid if B1z < 0.0 else flip sign of theta2 - - double byzvec1[3],quat2[4]; - MathExtra::copy3(bvec1,byzvec1); - byzvec1[0] = 0.0; - double byzvec1_len = MathExtra::len3(byzvec1); - double theta2 = acos(bvec1[1]/byzvec1_len); - if (bvec1[2] > 0.0) theta2 = -theta2; - MathExtra::axisangle_to_quat(xaxis,theta2,quat2); - - // quat_single = rotation via single quat = quat2 * quat1 - // quat_r2g = rotation from restricted to general - // rotate_g2r = general to restricted rotation matrix - // include flip of C vector if needed to obey right-hand rule - // rotate_r2g = restricted to general rotation matrix - // simply a transpose of rotate_g2r since orthonormal - - double quat_single[4]; - MathExtra::quatquat(quat2,quat1,quat_single); - MathExtra::quat_to_mat(quat_single,rotate_g2r); - MathExtra::transpose3(rotate_g2r,rotate_r2g); - - // rotate general ABC to restricted triclinic A'B'C' + // rotate_g2r = rotation matrix from general to restricted triclnic + // rotate_r2g = rotation matrix from restricted to general triclnic 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); + general_to_restricted_rotation(avec,bvec,cvec,rotate_g2r,aprime,bprime,cprime); + MathExtra::transpose3(rotate_g2r,rotate_r2g); // set restricted triclinic boxlo, boxhi, and tilt factors + boxlo[0] = origin_caller[0]; + boxlo[1] = origin_caller[1]; + boxlo[2] = origin_caller[2]; + boxhi[0] = boxlo[0] + aprime[0]; boxhi[1] = boxlo[1] + bprime[1]; boxhi[2] = boxlo[2] + cprime[2]; @@ -663,6 +599,7 @@ void Domain::define_general_triclinic(double *avec_caller, double *bvec_caller, // debug + /* printf("Theta1: %g\n",theta1); printf("Rotvec1: %g %g %g\n",rot1[0],rot1[1],rot1[2]); printf("Theta2: %g\n",theta2); @@ -680,6 +617,85 @@ void Domain::define_general_triclinic(double *avec_caller, double *bvec_caller, 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)); + */ +} + +/* ---------------------------------------------------------------------- + compute rotation matrix to transform from general to restricted triclinic + ABC = 3 general triclinic edge vectors + rotmat = rotation matrix + A`B`C` = 3 restricited triclinic edge vectors +------------------------------------------------------------------------- */ + +void Domain::general_to_restricted_rotation(double *a, double *b, double *c, + double rotmat[3][3], + double *aprime, double *bprime, double *cprime) +{ + // error checks + // A,B,C cannot be co-planar + // A x B must point in C direction (right-handed) + + double abcross[3]; + MathExtra::cross3(a,b,abcross); + double dot = MathExtra::dot3(abcross,c); + if (dot == 0.0) + error->all(FLERR,"General triclinic edge vectors are co-planar"); + if (dot < 0.0) + error->all(FLERR,"General triclinic edge vectors must be right-handed"); + + // 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 alen = MathExtra::len3(a); + MathExtra::cross3(a,xaxis,rot1); + MathExtra::norm3(rot1); + double theta1 = acos(a[0]/alen); + 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 b1[3]; + MathExtra::matvec(rotmat1,b,b1); + + // quat2 = rotation to convert B1 into B' in xy plane + // Byz1 = projection of B1 into yz plane + // +xaxis = unit vector to rotate B1 around + // theta2 = angle of rotation calculated from + // Byz1 dot yunit = B1y = |Byz1| cos(theta2) + // theta2 via acos() is positive (0 to PI) + // positive is valid if B1z < 0.0 else flip sign of theta2 + + double byzvec1[3],quat2[4]; + MathExtra::copy3(b1,byzvec1); + byzvec1[0] = 0.0; + double byzvec1_len = MathExtra::len3(byzvec1); + double theta2 = acos(b1[1]/byzvec1_len); + if (b1[2] > 0.0) theta2 = -theta2; + MathExtra::axisangle_to_quat(xaxis,theta2,quat2); + + // quat_single = rotation via single quat = quat2 * quat1 + // quat_r2g = rotation from restricted to general + // rotmat = general to restricted rotation matrix + + double quat_single[4]; + MathExtra::quatquat(quat2,quat1,quat_single); + MathExtra::quat_to_mat(quat_single,rotmat); + + // rotate general ABC to restricted triclinic A'B'C' + + MathExtra::matvec(rotate_g2r,a,aprime); + MathExtra::matvec(rotate_g2r,b,bprime); + MathExtra::matvec(rotate_g2r,c,cprime); } /* ---------------------------------------------------------------------- diff --git a/src/domain.h b/src/domain.h index 41994e2140..c54c08df3d 100644 --- a/src/domain.h +++ b/src/domain.h @@ -142,6 +142,9 @@ class Domain : protected Pointers { int ownatom(int, double *, imageint *, int); void define_general_triclinic(double *, double *, double *, double *); + void general_to_restricted_rotation(double *, double *, double *, + double [3][3], + double *, double *, double *); void general_to_restricted_coords(double *); void restricted_to_general_coords(double *); void restricted_to_general_coords(double *, double *); diff --git a/src/lattice.cpp b/src/lattice.cpp index 815700003d..d29a378589 100644 --- a/src/lattice.cpp +++ b/src/lattice.cpp @@ -138,6 +138,8 @@ Lattice::Lattice(LAMMPS *lmp, int narg, char **arg) : Pointers(lmp) // process optional args + int triclinic_general = 0; + int iarg = 2; while (iarg < narg) { if (strcmp(arg[iarg],"origin") == 0) { @@ -222,6 +224,11 @@ Lattice::Lattice(LAMMPS *lmp, int narg, char **arg) : Pointers(lmp) error->all(FLERR, "Invalid lattice basis argument: {}", z); add_basis(x,y,z); iarg += 4; + + } else if (strcmp(arg[iarg],"triclinic/general") == 0) { + triclinic_general = 1; + iarg++; + } else error->all(FLERR,"Unknown lattice keyword: {}", arg[iarg]); } @@ -256,6 +263,26 @@ 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() + // a123 prime are used to compute lattice spacings + + if (triclinic_general) { + if (style != CUSTOM) + error->all(FLERR,"Lattice triclnic/general must be style = CUSTOM"); + if (origin[0] != 0.0 || origin[1] != 0.0 || origin[2] != 0.0) + error->all(FLERR,"Lattice triclnic/general must have default origin"); + int oriented = 0; + if (orientx[0] != 1 || orientx[1] != 0 || orientx[2] != 0) oriented = 1; + if (orienty[0] != 0 || orienty[1] != 1 || orienty[2] != 0) oriented = 1; + if (orientz[0] != 0 || orientz[1] != 0 || orientz[2] != 1) oriented = 1; + if (oriented) + error->all(FLERR,"Lattice triclnic/general must have default orientation"); + + double rotmat[3][3]; + domain->general_to_restricted_rotation(a1,a2,a3,rotmat,a1_prime,a2_prime,a3_prime); + } + // reset scale for LJ units (input scale is rho*) // scale = (Nbasis/(Vprimitive * rho*)) ^ (1/dim) // use fabs() in case a1,a2,a3 are not right-handed for general triclinic @@ -267,18 +294,11 @@ Lattice::Lattice(LAMMPS *lmp, int narg, char **arg) : Pointers(lmp) scale = pow(nbasis/volume/scale,1.0/dimension); } - // set orientflag - // for general triclinic, create_box and create_atoms require orientflag = 0 - - oriented = 0; - if (orientx[0] != 1 || orientx[1] != 0 || orientx[2] != 0) oriented = 1; - if (orienty[0] != 0 || orienty[1] != 1 || orienty[2] != 0) oriented = 1; - if (orientz[0] != 0 || orientz[1] != 0 || orientz[2] != 1) oriented = 1; - // initialize lattice <-> box transformation matrices - setup_transform(); + setup_transform(a1,a2,a3); + // automatic calculation of lattice spacings // convert 8 corners of primitive unit cell from lattice coords to box coords // min to max = bounding box around the pts in box space // xlattice,ylattice,zlattice = extent of bbox in box space @@ -291,6 +311,14 @@ Lattice::Lattice(LAMMPS *lmp, int narg, char **arg) : Pointers(lmp) xmax = ymax = zmax = -BIG; xlattice = ylattice = zlattice = 0.0; + // for general triclinic, bounding box is around unit cell + // in restricted triclinic orientation, NOT general + // this enables lattice spacings to be used for other commands (e.g. region) + // after create_box and create_atoms create the restricted triclnic system + // reset transform used by bbox() to be based on rotated a123 prime vectors + + if (triclinic_general) setup_transform(a1_prime,a2_prime,a3_prime); + bbox(0,0.0,0.0,0.0,xmin,ymin,zmin,xmax,ymax,zmax); bbox(0,1.0,0.0,0.0,xmin,ymin,zmin,xmax,ymax,zmax); bbox(0,0.0,1.0,0.0,xmin,ymin,zmin,xmax,ymax,zmax); @@ -300,10 +328,16 @@ Lattice::Lattice(LAMMPS *lmp, int narg, char **arg) : Pointers(lmp) bbox(0,0.0,1.0,1.0,xmin,ymin,zmin,xmax,ymax,zmax); bbox(0,1.0,1.0,1.0,xmin,ymin,zmin,xmax,ymax,zmax); + // restore original general triclinic a123 transform + + if (triclinic_general) setup_transform(a2,a2,a3); + xlattice = xmax - xmin; ylattice = ymax - ymin; zlattice = zmax - zmin; + // user-defined lattice spacings + } else { xlattice *= scale; ylattice *= scale; @@ -325,24 +359,13 @@ Lattice::~Lattice() } /* ---------------------------------------------------------------------- - return 1 if style = CUSTOM, else 0 - queried by create_box and create_atoms for general triclinic + return 1 if lattice is for a general triclinic simulation box + queried by create_box and create_atoms ------------------------------------------------------------------------- */ -int Lattice::is_custom() +int Lattice::is_general_triclinic() { - if (style == CUSTOM) return 1; - return 0; -} - -/* ---------------------------------------------------------------------- - return 1 any orient vectors are non-default, else 0 - queried by create_box and create_atoms for general triclinic -------------------------------------------------------------------------- */ - -int Lattice::is_oriented() -{ - if (oriented) return 1; + if (triclinic_general) return 1; return 0; } @@ -395,21 +418,21 @@ int Lattice::collinear() initialize lattice <-> box transformation matrices ------------------------------------------------------------------------- */ -void Lattice::setup_transform() +void Lattice::setup_transform(double *a, double *b, double *c) { double length; // primitive = 3x3 matrix with primitive vectors as columns - primitive[0][0] = a1[0]; - primitive[1][0] = a1[1]; - primitive[2][0] = a1[2]; - primitive[0][1] = a2[0]; - primitive[1][1] = a2[1]; - primitive[2][1] = a2[2]; - primitive[0][2] = a3[0]; - primitive[1][2] = a3[1]; - primitive[2][2] = a3[2]; + primitive[0][0] = a[0]; + primitive[1][0] = a[1]; + primitive[2][0] = a[2]; + primitive[0][1] = b[0]; + primitive[1][1] = b[1]; + primitive[2][1] = b[2]; + primitive[0][2] = c[0]; + primitive[1][2] = c[1]; + primitive[2][2] = c[2]; // priminv = inverse of primitive diff --git a/src/lattice.h b/src/lattice.h index d91b5cc834..74f213bdb4 100644 --- a/src/lattice.h +++ b/src/lattice.h @@ -28,34 +28,37 @@ class Lattice : protected Pointers { int nbasis; // # of basis atoms in unit cell double **basis; // fractional coords of each basis atom // within unit cell (0 <= coord < 1) - + Lattice(class LAMMPS *, int, char **); ~Lattice() override; void lattice2box(double &, double &, double &); void box2lattice(double &, double &, double &); void bbox(int, double, double, double, double &, double &, double &, double &, double &, double &); - int is_custom(); - int is_oriented(); + int is_general_triclinic(); private: + int triclinic_general; // 1 if general triclinic, else 0 + int oriented; // 1 if non-default orient xyz, else 0 double scale; - double origin[3]; // lattice origin - int oriented; // 1 if non-default orient xyz, else 0 - int orientx[3]; // lattice orientation vecs - int orienty[3]; // orientx = what lattice dir lies - int orientz[3]; // along x dim in box - + double origin[3]; // lattice origin + int orientx[3]; // lattice orientation vecs + int orienty[3]; // orientx = what lattice dir lies + int orientz[3]; // along x dim in box - double primitive[3][3]; // lattice <-> box transform matrices + double primitive[3][3]; // lattice <-> box transformation matrices double priminv[3][3]; double rotaterow[3][3]; double rotatecol[3][3]; + double a1_prime[3]; // a123 rotated to restricted triclinic orientation + double a2_prime[3]; + double a3_prime[3]; + int orthogonal(); int right_handed(); int collinear(); - void setup_transform(); + void setup_transform(double *, double *, double *); void add_basis(double, double, double); double dot(double *, double *); void cross(double *, double *, double *);