ICODE-ADC/lammps
4
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

alter how lattice interacts with create_box and create_atoms for general triclinic

Browse Source
This commit is contained in:
Steve Plimpton
2023-11-17 17:10:45 -07:00
parent 4ccd59af80
commit 56b2c7ed46
11 changed files with 383 additions and 206 deletions
Download Patch File Download Diff File Expand all files Collapse all files

6
doc/src/Howto_triclinic.rst
View File

@ -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

112
doc/src/create_atoms.rst
View File

@ -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 <read_data>` or :doc:`read_restart
<read_restart>` command. A simulation box must already exist, which is
typically created via the :doc:`create_box <create_box>` command.
Before using this command, a lattice must also be defined using the
:doc:`lattice <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
<read_data>` or :doc:`read_restart <read_restart>` command. A
simulation box must already exist, which is typically created via the
:doc:`create_box <create_box>` command. Before using this command, a
lattice must typically also be defined using the :doc:`lattice
<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 <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 <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 <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 <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 <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
<create_atoms>` command are converted from general to restricted
triclinic form when the two commands are invoked. The
<Howto_triclinic>` 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 <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 <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.
----------

56
doc/src/create_box.rst
View File

@ -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 <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 <lattice>` command.
previously defind by the :doc:`lattice <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 <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 <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 <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 <lattice>` command default for a3 = (0,0,1).
.. note::

68
doc/src/lattice.rst
View File

@ -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
<Howto_triclnic>` 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 <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 <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 <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 <create_box>` and
:doc:`create_atoms <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 <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

2
src/atom_vec_line.cpp
View File

@ -374,7 +374,7 @@ void AtomVecLine::data_atom_bonus(int m, const std::vector<std::string> &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

52
src/create_atoms.cpp
View File

@ -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
// verify lattice was defined with triclinic/general option
int triclinic_general = domain->triclinic_general;
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];

12
src/create_box.cpp
View File

@ -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

160
src/domain.cpp
View File

@ -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);
}
/* ----------------------------------------------------------------------

3
src/domain.h
View File

@ -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 *);

91
src/lattice.cpp
View File

@ -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

15
src/lattice.h
View File

@ -35,27 +35,30 @@ class Lattice : protected Pointers {
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 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 *);