Introduction to atomman: Box class

Lucas M. Hale, lucas.hale@nist.gov, Materials Science and Engineering Division, NIST.

Disclaimers

1. Introduction

The atomman.Box class represents a generic parallelepiped that can be used as a simulation box/cell allowing for translational symmetry in all three dimensions. The underlying data structure consists of three 3D vectors:

  • avect is the Cartesian vector associated with the Box’s a lattice vector.

  • bvect is the Cartesian vector associated with the Box’s b lattice vector.

  • cvect is the Cartesian vector associated with the Box’s c lattice vector.

  • origin is the Cartesian position vector for the Box’s origin. The positions of the Box’s eight corners are given by adding combinations of avect, bvect and cvect to the origin.

Library Imports

[1]:
# Standard Python libraries
import datetime

# http://www.numpy.org/
import numpy as np

# https://github.com/usnistgov/atomman
import atomman as am
import atomman.unitconvert as uc

# Show atomman version
print('atomman version =', am.__version__)

# Show date of Notebook execution
print('Notebook executed on', datetime.date.today())
atomman version = 1.4.10
Notebook executed on 2023-07-28

Create a default box with vects that are normal unit vectors.

[2]:
box = am.Box()

print(box)
avect =  [ 1.000,  0.000,  0.000]
bvect =  [ 0.000,  1.000,  0.000]
cvect =  [ 0.000,  0.000,  1.000]
origin = [ 0.000,  0.000,  0.000]

2. Box parameters

The Box class supports a number of different parameters to represent the underlying box. All of these parameters can be fetched as class attributes.

  • avect, bvect, cvect are the Cartesian vectors associated with the Box’s lattice vectors.

  • origin is the Cartesian position vector for the Box’s origin.

  • vects is a 3x3 array containing [avect, bvect, cvect].

  • a, b, c are the Box’s lattice parameters (magnitudes of avect, bvect, cvect, respectively).

  • alpha, beta, gamma are the Box’s lattice angles in degrees.

  • xlo, xhi, ylo, yhi, zlo, zhi are the LAMMPS hi/lo box dimensions.

  • lx, ly, lz are the LAMMPS box dimensions (lx = xhi - xlo, etc.)

  • xy, xz, yz are the LAMMPS box tilt factors.

  • volume is the box’s volume (added version 1.2.5).

[3]:
# Individual box vectors
print('box.avect ->', box.avect)
print('box.bvect ->', box.bvect)
print('box.cvect ->', box.cvect)
print()

# Box origin position
print('box.origin ->', box.origin)
print()

# Array of box vectors
print('box.vects ->')
print(box.vects)
print()

# Crystal lattice parameters
print('box.a ->', box.a)
print('box.b ->', box.b)
print('box.c ->', box.c)
print('box.alpha ->', box.alpha)
print('box.beta  ->', box.beta)
print('box.gamma ->', box.gamma)
print()

# LAMMPS parameters
print('box.xlo ->', box.xlo)
print('box.xhi ->', box.xhi)
print('box.ylo ->', box.ylo)
print('box.yhi ->', box.yhi)
print('box.zlo ->', box.zlo)
print('box.zhi ->', box.zhi)
print('box.lx ->', box.lx)
print('box.ly ->', box.ly)
print('box.lz ->', box.lz)
print('box.xy ->', box.xy)
print('box.xz ->', box.xz)
print('box.yz ->', box.yz)
print()

# Box volume
print('box.volume ->', box.volume)
box.avect -> [1. 0. 0.]
box.bvect -> [0. 1. 0.]
box.cvect -> [0. 0. 1.]

box.origin -> [0. 0. 0.]

box.vects ->
[[1. 0. 0.]
 [0. 1. 0.]
 [0. 0. 1.]]

box.a -> 1.0
box.b -> 1.0
box.c -> 1.0
box.alpha -> 90.0
box.beta  -> 90.0
box.gamma -> 90.0

