Introduction to atomman: Stacking fault generator
Lucas M. Hale, lucas.hale@nist.gov, Materials Science and Engineering Division, NIST.
1. Introduction
This Notebook outlines the StackingFault class which generates atomic systems that can be used to evaluate stacking fault energies. See the 4.5 Gamma surface plotting Notebook for analysis and plotting of computed stacking fault energies.
Completely reworked for version 1.3.2 the new StackingFault class is a child of FreeSurface, and therefore has built-in tools for easier generation.
Documentation for the older StackingFault class can be found on github by exploring past versions of this Notebook.
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 stacking fault atomic configuration 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.
Specifying a geometrical fault plane that is parallel to the free surfaces and divides the atoms (roughly) in half.
Rigidly shifting all atoms on one side of the fault plane by the same shift vector, \(\mathbf{d}\).
The first two steps correspond with the construction of a planar-sliced free surface configuration. The corresponding theory for those steps is outlined in the Free surface configuration generator Notebook.
Once the free surface has been constructed, a stacking fault configuration can be defined by specifying the location to position the fault plane, and a fault shift \(\mathbf{d}\) to apply. For the plane position, it is usually convenient to select the halfway point in the system. If the system was constructed with an even number of replicas in the direction normal to the fault/surface planes, then placing the fault plane in the middle of the system will result in the fault plane matching the free surface plane.
As the fault plane can be explored in 2 dimensions, the applied shift, \(\mathbf{d}\), can be expressed using two non-parallel vectors within the fault plane, \(\mathbf{a_1}\) and \(\mathbf{a_2}\). Taking \(\mathbf{a_1}\) and \(\mathbf{a_2}\) as lattice slip vectors (i.e. vectors between two symmetrically identical lattice sites) will define a 2D cell encompassing all possible unique stacking fault shifts for the fault plane. \(\mathbf{d}\) can then be represented using fractional coordinates, \(a_1\) and \(a_2\), relative to \(\mathbf{a_1}\) and \(\mathbf{a_2}\)
If the free surface configuration was generated using a primitive unit cell, then the two in-plane box vectors of the intermediate rotated cell will correspond to the shortest in-plane lattice vectors. As such, using those box vectors will result in the smallest periodic 2D cell necessary to fully explore the stacking fault map.
3. Class basics
The class is set up to generate stacking fault configuration systems in three steps:
Initialize the class and specify the (hkl) fault/surface plane 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.
Construct stacking fault configurations by rigidly shifting all atoms on one side of the fault plane.
3.1. StackingFault initialization
Class initializer. Identifies the proper rotations for the given hkl plane and cutboxvector, and creates the rotated cell.
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.
a1vect_uvw (array-like object, optional) The crystal vector to use for one of the two shifting vectors. If not given, will be set to the shortest in-plane box vector of rcell.
a2vect_uvw (array-like object, optional) The crystal vector to use for one of the two shifting vectors. If not given, will be set to the other in-plane box vector of rcell.
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. StackingFault.surface()
Generates the free surface atomic system, which is used as the basis for generating the stacking fault configuration(s).
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.
faultpos_rel (float, optional) The position to place the slip plane within the system given as a relative coordinate along the out-of-plane direction. faultpos_rel and faultpos_cart cannot both be given. Default value is 0.5 if faultpos_cart is also not given.
faultpos_cart (float, optional) The position to place the slip plane within the system given as a Cartesian coordinate along the out-of-plane direction. faultpos_rel and faultpos_cart cannot both be given.
3.3. StackingFault.fault()
Generates a fault configuration by displacing all atoms above the slip plane.
Parameters
a1 (float, optional) The fractional coordinate of a1vect to shift by. Default value is 0.0.
a2 (float, optional) The fractional coordinate of a2vect to shift by. Default value is 0.0.
outofplane (float, optional) An out-of-plane shift, given in absolute units. Default value is 0.0.
faultshift (array-like object, optional) The full shifting vector to displace the atoms above the slip plane by. Cannot be given with a1, a2, or outofplane.
minimum_r (float, optional) Added version 1.2.7 Specifies the minimum allowed interatomic spacing across the slip plane. If any sets of atoms are closer than this value then the outofplane shift is increased. Default value is None, which performs no adjustment.
Returns
(atomman.System) The atomic configuration with stacking fault shift
3.4. StackingFault.iterfaultmap()
Iterates over generalized stacking fault configurations associated with a 2D map of equally spaced a1, a2 coordinates.
Parameters
num_a1 (int) The number of a1 values to generate systems for. Default value is 1 (only generate for a1=0.0).
num_a2 (int) The number of a2 values to generate systems for. Default value is 1 (only generate for a2=0.0).
outofplane (float, optional) An out-of-plane shift, given in absolute units. Default value is 0.0.
minimum_r (float, optional) Added version 1.2.7 Specifies the minimum allowed interatomic spacing across the slip plane. If any sets of atoms are closer than this value then the outofplane shift is increased. Default value is None, which performs no adjustment.
Yields
a1 (float) The a1 fractional coordinate of a1vect.
a2 (float) The a2 fractional coordinate of a2vect.
(atomman.System) The fault configuration associated with the a1, a2 shift.
4. Examples
4.1. fcc \(\{111\}\)
4.1.1. Define an fcc conventional cell and initialize a StackingFault
[2]:
a = uc.set_in_units(4.05, 'angstrom')
atoms = am.Atoms(pos=[[0.0, 0.0, 0.0], [0.5, 0.5, 0.0], [0.5, 0.0, 0.5], [0.0, 0.5, 0.5]])
ucell = am.System(atoms=atoms, box=am.Box.cubic(a), scale=True, symbols='Al')
[3]:
plane_hkl = [ 1, 1, 1]
sf = am.defect.StackingFault(plane_hkl, ucell)
4.1.2. Check the identified rotated system
The StackingFault class is a child of FreeSurface, so it has all the same attributes that you can check. In addition, default values of a1vect and a2vect are set during the class initialization.
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.
a1vect_uvw, a2vect_uvw are the two lattice fault vectors as crystallographic Miller vectors relative to ucell.
a1vect_cart, a2vect_cart are the two lattice fault vectors as Cartesian vectors relative to rcell.
[4]:
print(sf.uvws)
[[-1. 1. 0.]
[-1. 0. 1.]
[ 1. 1. 1.]]
[5]:
print(sf.rcell)
avect = [ 5.728, 0.000, 0.000]
bvect = [ 2.864, 4.960, 0.000]
cvect = [ 0.000, 0.000, 7.015]
origin = [ 0.000, 0.000, 0.000]
natoms = 12
natypes = 1
symbols = ('Al',)
pbc = [ True True True]
per-atom properties = ['atype', 'pos']
id atype pos[0] pos[1] pos[2]
0 1 5.728 4.960 0.000
1 1 7.159 2.480 0.000
2 1 4.296 2.480 0.000
3 1 8.591 4.960 0.000
4 1 1.432 0.827 2.338
5 1 4.296 0.827 2.338
6 1 2.864 3.307 2.338
7 1 5.728 3.307 2.338
8 1 5.728 1.653 4.677
9 1 4.296 4.134 4.677
10 1 7.159 4.134 4.677
11 1 2.864 1.653 4.677
The default a1vect_uvw and a2vect_uvw values are the two in-plane uvw rotation vectors.
[6]:
print('The crystal shift vectors:')
print(sf.a1vect_uvw)
print(sf.a2vect_uvw)
print('The Cartesian shift vectors:')
print(sf.a1vect_cart)
print(sf.a2vect_cart)
The crystal shift vectors:
[-1. 1. 0.]
[-1. 0. 1.]
The Cartesian shift vectors:
[ 5.72756493e+00 -5.57532102e-16 -8.29754091e-16]
[ 2.86378246e+00 4.96021673e+00 -8.29754091e-16]
As the StackingFault class was initialized with the conventional unit cell, the default shifting vectors are not the smallest possible lattice vectors. The shift values can be manually set either 1. during initialization, 2. by setting to the a1vect_uvw, a2vect_uvw attributes after initialization, or 3. when calling one of the fault configuration generation methods.
Note: the choice below is the most optimum for fcc {111} planes as the two shift vectors are the smallest non-parallel in-plane lattice vectors. The typical fcc {111} representation seen in papers with shifts of \(\frac{1}{2}[\bar{1}10]\) and \(\frac{1}{2}[11\bar{2}]\) vectors is actually a double cell and therefore requires sampling twice the area. The GammaSurface class accounts for periodicity, meaning that the standard representation can be plotted based on measured results from the optimum grid of shifts.
[7]:
sf.a1vect_uvw = [-0.5, 0.5, 0.0]
sf.a2vect_uvw = [-0.5, 0.0, 0.5]
print('The crystal shift vectors:')
print(sf.a1vect_uvw)
print(sf.a2vect_uvw)
print('The Cartesian shift vectors:')
print(sf.a1vect_cart)
print(sf.a2vect_cart)
The crystal shift vectors:
[-0.5 0.5 0. ]
[-0.5 0. 0.5]
The Cartesian shift vectors:
[ 2.86378246e+00 -2.78766051e-16 -4.14877046e-16]
[ 1.43189123e+00 2.48010836e+00 -4.14877046e-16]
4.1.3. Build the free surface configuration
The stacking fault configurations are modifications of the free surface configuration, so the next step is to build the free surface system. For stacking fault configurations, two different parallel geometrical planes need to be positioned inside the system: the free surface cut plane and the fault plane. Their positioning with respect to the atomic system is controlled by
shift rigidly shifts the atoms. This is best used to position the cut plane
faultpos_rel, faultpos_cart position the fault plane in the shifted system either at a relative coordinate or an absolute Cartesian coordinate, respectively.
[8]:
print(sf.shifts)
[[0. 0. 1.1691343 ]
[0. 0. 3.50740289]
[0. 0. 5.84567148]]
When building the free surface configuration, the default fault position corresponds to faultpos_rel = 0.5, i.e. it is placed in the center of the system. If the number of replicas in the cut box vector direction are even, then this default fault position will result in the cut and fault planes bisecting the rotated cell along atomically equivalent planes.
[9]:
surfacesystem = sf.surface(shift=sf.shifts[0], minwidth=15, even=True)
print(surfacesystem)
avect = [ 5.728, 0.000, 0.000]
bvect = [ 2.864, 4.960, 0.000]
cvect = [ 0.000, 0.000, 28.059]
origin = [ 0.000, 0.000, 0.000]
natoms = 48
natypes = 1
symbols = ('Al',)
pbc = [ True True False]
per-atom properties = ['atype', 'pos']
id atype pos[0] pos[1] pos[2]
0 1 5.728 4.960 1.169
1 1 7.159 2.480 1.169
2 1 4.296 2.480 1.169
3 1 8.591 4.960 1.169
4 1 1.432 0.827 3.507
5 1 4.296 0.827 3.507
6 1 2.864 3.307 3.507
7 1 5.728 3.307 3.507
8 1 5.728 1.653 5.846
9 1 4.296 4.134 5.846
10 1 7.159 4.134 5.846
11 1 2.864 1.653 5.846
12 1 5.728 4.960 8.184
13 1 7.159 2.480 8.184
14 1 4.296 2.480 8.184
15 1 8.591 4.960 8.184
16 1 1.432 0.827 10.522
17 1 4.296 0.827 10.522
18 1 2.864 3.307 10.522
19 1 5.728 3.307 10.522
20 1 5.728 1.653 12.860
21 1 4.296 4.134 12.860
22 1 7.159 4.134 12.860
23 1 2.864 1.653 12.860
24 1 5.728 4.960 15.199
25 1 7.159 2.480 15.199
26 1 4.296 2.480 15.199
27 1 8.591 4.960 15.199
28 1 1.432 0.827 17.537
29 1 4.296 0.827 17.537
30 1 2.864 3.307 17.537
31 1 5.728 3.307 17.537
32 1 5.728 1.653 19.875
33 1 4.296 4.134 19.875
34 1 7.159 4.134 19.875
35 1 2.864 1.653 19.875
36 1 5.728 4.960 22.214
37 1 7.159 2.480 22.214
38 1 4.296 2.480 22.214
39 1 8.591 4.960 22.214
40 1 1.432 0.827 24.552
41 1 4.296 0.827 24.552
42 1 2.864 3.307 24.552
43 1 5.728 3.307 24.552
44 1 5.728 1.653 26.890
45 1 4.296 4.134 26.890
46 1 7.159 4.134 26.890
47 1 2.864 1.653 26.890
4.1.4. Check the free surface system and fault position
Once the free surface configuration has been created, the quality of the fault position can be checked with the following class attributes:
abovefault is a boolean mask indicating which atoms are above the fault plane. A good fault position should have (roughly) the same number of atoms above and below the plane.
faultpos_rel is the relative position of the fault plane.
faultpos_cart is the absolute Cartesian position of the fault plane.
[10]:
print(f'{sf.system.natoms} total atoms')
print(f'{sf.abovefault.sum()} atoms above the fault')
print()
print(f'Relative fault plane position: {sf.faultpos_rel}')
print(f'Absolute fault plane position: {sf.faultpos_cart}')
48 total atoms
24 atoms above the fault
Relative fault plane position: 0.5
Absolute fault plane position: 14.029611541307906
A different fault position than the default faultpos_rel=0.5 can be set either 1. when surface() is called, 2. by setting to either the faultpos_rel or faultpos_cart after calling surface(), or 3. when calling one of the fault configuration generation methods.
[11]:
sf.faultpos_rel = 0.57
print(f'{sf.system.natoms} total atoms')
print(f'{sf.abovefault.sum()} atoms above the fault')
print()
print(f'Relative fault plane position: {sf.faultpos_rel}')
print(f'Absolute fault plane position: {sf.faultpos_cart}')
48 total atoms
20 atoms above the fault
Relative fault plane position: 0.57
Absolute fault plane position: 15.99375715709101
[12]:
sf.faultpos_cart = 14.029611541307906
print(f'{sf.system.natoms} total atoms')
print(f'{sf.abovefault.sum()} atoms above the fault')
print()
print(f'Relative fault plane position: {sf.faultpos_rel}')
print(f'Absolute fault plane position: {sf.faultpos_cart}')
48 total atoms
24 atoms above the fault
Relative fault plane position: 0.5
Absolute fault plane position: 14.029611541307906
4.1.5. Generate stacking fault configuration(s)
A single stacking fault configuration for a specific shift can be generated using the fault() method. A full 1D or 2D array of stacking fault configurations can be generated using the iterfaultmap() method.
[13]:
faultsystem = sf.fault(a1=0.5, a2=0.0)
print('atoms in faultsystem are displaced by:')
print(am.displacement(surfacesystem, faultsystem))
atoms in faultsystem are displaced by:
[[ 0.00000000e+00 0.00000000e+00 0.00000000e+00]
[ 0.00000000e+00 0.00000000e+00 0.00000000e+00]
[ 0.00000000e+00 0.00000000e+00 0.00000000e+00]
[ 0.00000000e+00 0.00000000e+00 0.00000000e+00]
[ 0.00000000e+00 0.00000000e+00 0.00000000e+00]
[ 0.00000000e+00 0.00000000e+00 0.00000000e+00]
[ 4.44089210e-16 0.00000000e+00 0.00000000e+00]
[ 0.00000000e+00 0.00000000e+00 0.00000000e+00]
[ 0.00000000e+00 0.00000000e+00 0.00000000e+00]
[ 0.00000000e+00 0.00000000e+00 0.00000000e+00]
[ 0.00000000e+00 0.00000000e+00 0.00000000e+00]
[ 0.00000000e+00 0.00000000e+00 0.00000000e+00]
[ 0.00000000e+00 0.00000000e+00 0.00000000e+00]
[ 0.00000000e+00 0.00000000e+00 0.00000000e+00]
[ 0.00000000e+00 0.00000000e+00 0.00000000e+00]
[ 0.00000000e+00 0.00000000e+00 0.00000000e+00]
[ 0.00000000e+00 0.00000000e+00 0.00000000e+00]
[ 0.00000000e+00 0.00000000e+00 0.00000000e+00]
[ 4.44089210e-16 0.00000000e+00 0.00000000e+00]
[ 0.00000000e+00 0.00000000e+00 0.00000000e+00]
[ 0.00000000e+00 0.00000000e+00 0.00000000e+00]
[ 0.00000000e+00 0.00000000e+00 0.00000000e+00]
[ 0.00000000e+00 0.00000000e+00 0.00000000e+00]
[ 0.00000000e+00 0.00000000e+00 0.00000000e+00]
[ 1.43189123e+00 0.00000000e+00 0.00000000e+00]
[ 1.43189123e+00 0.00000000e+00 0.00000000e+00]
[ 1.43189123e+00 0.00000000e+00 0.00000000e+00]
[ 1.43189123e+00 0.00000000e+00 0.00000000e+00]
[ 1.43189123e+00 -1.11022302e-16 0.00000000e+00]
[ 1.43189123e+00 -1.11022302e-16 0.00000000e+00]
[ 1.43189123e+00 0.00000000e+00 0.00000000e+00]
[ 1.43189123e+00 0.00000000e+00 0.00000000e+00]
[ 1.43189123e+00 -2.22044605e-16 0.00000000e+00]
[ 1.43189123e+00 0.00000000e+00 0.00000000e+00]
[ 1.43189123e+00 0.00000000e+00 0.00000000e+00]
[ 1.43189123e+00 -2.22044605e-16 0.00000000e+00]
[ 1.43189123e+00 0.00000000e+00 0.00000000e+00]
[ 1.43189123e+00 0.00000000e+00 0.00000000e+00]
[ 1.43189123e+00 0.00000000e+00 0.00000000e+00]
[ 1.43189123e+00 0.00000000e+00 0.00000000e+00]
[ 1.43189123e+00 -1.11022302e-16 0.00000000e+00]
[ 1.43189123e+00 -1.11022302e-16 0.00000000e+00]
[ 1.43189123e+00 0.00000000e+00 0.00000000e+00]
[ 1.43189123e+00 0.00000000e+00 0.00000000e+00]
[ 1.43189123e+00 -2.22044605e-16 0.00000000e+00]
[ 1.43189123e+00 0.00000000e+00 0.00000000e+00]
[ 1.43189123e+00 0.00000000e+00 0.00000000e+00]
[ 1.43189123e+00 -2.22044605e-16 0.00000000e+00]]
Generate configurations along the slip line with StackingFault.iterfaultmap(). Note that the displacements may wrap due to periodic boundaries.
[14]:
for a1, a2, faultsystem in sf.iterfaultmap(num_a1=5, num_a2=5):
d = am.displacement(surfacesystem, faultsystem)[-1]
print(f'a1 = {a1:.1f}, a2 = {a2:.1f}, d = [{d[0]: 5.3f} {d[1]: 5.3f} {d[2]: 5.3f}]')
a1 = 0.0, a2 = 0.0, d = [ 0.000 0.000 0.000]
a1 = 0.2, a2 = 0.0, d = [ 0.573 0.000 0.000]
a1 = 0.4, a2 = 0.0, d = [ 1.146 -0.000 0.000]
a1 = 0.6, a2 = 0.0, d = [ 1.718 -0.000 0.000]
a1 = 0.8, a2 = 0.0, d = [ 2.291 -0.000 0.000]
a1 = 0.0, a2 = 0.2, d = [ 0.286 0.496 0.000]
a1 = 0.2, a2 = 0.2, d = [ 0.859 0.496 0.000]
a1 = 0.4, a2 = 0.2, d = [ 1.432 0.496 0.000]
a1 = 0.6, a2 = 0.2, d = [ 2.005 0.496 0.000]
a1 = 0.8, a2 = 0.2, d = [ 2.577 0.496 0.000]
a1 = 0.0, a2 = 0.4, d = [ 0.573 0.992 0.000]
a1 = 0.2, a2 = 0.4, d = [ 1.146 0.992 0.000]
a1 = 0.4, a2 = 0.4, d = [ 1.718 0.992 0.000]
a1 = 0.6, a2 = 0.4, d = [ 2.291 0.992 0.000]
a1 = 0.8, a2 = 0.4, d = [-2.864 0.992 0.000]
a1 = 0.0, a2 = 0.6, d = [ 0.859 1.488 0.000]
a1 = 0.2, a2 = 0.6, d = [ 1.432 1.488 0.000]
a1 = 0.4, a2 = 0.6, d = [ 2.005 1.488 0.000]
a1 = 0.6, a2 = 0.6, d = [ 2.577 1.488 0.000]
a1 = 0.8, a2 = 0.6, d = [-2.577 1.488 0.000]
a1 = 0.0, a2 = 0.8, d = [ 1.146 1.984 0.000]
a1 = 0.2, a2 = 0.8, d = [ 1.718 1.984 0.000]
a1 = 0.4, a2 = 0.8, d = [-0.573 -2.976 0.000]
a1 = 0.6, a2 = 0.8, d = [ 0.000 -2.976 0.000]
a1 = 0.8, a2 = 0.8, d = [-2.291 1.984 0.000]
4.2. fcc {111} using primitive unit cell
Just like with the FreeSurface class, a primitive unit cell can be given with a conventional (hkl) slip plane by specifying what the conventional cell’s setting/centering is.
The advantages of doing this for StackingFaults are
the resulting systems will be smaller in the in-plane directions making DFT calculations easier, and
the default shifting vectors will be the smallest in-plane lattice vectors.
4.2.1. Define an fcc primitive cell and initialize a StackingFault
[15]:
a = uc.set_in_units(4.05, 'angstrom')
atoms = am.Atoms(pos=[[0.0, 0.0, 0.0]])
ucell = am.System(atoms=atoms, box=am.Box.trigonal(a/2**0.5, alpha=60), scale=True, symbols='Al')
[16]:
plane_hkl = [ 1, 1, 1]
conventional_setting = 'f'
sf = am.defect.StackingFault(plane_hkl, ucell, conventional_setting=conventional_setting)
4.2.2. Check the identified rotated system and shift vectors
[17]:
print(sf.uvws)
[[-0.5 0. 0.5]
[ 0. -0.5 0.5]
[ 1. 1. 1. ]]
[18]:
print(sf.rcell)
avect = [ 2.864, 0.000, 0.000]
bvect = [ 1.432, 2.480, 0.000]
cvect = [ 0.000, 0.000, 7.015]
origin = [ 0.000, 0.000, 0.000]
natoms = 3
natypes = 1
symbols = ('Al',)
pbc = [ True True True]
per-atom properties = ['atype', 'pos']
id atype pos[0] pos[1] pos[2]
0 1 4.296 2.480 0.000
1 1 2.864 1.653 2.338
2 1 1.432 0.827 4.677
[19]:
print('The crystal shift vectors:')
print(sf.a1vect_uvw)
print(sf.a2vect_uvw)
print('The Cartesian shift vectors:')
print(sf.a1vect_cart)
print(sf.a2vect_cart)
The crystal shift vectors:
[-0.5 0. 0.5]
[ 0. -0.5 0.5]
The Cartesian shift vectors:
[ 2.86378246e+00 -1.10922045e-17 9.57795207e-16]
[1.43189123e+00 2.48010836e+00 1.41222587e-16]
Now, we have an rcell with 1/4 the atoms and the algorithm automatically identified the same “best” shifts that we had to manually specify for the conventional cell!
4.2.3. Build the free surface configuration
The vacuumwidth parameter of surface() inserts a volume of vacuum at the free surface of this width. DFT calculations are more efficient for free surfaces introduced with a volume of vacuum rather than making one dimension non-periodic.
[20]:
surfacesystem = sf.surface(shift=sf.shifts[0], minwidth=15, even=True, vacuumwidth=20)
print(surfacesystem)
avect = [ 2.864, 0.000, 0.000]
bvect = [ 1.432, 2.480, 0.000]
cvect = [ 0.000, 0.000, 48.059]
origin = [ 0.000, 0.000, -10.000]
natoms = 12
natypes = 1
symbols = ('Al',)
pbc = [ True True False]
per-atom properties = ['atype', 'pos']
id atype pos[0] pos[1] pos[2]
0 1 1.432 2.480 1.169
1 1 2.864 1.653 3.507
2 1 1.432 0.827 5.846
3 1 1.432 2.480 8.184
4 1 2.864 1.653 10.522
5 1 1.432 0.827 12.860
6 1 1.432 2.480 15.199
7 1 2.864 1.653 17.537
8 1 1.432 0.827 19.875
9 1 1.432 2.480 22.214
10 1 2.864 1.653 24.552
11 1 1.432 0.827 26.890
4.2.4. Check the free surface system and fault position
The vacuum region changes the box dimensions but not the Cartesian atomic positions. As such, the Cartesian fault positions are the same with and without the vacuum region, but the relative fault positions are different for all but one value. Half of the vacuum region is added to the bottom of the system (by shifting origin) and the other half is added to the top. This means that a relative shift of 0.5 results in the same Cartesian location.
[21]:
print(f'{sf.system.natoms} total atoms')
print(f'{sf.abovefault.sum()} atoms above the fault')
print()
print(f'Relative fault plane position: {sf.faultpos_rel}')
print(f'Absolute fault plane position: {sf.faultpos_cart}')
12 total atoms
6 atoms above the fault
Relative fault plane position: 0.5
Absolute fault plane position: 14.029611541307904
[22]:
sf.faultpos_rel = 0.57
print(f'{sf.system.natoms} total atoms')
print(f'{sf.abovefault.sum()} atoms above the fault')
print()
print(f'Relative fault plane position: {sf.faultpos_rel}')
print(f'Absolute fault plane position: {sf.faultpos_cart}')
12 total atoms
5 atoms above the fault
Relative fault plane position: 0.57
Absolute fault plane position: 17.39375715709101
[23]:
sf.faultpos_cart = 14.029611541307906
print(f'{sf.system.natoms} total atoms')
print(f'{sf.abovefault.sum()} atoms above the fault')
print()
print(f'Relative fault plane position: {sf.faultpos_rel}')
print(f'Absolute fault plane position: {sf.faultpos_cart}')
12 total atoms
6 atoms above the fault
Relative fault plane position: 0.5000000000000001
Absolute fault plane position: 14.029611541307906
4.2.5. Generate stacking fault configuration(s)
Further support for DFT is that the fault() and iterfaultmap() methods also have a minimum_r parameter that can be used to ensure all atoms across the fault plane remain at least a set distance apart. This can help with relaxations when certain faults may otherwise position atoms very close together.
Here, minimum_r is set to be 0.85 of the unit cell’s smallest interatomic vector r0. This results in some shifts in the z direction for certain fault configurations.
[24]:
for a1, a2, faultsystem in sf.iterfaultmap(num_a1=5, num_a2=5, minimum_r=0.85*ucell.r0()):
d = am.displacement(surfacesystem, faultsystem)[-1]
print(f'a1 = {a1:.1f}, a2 = {a2:.1f}, d = [{d[0]: 5.3f} {d[1]: 5.3f} {d[2]: 5.3f}]')
a1 = 0.0, a2 = 0.0, d = [ 0.000 0.000 0.000]
a1 = 0.2, a2 = 0.0, d = [ 0.573 0.000 0.000]
a1 = 0.4, a2 = 0.0, d = [ 1.146 0.000 0.000]
a1 = 0.6, a2 = 0.0, d = [-1.146 0.000 0.000]
a1 = 0.8, a2 = 0.0, d = [-0.573 0.000 0.000]
a1 = 0.0, a2 = 0.2, d = [ 0.286 0.496 0.000]
a1 = 0.2, a2 = 0.2, d = [ 0.859 0.496 0.004]
a1 = 0.4, a2 = 0.2, d = [-1.432 0.496 0.073]
a1 = 0.6, a2 = 0.2, d = [-0.859 0.496 0.004]
a1 = 0.8, a2 = 0.2, d = [-0.286 0.496 0.000]
a1 = 0.0, a2 = 0.4, d = [ 0.573 0.992 0.000]
a1 = 0.2, a2 = 0.4, d = [-0.286 -1.488 0.073]
a1 = 0.4, a2 = 0.4, d = [-1.146 0.992 0.073]
a1 = 0.6, a2 = 0.4, d = [-0.573 0.992 0.000]
a1 = 0.8, a2 = 0.4, d = [ 0.000 0.992 0.000]
a1 = 0.0, a2 = 0.6, d = [-0.573 -0.992 0.000]
a1 = 0.2, a2 = 0.6, d = [ 0.000 -0.992 0.004]
a1 = 0.4, a2 = 0.6, d = [ 0.573 -0.992 0.000]
a1 = 0.6, a2 = 0.6, d = [-0.286 1.488 0.000]
a1 = 0.8, a2 = 0.6, d = [-1.146 -0.992 0.000]
a1 = 0.0, a2 = 0.8, d = [-0.286 -0.496 0.000]
a1 = 0.2, a2 = 0.8, d = [ 0.286 -0.496 0.000]
a1 = 0.4, a2 = 0.8, d = [ 0.859 -0.496 0.000]
a1 = 0.6, a2 = 0.8, d = [ 1.432 -0.496 0.000]
a1 = 0.8, a2 = 0.8, d = [-0.859 -0.496 0.000]
4.3. hcp {0001} slip plane
The system rotation vectors and the shifting vectors can be given in Miller-Bravais indices if the reference unit cell is hexagonal.
4.3.1. Define an fcc primitive cell and initialize a StackingFault
[25]:
a = uc.set_in_units(3.18, 'angstrom')
c = uc.set_in_units(5.17, 'angstrom')
atoms = am.Atoms(pos=[[0.0, 0.0, 0.0], [1/3, 2/3, 0.5]])
ucell = am.System(atoms=atoms, box=am.Box.hexagonal(a, c), scale=True)
[26]:
plane_hkl = [ 0, 0, 0, 1]
sf = am.defect.StackingFault(plane_hkl, ucell)
4.3.2. Check the identified rotated system and shift vectors
If plane_hkl was given in the 4-index representation, then uvws, a1vect_uvw and a2vect_uvw will be saved in the 4-index representation as well.
[27]:
print(sf.uvws)
[[ 0.66666667 -0.33333333 -0.33333333 0. ]
[ 0.33333333 0.33333333 -0.66666667 0. ]
[ 0. 0. -0. 1. ]]
[28]:
print(sf.rcell)
avect = [ 3.180, 0.000, 0.000]
bvect = [ 1.590, 2.754, 0.000]
cvect = [ 0.000, 0.000, 5.170]
origin = [ 0.000, 0.000, 0.000]
natoms = 2
natypes = 1
symbols = (None,)
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 3.180 1.836 2.585
[29]:
print('The crystal shift vectors:')
print(sf.a1vect_uvw)
print(sf.a2vect_uvw)
print('The Cartesian shift vectors:')
print(sf.a1vect_cart)
print(sf.a2vect_cart)
The crystal shift vectors:
[ 0.66666667 -0.33333333 -0.33333333 0. ]
[ 0.33333333 0.33333333 -0.66666667 0. ]
The Cartesian shift vectors:
[ 3.18000000e+00 -2.73153518e-16 0.00000000e+00]
[1.59 2.75396078 0. ]
The remaining steps are similar to those mentioned for the other systems above.