Introduction to atomman: Free surface generator
Lucas M. Hale, lucas.hale@nist.gov, Materials Science and Engineering Division, NIST.
1. Introduction
This Notebook outlines the FreeSurface class that allows for the easy generation of planar-sliced free surface atomic configurations based on a unit cell and Miller (hkl) or Miller-Bravais (hkil) crystal plane indices for any crystal unit cell.
New as of version 1.3.2
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
2. Theory
A planar-sliced free surface atomic system can be constructed by
Generating a bulk crystalline system with all periodic boundaries.
Slicing the system along a geometrical plane creating two ideal planar-sliced free surfaces.
While these steps are conceptually simple, it is still tricky to create a routine that can successfully work for any crystal plane in any crystal. The main challenge in developing a generalized free surface routine is working with the multiple different axes systems that describe the defect and the atomic configurations.
Planar free surfaces are typically defined as crystallographic \((hkl)\) Miller or \((hkil)\) Miller-Bravais planes relative to conventional unit cells with lattice box vectors of \(a_{unit}\), \(b_{unit}\), and \(c_{unit}\). Sun and Cedar showed an algorithm that can be easily used to identify the Cartesian slip plane normal for any \((hkl)\) plane in any crystal. In atomman, this is handled by atomman.tools.miller.plane_crystal_to_cartesian.
Using the Cartesian plane normal, a search is then performed to identify three \([uvw]\) crystal vectors corresponding to
The crystal vector with absolute \(h\), \(k\), and \(l\) values less than some search maximum that is closest to the plane normal.
A crystal vector in the slip plane that has the shortest length.
A second crystal vector in the slip plane that has the shortest length without being parallel to the previous vector.
Uniqueness is supported by only picking vectors that form a right-handed system and which have the smallest angle between the two in-plane vectors.
A rotated cell is then constructed in which the rotated cell’s box vectors \(a_{rot}\), \(b_{rot}\), and \(c_{rot}\) correspond to the three \([uvw]\) Miller vectors of the unit cell. The specific alignment of the rotated cell depends on which box vector the free surface cut will be made along, i.e. the “cut box vector”
The cut box vector is aligned with the out-of-plane crystal vector.
The Cartesian system for the rotated cell is defined such that the cut box vector is the only box vector with a component along one of the three x, y, or z Cartesian axes.
This choice means that the free surface can be inserted into the system easily by making the boundary conditions along the cut box vector non-periodic. As the other two box vectors are in-plane vectors, the resulting free surface plane is the plane of interest. Also, the choice of Cartesian system mentioned will result in the plane normal coinciding with one of the three Cartesian axes.
Note: Rotations in atomman create systems that are compatible with the box vector limitations imposed by LAMMPS
\(a = [lx, 0, 0]\)
\(b = [xy, ly, 0]\)
\(c = [xz, yz, lz]\)
As such, \(c\) is the preferred cut box vector as it is the only one with a component in the z direction, and therefore will work for any crystal system.
The resulting free surface system can be constructed from the rotated cell by
Creating a supercell by replicating along the rotated cell box vectors
Rigidly shifting all atoms in the system so that the box edge where the cut will occur falls between two atomic planes
Creating the free surface by either making the cut box vector boundary condition non-periodic, or by increasing the cut box vector to effectively insert a region of vacuum.
3. Class basics
The class is set up to generate a free surface system in two steps:
Initialize the class and specify the (hkl) surface and unit cell structure to use. The proper rotated cell will be generated as well as a list of shifts that can be applied.
Build the free surface system by selecting the shift value corresponding to the termination planes of interest and the size of the final system.
3.1. FreeSurface initialization
hkl (array-like object) The free surface plane to generate expressed in either 3 indices Miller (hkl) format or 4 indices Miller-Bravais (hkil) format.
ucell (atomman.System) The unit cell to use in generating the system.
cutboxvector (str, optional) Specifies which of the three box vectors corresponds to the out-of-plane vector. Default value is c.
maxindex (int, optional) Max uvw index value to use in identifying the best uvw set for the out-of-plane vector. If not given, will use the largest absolute index between the given hkl and the initial in-plane vector guesses.
conventional_setting (str, optional) Allows for rotations of a primitive unit cell to be determined from (hkl) indices specified relative to a conventional unit cell. Allowed settings: ‘p’ for primitive (no conversion), ‘f’ for face-centered, ‘i’ for body-centered, and ‘a’, ‘b’, or ‘c’ for side-centered. Default behavior is to perform no conversion, i.e. take (hkl) relative to the given ucell.
tol (float, optional) Tolerance parameter used to round off near-zero values. Default value is 1e-8.
3.2. FreeSurface.surface()
shift (array-like object, optional) Applies a shift to all atoms. Different values allow for free surfaces with different termination planes to be selected.
vacuumwidth (float, optional) If given, the free surface is created by modifying the system’s box to insert a region of vacuum with this width. This is typically used for DFT calculations where it is computationally preferable to insert a vacuum region and keep all dimensions periodic.
sizemults (list or tuple, optional) The three System.supersize multipliers [a_mult, b_mult, c_mult] to use on the rotated cell to build the final system. Note that the cutboxvector sizemult must be an integer and not a tuple. Default value is [1, 1, 1].
minwidth (float, optional) If given, the sizemult along the cutboxvector will be selected such that the width of the resulting final system in that direction will be at least this value. If both sizemults and minwidth are given, then the larger of the two in the cutboxvector direction will be used.
even (bool, optional) A True value means that the sizemult for cutboxvector will be made an even number by adding 1 if it is odd. Default value is False.
4. Examples
4.1. Conventional bcc unit cell
4.1.1. Define a bcc conventional cell and initialize a FreeSurface
[2]:
a = uc.set_in_units(2.866, 'angstrom')
box = am.Box.cubic(a)
atoms = am.Atoms(pos=[[0.0, 0.0, 0.0], [0.5, 0.5, 0.5]])
ucell = am.System(atoms=atoms, box=box, scale=True, symbols='Fe')
[3]:
plane_hkl = [1, 1, 0]
fs = am.defect.FreeSurface(plane_hkl, ucell)
4.1.2. Check the identified rotated system
The FreeSurface class has a number of attributes that can be used to check the solution.
uvws are the three conventional Miller vectors used as the box vectors of the rotated system.
rcell is the rotated cell.
transform is the Cartesian transformation matrix associated with going from ucell to rcell.
rcellwidth is the width of the rotated cell in the direction normal to the hkl plane.
[4]:
print(fs.uvws)
[[ 0. 0. 1.]
[ 1. -1. 0.]
[ 1. 1. 0.]]
[5]:
print(fs.rcell)
avect = [ 2.866, 0.000, 0.000]
bvect = [ 0.000, 4.053, 0.000]
cvect = [ 0.000, 0.000, 4.053]
origin = [ 0.000, 0.000, 0.000]
natoms = 4
natypes = 1
symbols = ('Fe',)
pbc = [ True True True]
per-atom properties = ['atype', 'pos']
id atype pos[0] pos[1] pos[2]
0 1 1.433 2.027 0.000
1 1 0.000 0.000 0.000
2 1 1.433 0.000 2.027
3 1 0.000 2.027 2.027
[6]:
print(fs.transform)
[[ 0.00000000e+00 4.46536745e-18 1.00000000e+00]
[ 7.07106781e-01 -7.07106781e-01 3.19640147e-16]
[ 7.07106781e-01 7.07106781e-01 -9.73839539e-18]]
[7]:
print(fs.rcellwidth)
4.05313606976129
4.1.3. Build the free surface configuration
Building the final free surface configuration requires specifying an atomic shift to position the fault plane between two atomic planes. The FreeSurface.shifts attribute lists all identified shift positions for rcell that will place the geometrical plane halfway between atomic planes. A shift value is not automatically used as a default as different shifts can result in distinctly different termination planes being created.
[8]:
print(fs.shifts)
[[0. 0. 1.01328402]
[0. 0. 3.03985205]]
For a given shift, a free surface configuration can then be retrieved using the surface method.
[9]:
print(fs.surface(shift=fs.shifts[0], minwidth=10))
avect = [ 2.866, 0.000, 0.000]
bvect = [ 0.000, 4.053, 0.000]
cvect = [ 0.000, 0.000, 12.159]
origin = [ 0.000, 0.000, 0.000]
natoms = 12
natypes = 1
symbols = ('Fe',)
pbc = [ True True False]
per-atom properties = ['atype', 'pos']
id atype pos[0] pos[1] pos[2]
0 1 1.433 2.027 1.013
1 1 0.000 0.000 1.013
2 1 1.433 0.000 3.040
3 1 0.000 2.027 3.040
4 1 1.433 2.027 5.066
5 1 0.000 0.000 5.066
6 1 1.433 0.000 7.093
7 1 0.000 2.027 7.093
8 1 1.433 2.027 9.120
9 1 0.000 0.000 9.120
10 1 1.433 0.000 11.146
11 1 0.000 2.027 11.146
After a free surface configuration has been created, the FreeSurface object will also have a surfacearea attribute which gives the surface area for one of the two free surfaces created.
Note: For free surface formation energy calculations you must divide by 2*surfacearea as two free surfaces are created.
[10]:
print(fs.surfacearea)
11.616287975935858
4.2. Primitive bcc unit cell
The algorithm relies on identifying the smallest in-plane Miller vectors related to the unit cell given. If the conventional unit cell is not primitive, then the identified vectors may not be the smallest lattice vectors possible, and therefore the system may be larger than needed. Alternately, the primitive unit cell can be used to construct the free surface configuration by specifying the conventional_setting parameter
Give ucell as the primitive unit cell.
Use conventional_setting to indicate the Bravais lattie centering: ‘p’ for primitive (no conversion), ‘f’ for face-centered, ‘i’ for body-centered, and ‘a’, ‘b’, or ‘c’ for side-centered.
Then, (hkl) can be specified relative to the conventional unit cell, but used to operate on the primitive unit cell.
4.2.1. Define a primitive bcc unit cell and initialize a FreeSurface
[11]:
box = am.Box.trigonal(a = a * 3**0.5 / 2, alpha = 109.47122063449069)
atoms = am.Atoms(pos=[[0.0, 0.0, 0.0]])
ucell = am.System(atoms=atoms, box=box, scale=True, symbols='Fe')
Note we are giving the same plane_hkl value as above (relative to the conventional cubic cell), but specifying the conventional setting.
[12]:
plane_hkl = [1, 1, 0]
setting = 'i'
fs = am.defect.FreeSurface(plane_hkl, ucell, conventional_setting=setting)
4.1.2. Check the identified rotated system
Using the primitive unit cell allows for the algorithm to find shorter lattice vectors than can be identified with the conventional unit cell. The identified uvws (relative to the conventional cell) may then be fractions instead of integers.
[13]:
print(fs.uvws)
[[-0.5 0.5 -0.5]
[-0.5 0.5 0.5]
[ 1. 1. 0. ]]
For rcell, the system is half the size but still has the same width.
[14]:
print(fs.rcell)
avect = [ 2.482, 0.000, 0.000]
bvect = [ 0.827, 2.340, 0.000]
cvect = [ 0.000, 0.000, 4.053]
origin = [ 0.000, 0.000, 0.000]
natoms = 2
natypes = 1
symbols = ('Fe',)
pbc = [ True True True]
per-atom properties = ['atype', 'pos']
id atype pos[0] pos[1] pos[2]
0 1 0.000 0.000 0.000
1 1 1.655 1.170 2.027
4.1.3. Build the free surface configuration
The resulting system is then more friendly to DFT by having fewer atoms in the periodic directions. Here, we also use the vacuumwidth parameter to insert a region of vacuum as that is also computationally favorable for DFT calculations.
[15]:
print(fs.surface(shift=fs.shifts[0], minwidth=10, vacuumwidth=15))
avect = [ 2.482, 0.000, 0.000]
bvect = [ 0.827, 2.340, 0.000]
cvect = [ 0.000, 0.000, 27.159]
origin = [ 0.000, 0.000, -7.500]
natoms = 6
natypes = 1
symbols = ('Fe',)
pbc = [ True True False]
per-atom properties = ['atype', 'pos']
id atype pos[0] pos[1] pos[2]
0 1 0.000 0.000 1.013
1 1 1.655 1.170 3.040
2 1 0.000 0.000 5.066
3 1 1.655 1.170 7.093
4 1 0.000 0.000 9.120
5 1 1.655 1.170 11.146
[16]:
print(fs.surfacearea)
5.808143987967931
4.3. hcp unit cell
Hexagonal systems can also be worked with using Miller-Bravais planes and vectors.
4.3.1. Define an hcp unit cell and initialize a FreeSurface
[17]:
a = uc.set_in_units(3.1853, 'angstrom')
c = uc.set_in_units(5.1853, 'angstrom')
box = am.Box.hexagonal(a, c)
atoms = am.Atoms(pos = [[0.0, 0.0, 0.0], [1/3, 2/3, 0.5]])
ucell = am.System(atoms=atoms, box=box, scale=True, symbols='Mg')
print(ucell)
avect = [ 3.185, 0.000, 0.000]
bvect = [-1.593, 2.759, 0.000]
cvect = [ 0.000, 0.000, 5.185]
origin = [ 0.000, 0.000, 0.000]
natoms = 2
natypes = 1
symbols = ('Mg',)
pbc = [ True True True]
per-atom properties = ['atype', 'pos']
id atype pos[0] pos[1] pos[2]
0 1 0.000 0.000 0.000
1 1 0.000 1.839 2.593
[18]:
plane_hkl = [1, 0, -1, 0]
fs = am.defect.FreeSurface(plane_hkl, ucell)
4.1.2. Check the identified rotated system
If the plane is defined using Miller-Bravais indices, then the uvws will be converted to them as well.
[19]:
print(fs.uvws)
[[-0.33333333 0.66666667 -0.33333333 0. ]
[ 0. 0. -0. 1. ]
[ 0.66666667 -0.33333333 -0.33333333 0. ]]
[20]:
print(fs.rcell)
avect = [ 3.185, 0.000, 0.000]
bvect = [ 0.000, 5.185, 0.000]
cvect = [-1.593, 0.000, 2.759]
origin = [ 0.000, 0.000, 0.000]
natoms = 2
natypes = 1
symbols = ('Mg',)
pbc = [ True True True]
per-atom properties = ['atype', 'pos']
id atype pos[0] pos[1] pos[2]
0 1 1.593 0.000 2.759
1 1 1.593 2.593 0.920
[21]:
print(fs.shifts)
[[0. 0. 0.91951691]
[0. 0. 2.29879227]]
4.1.3. Build the free surface configuration
If you give both sizemults and minwidth, then the larger of the two will be used for the determining the final system width.
[22]:
print(fs.surface(shift=fs.shifts[0], sizemults=(1, 1, 6), minwidth=10))
avect = [ 3.185, 0.000, 0.000]
bvect = [ 0.000, 5.185, 0.000]
cvect = [-9.556, 0.000, 16.551]
origin = [ 0.000, 0.000, 0.000]
natoms = 12
natypes = 1
symbols = ('Mg',)
pbc = [ True True False]
per-atom properties = ['atype', 'pos']
id atype pos[0] pos[1] pos[2]
0 1 -1.593 0.000 3.678
1 1 1.593 2.593 1.839
2 1 -3.185 0.000 6.437
3 1 0.000 2.593 4.598
4 1 -4.778 0.000 9.195
5 1 -1.593 2.593 7.356
6 1 -6.371 0.000 11.954
7 1 -3.185 2.593 10.115
8 1 -7.963 0.000 14.712
9 1 -4.778 2.593 12.873
10 1 0.000 0.000 0.920
11 1 -6.371 2.593 15.632
Also, note that the two shifts result in different termination planes. Both can be retrived from the same FreeSurface object by simply giving surface() a different shift value.
[23]:
print(fs.surface(shift=fs.shifts[1], sizemults=(1, 1, 6), minwidth=10))
avect = [ 3.185, 0.000, 0.000]
bvect = [ 0.000, 5.185, 0.000]
cvect = [-9.556, 0.000, 16.551]
origin = [ 0.000, 0.000, 0.000]
natoms = 12
natypes = 1
symbols = ('Mg',)
pbc = [ True True False]
per-atom properties = ['atype', 'pos']
id atype pos[0] pos[1] pos[2]
0 1 -1.593 0.000 5.057
1 1 -1.593 2.593 3.218
2 1 -3.185 0.000 7.816
3 1 -3.185 2.593 5.977
4 1 -4.778 0.000 10.574
5 1 -4.778 2.593 8.735
6 1 -6.371 0.000 13.333
7 1 -6.371 2.593 11.494
8 1 -7.963 0.000 16.092
9 1 -7.963 2.593 14.253
10 1 0.000 0.000 2.299
11 1 0.000 2.593 0.460