box.xlo -> 0.0
box.xhi -> 1.0
box.ylo -> 0.0
box.yhi -> 1.0
box.zlo -> 0.0
box.zhi -> 1.0
box.lx -> 1.0
box.ly -> 1.0
box.lz -> 1.0
box.xy -> 0.0
box.xz -> 0.0
box.yz -> 0.0

box.volume -> 1.0

3. Initializing and setting Box parameters

Only a few of the parameters listed in Section #2 can be directly set. This is done as setting some of the parameters independent of others can lead to ambiguous answers. For better control, set() functions are defined that allow for the setting of complete parameter sets for defining the box.

3.1. Direct setting box vects and origin

Only vects and origin can be directly set.

[4]:
# Set avect, bvect, cvect and origin with random vectors
box.vects = 4 * np.random.rand(3,3)
box.origin = np.random.rand(3)
print(box)
avect =  [ 0.758,  0.421,  3.528]
bvect =  [ 1.131,  2.927,  0.314]
cvect =  [ 0.966,  2.736,  3.084]
origin = [ 0.540,  0.499,  0.814]

3.2. Box initializing and Box.set() method

A new Box can be initialized or an existing Box can be changed using the Box.set() method by using one of the following sets of parameters (optional terms in parenthesis):

  • vects (and origin)

  • avect, bvect, cvect (and origin)

  • a, b, c, (alpha, beta, gamma and origin)

  • xlo, xhi, ylo, yhi, zlo, zhi, (xy, xz and yz)

  • lx, ly, lz, (xy, xz, yz, and origin)

[5]:
# Use set with vects (default origin is [0,0,0])
box.set(vects=[[1,2,3], [2,3,1], [3,1,2]])
print(box)
avect =  [ 1.000,  2.000,  3.000]
bvect =  [ 2.000,  3.000,  1.000]
cvect =  [ 3.000,  1.000,  2.000]
origin = [ 0.000,  0.000,  0.000]
[6]:
# Use set with avect, bvect, cvect (default origin is [0,0,0])
box.set(avect=[3.2, 0.0, 0.0], bvect=[0.0, 3.2, 0.0], cvect=[0.0, 0.0, 3.2])
print(box)
avect =  [ 3.200,  0.000,  0.000]
bvect =  [ 0.000,  3.200,  0.000]
cvect =  [ 0.000,  0.000,  3.200]
origin = [ 0.000,  0.000,  0.000]
[7]:
# Use set with a, b, c and alpha (default angles are 90, origin is [0,0,0])
box.set(a=4.3, b=3.2, c=8.1, alpha=110)
print(box)
print()

# Show that box definition coincides with parameters set
print('box.a ->', box.a)
print('box.b ->', box.b)
print('box.c ->', box.c)
print('box.alpha ->', box.alpha)
print('box.beta  ->', box.beta)
print('box.gamma ->', box.gamma)
avect =  [ 4.300,  0.000,  0.000]
bvect =  [ 0.000,  3.200,  0.000]
cvect =  [ 0.000, -2.770,  7.612]
origin = [ 0.000,  0.000,  0.000]

box.a -> 4.3
box.b -> 3.2
box.c -> 8.1
box.alpha -> 110.0
box.beta  -> 90.0
box.gamma -> 90.0
[8]:
# Use set with xlo, xhi, ylo, yhi, zlo, zhi and xy (default tilts are 0)
box.set(xlo=-1, xhi=5, ylo=-2.1, yhi=5, zlo=0.1, zhi=3.1, xy=0.5)
print(box)
print()

# Show that box definition coincides with parameters set
print('box.xlo ->', box.xlo)
print('box.xhi ->', box.xhi)
print('box.ylo ->', box.ylo)
print('box.yhi ->', box.yhi)
print('box.zlo ->', box.zlo)
print('box.zhi ->', box.zhi)
print('box.xy ->', box.xy)
print('box.xz ->', box.xz)
print('box.yz ->', box.yz)
avect =  [ 6.000,  0.000,  0.000]
bvect =  [ 0.500,  7.100,  0.000]
cvect =  [ 0.000,  0.000,  3.000]
origin = [-1.000, -2.100,  0.100]

