Introduction to atomman: Stacking fault configuration generator

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

Disclaimers

1. Introduction

This Notebook outlines the StackingFault class which generates atomic systems that can be used to evaluate stacking fault energies. See the 04.6 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.3.2
Notebook executed on 2020-04-14

2. Theory

A stacking fault atomic configuration can be constructed by

  1. Generating a bulk crystalline system with all periodic boundaries.

  2. Slicing the system along a geometrical plane creating two ideal planar-sliced free surfaces.

  3. Specifying a geometrical fault plane that is parallel to the free surfaces and divides the atoms (roughly) in half.

  4. 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}\)

\[\mathbf{d} = a_1 \mathbf{a_1} + a_2 \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:

  1. 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.

  2. 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. 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]
 [ 0.00000000e+00  0.00000000e+00  0.00000000e+00]
 [ 0.00000000e+00  0.00000000e+00  0.00000000e+00]
 [ 8.88178420e-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]
 [ 8.88178420e-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]
 [ 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 -2.22044605e-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 -2.22044605e-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 = [ 2.291  1.984  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.  -0.5  0.5]
 [ 0.5 -0.5  0. ]
 [-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   1.432   0.827   4.677
      1       1   2.864   1.653   2.338
      2       1   1.432   2.480   0.000
[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.  -0.5  0.5]
[ 0.5 -0.5  0. ]
The Cartesian shift vectors:
[ 2.86378246e+00  2.86439546e-16 -2.38358656e-16]
[ 1.43189123e+00  2.48010836e+00 -2.46533587e-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   0.827   5.846
      1       1   2.864   1.653   3.507
      2       1   1.432   2.480   1.169
      3       1   1.432   0.827  12.860
      4       1   2.864   1.653  10.522
      5       1   1.432   2.480   8.184
      6       1   1.432   0.827  19.875
      7       1   2.864   1.653  17.537
      8       1   1.432   2.480  15.199
      9       1   1.432   0.827  26.890
     10       1   2.864   1.653  24.552
     11       1   1.432   2.480  22.214

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 = [ 1.146  0.992  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 = [ 1.146 -0.992  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.