"""mBuild box module."""
from warnings import warn
import numpy as np
from mbuild.exceptions import MBuildError
__all__ = ["Box"]
[docs]class Box(object):
"""A box representing the bounds of the system.
Parameters
----------
lengths : list-like, shape=(3,), dtype=float
Lengths of the edges of the box.
angles : list-like, shape=(3,), dtype=float, default=None
Angles (in degrees) that define the tilt of the edges of the box. If
None is given, angles are assumed to be [90.0, 90.0, 90.0]. These are
also known as alpha, beta, gamma in the crystallography community.
precision : int, optional, default=None
Control the precision of the floating point representation of box
attributes. If none provided, the default is 6 decimals.
Attributes
----------
vectors : np.ndarray, shape=(3,3), dtype=float
Vectors that define the parallelepiped (Box).
lengths : tuple, shape=(3,), dtype=float
Lengths of the box in x,y,z
angles : tuple, shape=(3,), dtype=float
Angles defining the tilt of the box.
Lx : float
Length of the Box in the x dimension
Ly : float
Length of the Box in the y dimension
Lz : float
Length of the Box in the z dimension
xy : float
Tilt factor needed to displace an orthogonal box's xy face to its
parallelepiped structure.
xz : float
Tilt factor needed to displace an orthogonal box's xz face to its
parallelepiped structure.
yz : float
Tilt factor needed to displace an orthogonal box's yz face to its
parallelepiped structure.
precision : int
Precision of the floating point numbers when accessing values.
Notes
-----
Box vectors are expected to be provided in row-major format.
"""
def __init__(self, lengths, angles=None, precision=None):
if precision is not None:
self._precision = int(precision)
else:
self._precision = 6
if angles is None:
angles = [90.0, 90.0, 90.0]
self._vectors = _lengths_angles_to_vectors(
lengths=lengths, angles=angles, precision=self.precision
)
(Lx, Ly, Lz, xy, xz, yz) = self._from_vecs_to_lengths_tilt_factors()
self._Lx = Lx
self._Ly = Ly
self._Lz = Lz
self._xy = xy
self._xz = xz
self._yz = yz
[docs] @classmethod
def from_lengths_angles(cls, lengths, angles, precision=None):
"""Generate a box from lengths and angles."""
return cls(lengths=lengths, angles=angles, precision=precision)
[docs] @classmethod
def from_uvec_lengths(cls, uvec, lengths, precision=None):
"""Generate a box from unit vectors and lengths."""
uvec = np.asarray(uvec)
uvec.reshape(3, 3)
if not np.allclose(np.linalg.norm(uvec, axis=1), 1.0):
raise MBuildError(
"Unit vector magnitudes provided are not close to 1.0, "
f"magnitudes: {np.linalg.norm(uvec, axis=1)}"
)
lengths = np.asarray(lengths)
lengths.reshape(1, 3)
_validate_box_vectors(uvec)
scaled_vec = (uvec.T * lengths).T
(alpha, beta, gamma) = _calc_angles(scaled_vec)
return cls(
lengths=lengths, angles=(alpha, beta, gamma), precision=precision
)
[docs] @classmethod
def from_mins_maxs_angles(cls, mins, maxs, angles, precision=None):
"""Generate a box from min/max distance calculations and angles."""
(x_min, y_min, z_min) = mins
(x_max, y_max, z_max) = maxs
lengths = (x_max - x_min, y_max - y_min, z_max - z_min)
return cls(lengths=lengths, angles=angles, precision=precision)
[docs] @classmethod
def from_vectors(cls, vectors, precision=None):
"""Generate a box from box vectors."""
vectors = _validate_box_vectors(vectors)
(alpha, beta, gamma) = _calc_angles(vectors)
v1 = vectors[0, :]
v2 = vectors[1, :]
v3 = vectors[2, :]
Lx = np.linalg.norm(v1)
Ly = np.linalg.norm(v2)
Lz = np.linalg.norm(v3)
lengths = (Lx, Ly, Lz)
return cls(
lengths=lengths, angles=(alpha, beta, gamma), precision=precision
)
[docs] @classmethod
def from_lengths_tilt_factors(
cls, lengths, tilt_factors=None, precision=None
):
"""Generate a box from box lengths and tilt factors."""
(Lx, Ly, Lz) = lengths
if tilt_factors is None:
(xy, xz, yz) = (0.0, 0.0, 0.0)
else:
(xy, xz, yz) = tilt_factors
vecs = np.asarray(
[[Lx, 0.0, 0.0], [Ly * xy, Ly, 0.0], [Lz * xz, Lz * yz, Lz]]
)
(alpha, beta, gamma) = _calc_angles(vecs)
return cls(
lengths=lengths, angles=[alpha, beta, gamma], precision=precision
)
[docs] @classmethod
def from_lo_hi_tilt_factors(cls, lo, hi, tilt_factors, precision=None):
"""Generate a box from a lo, hi convention and tilt factors."""
(xlo, ylo, zlo) = lo
(xhi, yhi, zhi) = hi
(xy, xz, yz) = tilt_factors
xlo_bound = xlo + min([0.0, xy, xz, xy + xz])
xhi_bound = xhi + max([0.0, xy, xz, xy + xz])
ylo_bound = ylo + min([0.0, yz])
yhi_bound = yhi + max([0.0, yz])
lengths = [xhi_bound - xlo_bound, yhi_bound - ylo_bound, zhi - zlo]
return cls.from_lengths_tilt_factors(
lengths=lengths, tilt_factors=tilt_factors
)
@property
def vectors(self):
"""Box representation as a 3x3 matrix."""
return self._vectors
@property
def box_parameters(self):
"""Lengths and tilt factors of the box."""
return self.Lx, self.Ly, self.Lz, self.xy, self.xz, self.xy
@property
def Lx(self):
"""Length in the x direction."""
return round(self._Lx, self.precision)
@property
def Ly(self):
"""Length in the y direction."""
return round(self._Ly, self.precision)
@property
def Lz(self):
"""Length in the z direction."""
return round(self._Lz, self.precision)
@property
def lengths(self):
"""Lengths of the box."""
return self.Lx, self.Ly, self.Lz
@property
def xy(self):
"""Tilt factor xy of the box."""
return round(self._xy, self.precision)
@property
def xz(self):
"""Tilt factor xz of the box."""
return round(self._xz, self.precision)
@property
def yz(self):
"""Tilt factor yz of the box."""
return round(self._yz, self.precision)
@property
def tilt_factors(self):
"""Return the 3 tilt_factors (xy, xz, yz) of the box."""
return self.xy, self.xz, self.yz
@property
def angles(self):
"""Angles defining the tilt of the box (alpha, beta, gamma)."""
(alpha, beta, gamma) = self._get_angles()
alpha = round(alpha, self.precision)
beta = round(beta, self.precision)
gamma = round(gamma, self.precision)
return alpha, beta, gamma
@property
def precision(self):
"""Amount of decimals to represent floating point values."""
return self._precision
@precision.setter
def precision(self, value):
"""Decimal point precision, if None use 16, else cast as int."""
if not value:
precision = 16
else:
precision = int(value)
self._precision = precision
@property
def bravais_parameters(self):
"""Return the Box representation as Bravais lattice parameters.
Based on the box vectors, return the parameters to describe the box in
terms of the Bravais lattice parameters:
a,b,c = the edges of the Box
alpha, beta, gamma = angles(tilt) of the parallelepiped, in degrees
Returns
-------
parameters : tuple of floats,
(a, b, c, alpha, beta, gamma)
"""
(alpha, beta, gamma) = self.angles
(Lx, Ly, Lz) = self.lengths
return Lx, Ly, Lz, alpha, beta, gamma
def __repr__(self):
"""Return a string representation of the box."""
(Lx, Ly, Lz, xy, xz, yz) = self.box_parameters
format_precision = f".{self._precision}f" if self._precision else ""
desc = (
f"Box: Lx={Lx:{format_precision}}, "
f"Ly={Ly:{format_precision}}, "
f"Lz={Lz:{format_precision}}, "
f"xy={xy:{format_precision}}, "
f"xz={xz:{format_precision}}, "
f"yz={yz:{format_precision}}, "
)
return desc
def _from_vecs_to_lengths_tilt_factors(self):
# vectors should already be aligned by _normalize_box
v = np.zeros((3, 3))
v[0, :] = self._vectors[0, :]
v[1, :] = self._vectors[1, :]
v[2, :] = self._vectors[2, :]
Lx = np.sqrt(np.dot(v[0], v[0]))
a2x = np.dot(v[0], v[1]) / Lx
Ly = np.sqrt(np.dot(v[1], v[1]) - a2x * a2x)
xy = a2x / Ly
v0xv1 = np.cross(v[0], v[1])
v0xv1mag = np.sqrt(np.dot(v0xv1, v0xv1))
Lz = np.dot(v[2], v0xv1) / v0xv1mag
a3x = np.dot(v[0], v[2]) / Lx
xz = a3x / Lz
yz = (np.dot(v[1], v[2]) - a2x * a3x) / (Ly * Lz)
len_x = np.sqrt(np.dot(v[0], v[0]))
len_y = np.sqrt(np.dot(v[1], v[1]))
len_z = np.sqrt(np.dot(v[2], v[2]))
return len_x, len_y, len_z, xy, xz, yz
def _get_angles(self):
return _calc_angles(self.vectors)
def _validate_box_vectors(box_vectors):
"""Determine if the vectors are in the convention we use.
This method will parse the provided box vectors, determine if the vectors
follow the conventions the Box class adheres to. In this case:
1. It is a 3x3 matrix that can be coerced into a numpy array of floats
2. Vectors are in a right-handed basis (determinant of matrix is +)
3. The first vector is aligned along the [1,0,0] direction
4. The second vector in aligned along the xy plane
5. The third vector can align freely in the x,y, and +z direction
If the vectors are not right-handed, a warning will be raised, and the
vectors will be transformed into the right-handed coordinate system.
If the three vectors are not following conventions 3-5, the matrix will be
transformed to comply with them, and will also raise a warning.
"""
vecs = np.asarray(box_vectors, dtype=np.float64)
vecs.reshape(3, 3)
return _normalize_box(vecs)
def _lengths_angles_to_vectors(lengths, angles, precision):
(a, b, c) = lengths
(alpha, beta, gamma) = np.deg2rad(angles)
cos_a = np.clip(np.cos(alpha), -1.0, 1.0)
cos_b = np.clip(np.cos(beta), -1.0, 1.0)
cos_g = np.clip(np.cos(gamma), -1.0, 1.0)
sin_a = np.clip(np.sin(alpha), -1.0, 1.0)
sin_b = np.clip(np.sin(beta), -1.0, 1.0)
sin_g = np.clip(np.sin(gamma), -1.0, 1.0)
a_vec = np.asarray([a, 0.0, 0.0])
b_x = b * cos_g
b_y = b * sin_g
b_vec = np.asarray([b_x, b_y, 0.0])
c_x = c * cos_b
c_cos_y_term = (cos_a - (cos_b * cos_g)) / sin_g
c_y = c * c_cos_y_term
c_z = c * np.sqrt(1 - np.square(cos_b) - np.square(c_cos_y_term))
c_vec = np.asarray([c_x, c_y, c_z])
box_vectors = np.asarray((a_vec, b_vec, c_vec))
box_vectors.reshape(3, 3)
_validate_box_vectors(box_vectors=box_vectors)
return box_vectors.round(precision)
def _normalize_box(vectors):
"""Align the box matrix into a right-handed coordinate frame.
NOTE: This assumes that the matrix is in a row-major format.
NOTE: Inspiration and logic are from the Glotzer group package, Garnett;
which is provided under a BSD 3-clause License.
For additional information, refer to the License file provided with this
package.
"""
det = np.linalg.det(vectors)
if np.isclose(det, 0.0, atol=1e-5):
raise MBuildError(
"The vectors to define the box are co-linear, this does not form a "
f"3D region in space.\n Box vectors evaluated: {vectors}"
)
if det < 0.0:
warn(
"Box vectors provided for a left-handed basis, these will be "
"transformed into a right-handed basis automatically."
)
# transpose to column-major for the time being
Q, R = np.linalg.qr(vectors.T)
# left or right handed: det<0 left, >0, right
sign = np.linalg.det(Q)
R = R * sign
signs = np.diag(
np.diag(np.where(R < 0, -np.ones(R.shape), np.ones(R.shape)))
)
transformed_vecs = R.dot(signs)
return _reduced_form_vectors(transformed_vecs.T)
def _reduced_form_vectors(box_vectors):
"""Get reduced vectors from vectors.
Adapted from HOOMD-Blue's documentation on periodic boundary conditions:
https://hoomd-blue.readthedocs.io/en/stable/box.html
"""
v1 = box_vectors[0, :]
v2 = box_vectors[1, :]
v3 = box_vectors[2, :]
lx = np.linalg.norm(v1)
a_2x = np.dot(v1, v2) / lx
ly = np.sqrt(np.dot(v2, v2) - a_2x * a_2x)
xy = a_2x / ly
v1_x_v2 = np.cross(v1, v2)
lz = np.dot(v3, (v1_x_v2 / np.linalg.norm(v1_x_v2)))
a_3x = np.dot(v1, v3) / lx
xz = a_3x / lz
yz = (np.dot(v2, v3) - a_2x * a_3x) / (ly * lz)
reduced_vecs = np.asarray(
[[lx, 0.0, 0.0], [xy * ly, ly, 0.0], [xz * lz, yz * lz, lz]]
)
return reduced_vecs
def _calc_angles(vectors):
"""Calculate the angles between the vectors that define the box.
Calculates the angles alpha, beta, and gamma from the Box object
attribute box_vectors, rounded to 'precision' number of decimal points.
"""
vector_magnitudes = np.linalg.norm(vectors, axis=1)
v = np.zeros((3, 3))
v[0, :] = vectors[0, :]
v[1, :] = vectors[1, :]
v[2, :] = vectors[2, :]
a_dot_b = np.dot(vectors[0], vectors[1])
b_dot_c = np.dot(vectors[1], vectors[2])
a_dot_c = np.dot(vectors[0], vectors[2])
alpha_raw = b_dot_c / (vector_magnitudes[1] * vector_magnitudes[2])
beta_raw = a_dot_c / (vector_magnitudes[0] * vector_magnitudes[2])
gamma_raw = a_dot_b / (vector_magnitudes[0] * vector_magnitudes[1])
(alpha, beta, gamma) = np.rad2deg(
np.arccos(np.clip([alpha_raw, beta_raw, gamma_raw], -1.0, 1.0))
)
return alpha, beta, gamma