box.xlo -> -1.0
box.xhi -> 5.0
box.ylo -> -2.1
box.yhi -> 5.0
box.zlo -> 0.1
box.zhi -> 3.1
box.xy -> 0.5
box.xz -> 0.0
box.yz -> 0.0
[9]:
# Use set with lx, ly, lz and xz (default tilts are 0, origin is [0,0,0])
box.set(lx=42, ly=57, lz=112, xz=15, origin=[1,2,3])
print(box)
print()

# Show that box definition coincides with parameters set
print('box.lx ->', box.lx)
print('box.ly ->', box.ly)
print('box.lz ->', box.lz)
print('box.xy ->', box.xy)
print('box.xz ->', box.xz)
print('box.yz ->', box.yz)
avect =  [42.000,  0.000,  0.000]
bvect =  [ 0.000, 57.000,  0.000]
cvect =  [15.000,  0.000, 112.000]
origin = [ 1.000,  2.000,  3.000]

box.lx -> 42.0
box.ly -> 57.0
box.lz -> 112.0
box.xy -> 0.0
box.xz -> 15.0
box.yz -> 0.0

3.3. Crystal family static methods

Added version 1.2.5

There are also methods for each of the seven crystal families that are convenient for the construction of unit cell systems in standard representations:

  • cubic(a): \(a = b = c; \alpha = \beta = \gamma = 90^\circ\)

  • hexagonal(a, c): \(a = b \ne c; \alpha = \beta = 90^\circ; \gamma = 120^\circ\)

  • tetragonal(a, c): \(a = b \ne c; \alpha = \beta = \gamma = 90^\circ\)

  • trigonal(a, alpha): \(a = b = c; \alpha = \beta = \gamma < 120^\circ\)

  • orthorhombic(a, b, c): \(a \ne b \ne c; \alpha = \beta = \gamma = 90^\circ\)

  • monoclinic(a, b, c, beta): \(a \ne b \ne c; \alpha = \gamma = 90^\circ; \beta > 90^\circ\)

  • triclinic(a, b, c, alpha, beta, gamma): \(a \ne b \ne c; \alpha \ne \beta \ne \gamma\)

[10]:
print('am.Box.cubic(4.25) ->')
box = am.Box.cubic(4.25)
print('box.a ->', box.a)
print('box.b ->', box.b)
print('box.c ->', box.c)
print('box.alpha ->', box.alpha)
print('box.beta  ->', box.beta)
print('box.gamma ->', box.gamma)
am.Box.cubic(4.25) ->
box.a -> 4.25
box.b -> 4.25
box.c -> 4.25
box.alpha -> 90.0
box.beta  -> 90.0
box.gamma -> 90.0
[11]:
print('am.Box.hexagonal(3.12, 5.14) ->')
box = am.Box.hexagonal(3.12, 5.14)
print('box.a ->', box.a)
print('box.b ->', box.b)
print('box.c ->', box.c)
print('box.alpha ->', box.alpha)
print('box.beta  ->', box.beta)
print('box.gamma ->', box.gamma)
am.Box.hexagonal(3.12, 5.14) ->
box.a -> 3.12
box.b -> 3.12
box.c -> 5.14
box.alpha -> 90.0
box.beta  -> 90.0
box.gamma -> 119.99999999999999

4. Box model

Added version 1.2.7

A JSON/XML equivalent data model representation of the Box object can be generated using the model() method.

The units of length used in the model can be set using the ‘length_unit’ parameter.

[12]:
model = box.model()
print(model.json())
print()
print(model.xml())
{"box": {"avect": {"value": [3.12, 0.0, 0.0], "unit": "angstrom"}, "bvect": {"value": [-1.5599999999999994, 2.701999259807449, 0.0], "unit": "angstrom"}, "cvect": {"value": [0.0, 0.0, 5.14], "unit": "angstrom"}, "origin": {"value": [0.0, 0.0, 0.0], "unit": "angstrom"}}}

<?xml version="1.0" encoding="utf-8"?>
<box><avect><value>3.12</value><value>0.0</value><value>0.0</value><unit>angstrom</unit></avect><bvect><value>-1.5599999999999994</value><value>2.701999259807449</value><value>0.0</value><unit>angstrom</unit></bvect><cvect><value>0.0</value><value>0.0</value><value>5.14</value><unit>angstrom</unit></cvect><origin><value>0.0</value><value>0.0</value><value>0.0</value><unit>angstrom</unit></origin></box>
[13]:
model = box.model(length_unit='nm')
print(model.json())
{"box": {"avect": {"value": [0.312, 0.0, 0.0], "unit": "nm"}, "bvect": {"value": [-0.15599999999999994, 0.27019992598074494, 0.0], "unit": "nm"}, "cvect": {"value": [0.0, 0.0, 0.514], "unit": "nm"}, "origin": {"value": [0.0, 0.0, 0.0], "unit": "nm"}}}

Any stored model information can then be reloaded in as a Box object by passing the ‘model’ parameter to either the class initializer or the model() method.

[14]:
print(am.Box(model=model))
avect =  [ 3.120,  0.000,  0.000]
bvect =  [-1.560,  2.702,  0.000]
cvect =  [ 0.000,  0.000,  5.140]
origin = [ 0.000,  0.000,  0.000]
[15]:
box.model(model=model)
print(box)
avect =  [ 3.120,  0.000,  0.000]
bvect =  [-1.560,  2.702,  0.000]
cvect =  [ 0.000,  0.000,  5.140]
origin = [ 0.000,  0.000,  0.000]

5. Characterization and utility methods

This section describes some additional useful methods built into the Box class.

5.1. LAMMPS-compatible Boxes

For boxes to be compatible with LAMMPS, they have to adhere to the following conditions:

  • avect, bvect, cvect are right-handed.

  • Only certain components of the lattice vectors are allowed to be non-zero:

    avect = [lx, 0.0, 0.0]
    bvect = [xy,  ly, 0.0]
    cvect = [xz,  yz,  lz]
    
  • The tilt factors are limited to being between -1/2 and 1/2 the corresponding length terms.

The first two conditions are automatically adhered to if the box is set with LAMMPS terms or crystal lattice parameters, but may not be true if the box was set using the crystal vectors. The third condition is not checked by atomman. Large tilt factors are allowed by LAMMPS if the “box tilt large” command is used.

5.1.1. Box.is_lammps_norm()

This function returns True if the Box is compatible with the first two LAMMPS condtions and False otherwise.

[16]:
# The current box was defined with LAMMPS parameters, therefore is compatible
print('box.is_lammps_norm() ->', box.is_lammps_norm())
box.is_lammps_norm() -> True
[17]:
# Define a non-lammps compatible box using set(avect, bvect, cvect)
box.set(avect=[10, 0.1, 0.0], bvect=[0.2, 9.0, 0.0], cvect=[-0.2, 0.5, 14])

print('box.is_lammps_norm() ->', box.is_lammps_norm())
box.is_lammps_norm() -> False

Trying to access LAMMPS box parameters for incompatible boxes will issue an error.

[18]:
try:
    box.lx
except AssertionError as e:
    print('AssertionError:', e)
AssertionError: Box is not normalized for LAMMPS style parameters

5.2. Inside/outside

Added version 1.3.0

The Box is treated as a regional shape object (see 3.3._Region_selectors.html) and has methods inside() and outside() that indicate whether any 3D coordinate(s) are located inside or outside the shape.

[19]:
# Build random list of points
points = np.random.rand(20,3) * 10
print(points)
[[0.55155781 2.49380003 9.3967093 ]
 [3.3184798  1.20799191 7.13164885]
 [9.01978498 6.83706011 4.93770814]
 [2.28427281 8.89576691 1.36390599]
 [9.13244568 8.95582225 0.6753787 ]
 [8.82825304 5.60954094 4.37509114]
 [8.72884521 6.66389441 0.27514523]
 [7.71076227 9.11630048 7.75768688]
 [9.48224726 2.40686471 8.43615747]
 [3.19980874 3.80016388 2.24761004]
 [4.48075767 6.73926091 1.42195919]
 [1.61078464 5.56715807 7.12491481]
 [7.72580749 7.54518157 4.28075216]
 [1.18750119 3.26711573 4.83214076]
 [7.90676427 8.19921855 8.48659655]
 [7.98902068 0.54456826 0.9633501 ]
 [4.12189515 2.0876842  8.70759034]
 [5.96837357 2.66878026 8.14396448]
 [9.64936991 5.17531275 8.63342089]
 [6.50309524 6.88621626 1.06736169]]
[20]:
print('Points inside:')
print(box.inside(points))
print('Points outside:')
print(box.outside(points))
Points inside:
[ True  True  True  True  True  True  True  True  True  True  True  True
  True  True  True  True  True  True  True  True]
Points outside:
[False False False False False False False False False False False False
 False False False False False False False False]

5.3. Vector/plane conversions

Added version 1.4.4:

It is often convenient to define positions, vectors and planes relative to a Box. These methods perform the appropriate conversions.

  • vector_crystal_to_cartesian(indices) Converts crystal indices to Cartesian vectors relative to the box’s lattice vectors.

  • plane_crystal_to_cartesian(indices) Converts crystal planar indices to Cartesian plane normal vectors based on the box’s lattice vectors. Note: the algorithm used requires that the planar indices be integers.

  • position_relative_to_cartesian(relpos) Converts position vectors from relative box coordinates to absolute Cartesian coordinates based on the box’s vects and origin.

  • position_cartesian_to_relative(cartpos) Converts position vectors from absolute Cartesian coordinates to relative box coordinates based on the box’s vects and origin.

[21]:
# Define an orthogonal box
a = uc.set_in_units(2.51, 'angstrom')
b = uc.set_in_units(3.13, 'angstrom')
c = uc.set_in_units(4.07, 'angstrom')
box = am.Box(a=a, b=b, c=c)
[22]:
box.vector_crystal_to_cartesian([1,1,1])
[22]:
array([2.51, 3.13, 4.07])
[23]:
box.plane_crystal_to_cartesian([1,1,1])
[23]:
array([0.70300632, 0.56375267, 0.43354935])
[24]:
relpos = [0.5, 0.5, 0.5]
cartpos = box.position_relative_to_cartesian(relpos)
print(cartpos)

relpos2 = box.position_cartesian_to_relative(cartpos)
print(relpos2)
[1.255 1.565 2.035]
[0.5 0.5 0.5]

5.4. Crystal lattice identification

Added version 1.4.4

There are also a few tests for identifying if the box is consistent with a standard representation of a crystal family unit cell.

  • identifyfamily() returns str crystal family if box corresponds to a standard crystal representation. Otherwise, returns None.

  • iscubic() returns bool indicating if box is a standard cubic box.

  • ishexagonal() returns bool indicating if box is a standard hexagonal box.

  • istetragonal() returns bool indicating if box is a standard tetragonal box.

  • isrhombohedral() returns bool indicating if box is a standard rhombohedral box.

  • isorthorhombic() returns bool indicating if box is a standard orthorhombic box.

  • ismonoclinic() returns bool indicating if box is a standard monoclinic box.

  • istriclinic() returns bool indicating if box is a standard triclinic box.

[25]:
print('identifyfamily =', box.identifyfamily())
print('iscubic =       ', box.iscubic())
print('ishexagonal =   ', box.ishexagonal())
print('istetragonal =  ', box.istetragonal())
print('isrhombohedral =', box.isrhombohedral())
print('isorthorhombic =', box.isorthorhombic())
print('ismonoclinic =  ', box.ismonoclinic())
print('istriclinic =   ', box.istriclinic())
identifyfamily = orthorhombic
iscubic =        False
ishexagonal =    False
istetragonal =   False
isrhombohedral = False
isorthorhombic = True
ismonoclinic =   False
istriclinic =    False