Introduction to atomman: Dislocation analysis tools
Lucas M. Hale, lucas.hale@nist.gov, Materials Science and Engineering Division, NIST.
1. Introduction
This Notebook describes the analysis tools contained in atomman.defect for finding and characterizing dislocations.
Library Imports
[1]:
# Standard Python libraries
from copy import deepcopy
import datetime
# http://www.numpy.org/
import numpy as np
# https://matplotlib.org/
import matplotlib.pyplot as plt
%matplotlib inline
# 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
1.1. Load atomic systems
A system containing a single Al \(\frac{a}{2}\)<111> edge dislocation was created by
Using the Dislocation class to solve the elasticity solution for the dislocation.
Generating a dislocation monopole and base reference system using Dislocation’s monopole method.
Performing a short low temerature anneal and energy minimization to relax the atomic configuration. During relaxation, the system is kept periodic along the dislocation line, but is non-periodic in the other two directions. Atoms at the edge of the non-periodic directions are held fixed.
All atomic dislocation analysis tools identify and characterize the presence of defects in an atomic system by comparing the local atomic environment around every atom to a reference state. Some of the tools in atomman do this by directly mapping each atom in the defect system to a corresponding defect-free base system. Using the Dislocation class allows for the appropriate base system to be constructed along with the unrelaxed dislocation configuration.
[2]:
base_system = am.load('atom_dump', 'files/fcc_Al_base.dump')
disl_system = am.load('atom_dump', 'files/fcc_Al_disl.dump')
The dislocation system is oriented with the line direction parallel to the x-axis and the slip plane perpendicular to the z-axis. As such, the Burgers vector of the full edge dislocation is entirely in the y direction.
[3]:
alat = uc.set_in_units(4.05, 'Å')
burgers = np.array([0.0, alat / 2**0.5, 0.0])
burgers
[3]:
array([0. , 2.86378246, 0. ])
[4]:
# plot original positions of atoms in system
fig = plt.figure(figsize=(6,6))
plt.plot(base_system.atoms.pos[:, 1], base_system.atoms.pos[:, 2], 'o')
plt.xlim([-25, 25])
plt.ylim([-25, 25])
plt.show()
[5]:
# plot positions of atoms in dislocation system
fig = plt.figure(figsize=(6,6))
plt.plot(disl_system.atoms.pos[:, 1], disl_system.atoms.pos[:, 2], 'o')
plt.xlim([-25, 25])
plt.ylim([-25, 25])
plt.show()
2. Slip vector
The slip vector is useful for characterizing the path that dislocations have moved between a current state and a reference state. It characterizes how much crystal slip happens between every atom and their nearest neighbors. The burgers vector of any moving dislocation can be estimated by looking at the magnitude and direction of the slip vector of any atoms that the dislocation passed by. The slip vector was originally outlined in Zimmerman, et. al, Phys Rev Lett 87, 165507.
where \(n_{\alpha}\) is the number of nearest neighbor atoms \(\beta\) of atom \(\alpha\), \(x_{i \alpha \beta}\) is the vector difference between atoms \(\alpha\) and \(\beta\) in the current configuration, and \(X_{i \alpha \beta}\) is the vector difference between atoms \(\alpha\) and \(\beta\) in the reference configuration.
Note 1: To be general with any crystal and dislocation, the slip vector computed here is not scaled by the number of slipped neighbors like it is in the original formulation. This means that the value of the slip vector returned will be a multiple of the Burgers vector that caused the slip. To estimate the Burgers vector of a dislocation that moved along a slip plane, the slip vector value needs to be divided by the total number of neighbor atoms that are across the slip plane.
Note 2: The reference system does not need to be of a perfect crystal and can simply be an earlier state. With this choice, the slip vector will reveal the slip deformation that occurred between the two states of the system. This can be useful for identifying how dislocations move in highly deformed/defective systems.
2.1. slip_vector()
Parameters
system_0 (atomman.system) The base/reference system to use.
system_1 (atomman.system) The defect/current system to use.
neighbors (atomman.NeighborList, optional) The neighbor list associated with system_0 to use. Either neighbors or cutoff must be given, or system_0 must have a neighbors attribute.
cutoff (float, optional) Cutoff distance for computing a neighbor list for system_0. Either neighbors or cutoff must be given, or system_0 have a neighbors attribute.
Returns
(numpy.ndarray) The computed slip vectors.
[6]:
# Make reference system's pbc match defect system's pbc
# This avoids free surfaces having large slip vectors
base_system.pbc = (True, False, False)
# Build the neighbor list using the reference system
neighbors0 = base_system.neighborlist(cutoff=0.9*alat)
# Compute the slip vector for all atoms
slip = am.defect.slip_vector(base_system, disl_system, neighbors=neighbors0)
As the reference system is the defect free system used to construct the dislocation system, the slip vector computed here shows the dislocation slip needed to insert the dislocation into the system. Note that the slip vector values on the two sides of the plane are equal but opposite.
[7]:
# plot slip vector components for all atoms
fig, ax = plt.subplots(ncols=3, figsize=(17,4.5))
im = ax[0].scatter(disl_system.atoms.pos[:, 1], disl_system.atoms.pos[:, 2], c=slip[:, 0], vmin=-10, vmax=10, cmap='jet')
plt.colorbar(im, ax=ax[0])
ax[0].set_xlim(-30, 30)
ax[0].set_ylim(-30, 30)
im = ax[1].scatter(disl_system.atoms.pos[:, 1], disl_system.atoms.pos[:, 2], c=slip[:, 1], vmin=-10, vmax=10, cmap='jet')
plt.colorbar(im, ax=ax[1])
ax[1].set_xlim(-30, 30)
ax[1].set_ylim(-30, 30)
im = ax[2].scatter(disl_system.atoms.pos[:, 1], disl_system.atoms.pos[:, 2], c=slip[:, 2], vmin=-10, vmax=10, cmap='jet')
plt.colorbar(im, ax=ax[2])
ax[2].set_xlim(-30, 30)
ax[2].set_ylim(-30, 30)
plt.show()
The Burgers vector for the dislocation can then be estimated from slip vector values for atoms along the slipped plane.
[8]:
# Find atoms with coordinates -30 < y < -25 and 0.0 < z < 3.0 and
matches = ((disl_system.atoms.pos[:, 1] > -30.0) & (disl_system.atoms.pos[:, 1] < -25.0)
&(disl_system.atoms.pos[:, 2] > 0.0) & (disl_system.atoms.pos[:, 2] < 3.0))
# Get average slip for the atoms
aveslip = slip[matches].mean(axis=0)
# Divide the slip by 3 to estimate the Burgers vector
# as fcc (111) planes have three cross-plane neighbors
burgers_est = aveslip / 3
print(f'actual Burgers vector: [{burgers[0]: 7.4f}, {burgers[1]: 7.4f}, {burgers[2]: 7.4f}]')
print(f'estimated Burgers vector: [{burgers_est[0]: 7.4f}, {burgers_est[1]: 7.4f}, {burgers_est[2]: 7.4f}]')
actual Burgers vector: [ 0.0000, 2.8638, 0.0000]
estimated Burgers vector: [-0.0016, 2.8662, 0.0060]
Try estimating the partial dislocation by looking at slip in the stacking fault area
[9]:
# Find atoms with coordinates -1 < y < 1 and 0.0 < z < 3.0 and
matches = ((disl_system.atoms.pos[:, 1] > -1.0) & (disl_system.atoms.pos[:, 1] < 1.0)
&(disl_system.atoms.pos[:, 2] > 0.0) & (disl_system.atoms.pos[:, 2] < 3.0))
# Get average slip for the atoms
aveslip = slip[matches].mean(axis=0)
# Divide the slip by 3 to estimate the Burgers vector
# as fcc (111) planes have three cross-plane neighbors
burgers_est = aveslip / 3
# Define the Shockley Burgers vector relative to the system definintion
shockley = alat * np.array([ 6**0.5 / 12, 2**0.5 / 4, 0.0])
print(f'actual Burgers vector: [{shockley[0]: 7.4f}, {shockley[1]: 7.4f}, {shockley[2]: 7.4f}]')
print(f'estimated Burgers vector: [{burgers_est[0]: 7.4f}, {burgers_est[1]: 7.4f}, {burgers_est[2]: 7.4f}]')
actual Burgers vector: [ 0.8267, 1.4319, 0.0000]
estimated Burgers vector: [ 0.7035, 1.4935, -0.0257]
3. Disregistry
The disregistry, \(\delta_i\) characterizes the planar spreading of a dislocation in direction \(\xi\) along a slip plane by measuring the difference in the displacements between the atomic planes of atoms just above, \(u_i^+\), and below, \(u_i^-\), the mathematical slip plane. The displacements are taken relative to a perfect crystal base configuration.
The algorithm used by atomman does the following:
Uses the base configuration to identify all atoms in the two atomic planes neighboring the slip plane.
For all atoms identified in #1, finds the unique \(\xi\) coordinates from the base configuration and computes \(u_i^+\) and \(u_i^-\) for the defect configuration relative to the base configuration.
Any \(u_i^+\) or \(u_i^-\) corresponding to the same \(\xi\) coordinates are averaged, and the values of \(u_i^+\) and \(u_i^-\) are linearly interpolated to all unique \(\xi\) coordinates. This is done as the \(\xi\) coordinates may differ for the two atomic planes.
The disregistry is computed for all unique \(\xi\) coordinates by finding the difference in the interpolated \(u_i^+\) and \(u_i^-\) displacements.
3.1. disregistry()
Parameters
basesystem (atomman.System) A perfect reference system with atoms directly corresponding to atoms in dislsystem.
dislsystem (atomman.System) A dislocation-containing system.
m (array-like object, optional) Unit vector defining the direction associated with the planar spreading along the slip plane. Default value is [1, 0, 0] (Cartesian x-coordinates).
n (array-like object, optional) Unit vector defining the normal to the slip plane. Must be perpendicular to xvector. Default value is [0, 1, 0] (Cartesian y-axis).
planepos (array-like object, optional) A position on the slip plane so that the plane can be fully defined. The slip plane position should fall between two planes of atoms. Default value is [0,0,0].
Returns
coord (numpy.ndarray) The (N,) array of unique coordinates (atomic columns) neighboring the slip plane.
disregistry (numpy.ndarray) A (N, 3) array of the dislocation’s disregistry at each x.
[10]:
xi, disreg = am.defect.disregistry(base_system, disl_system, m=[0,1,0], n=[0,0,1])
fig = plt.figure(figsize=(8,6))
plt.plot(xi, disreg[:, 0], label='x-disregistry', linewidth=3)
plt.plot(xi, disreg[:, 1], label='y-disregistry', linewidth=3)
plt.plot(xi, disreg[:, 2], label='z-disregistry', linewidth=3)
plt.xlabel('$\\xi$-coordinate ($\AA$)', size='x-large')
plt.ylabel('disregistry ($\AA$)', size='x-large')
plt.legend(fontsize='xx-large')
plt.legend()
plt.show()
4. Nye tensor
Update version 1.3.7: The Strain class has been added that provides a faster, more user-friendly interface for computing the same properties as the nye_tensor function. It is recommended to use it instead - See the 4.8. Strain class Juptyer Notebook for a description of the class and examples.
5. Differential displacement maps
Dislocation core structures can also be characterized using differential displacement maps as first used by Vitek, Perrin and Bowen.
For every pair of neighboring atoms \(\alpha\) and \(\beta\), a differential displacement vector, \(d_{i \alpha \beta}\) is computed as
where \(x_{i \alpha \beta}\) is the vector difference between atoms \(\alpha\) and \(\beta\) in the current (defect) configuration, and \(X_{i \alpha \beta}\) is the vector difference between atoms \(\alpha\) and \(\beta\) in the reference (perfect crystal) configuration.
Plots are then constructed where atomic positions are shown as circles with arrows representing components of \(d_{i \alpha \beta}\). The arrows are centered halfway between the two atoms \(\alpha\) and \(\beta\) for which that particular \(d_{i \alpha \beta}\) was computed. For most variations, the direction of the arrows indicate the two atoms \(\alpha\) and \(\beta\) that \(d_{i \alpha \beta}\) was computed for. Various plotting options and variations are described in more detail below.
Update version 1.3.2 A DifferentialDisplacement class was introduced to make it easier to generate differential displacement maps and to allow for more options with different materials systems by separating the differential displacement calculation from the plotting. The class is meant to replace the old differential_displacement function.
5.1. DifferentialDisplacement initialization/solve
Computing the differential displacement requires comparing the relative position vectors between a defect-free base system and a defect system for all pairs of neighbor atoms. Therefore, the calculation requires two compatible systems and a list of neighboring atoms.
system_0 (atomman.system) The defect-free base system to use.
system_1 (atomman.system) The defect system to use.
neighbors (atomman.NeighborList, optional) The neighbor list to use.
cutoff (float, optional) Cutoff distance for computing a neighbor list. If reference = 0, then system_0 will be used to generate the list. If reference = 1, then system_1 will be used to generate the list.
reference (int, optional) Indicates which of the two systems should be used for the plotting reference: 0 or 1. If 0, then system_0’s atomic positions will be used for the calculation and neighbors should be for system_0. If 1 (default), then system_1’s atomic positions will be used for the calculation and neighbors should be for system_1.
Either system_0 or system_1 can be used as the reference state for identifying neighbors and plotting atomic positions.
Using system_1 as the reference state will overlay the differential displacement vectors onto the relaxed atomic positions. This provides insight for the dislocation both in terms of its atomic nature and core spreading. This is the more common representation, especially for dislocations with large edge components.
Using system_0 as the reference state will overlay the differential displacement vectors onto a perfect crystal. This indicates how the crystal needs to be deformed to introduce the dislocation.
[11]:
# Note: the reference choice corresponds to the system for which the neighbor list was computed!
neighbors = disl_system.neighborlist(cutoff = 0.9*alat)
dd = am.defect.DifferentialDisplacement(base_system, disl_system, neighbors=neighbors, reference=1)
# or
#neighbors = base_system.neighborlist(cutoff = 0.9*alat)
#dd = am.defect.DifferentialDisplacement(base_system, disl_system, neighbors=neighbors, reference=0)
5.2. DifferentialDisplacement plotting
DifferentialDisplacement.plot() creates a differential displacement map as a matplotlib figure. There are many options to help control what is plotted as well as give you options for tweaking how the plots look.
component (str or array-like object) Indicates the component(s) of the differential displacement to plot. Values of ‘x’, ‘y’, or ‘z’ will plot the component along that Cartesian direction. A value of ‘projection’ will plot the differential displacement vectors as projected onto the plotting plane, thereby showing the two components perpendicular to the line direction. If a 3D vector is given, then the component parallel to that direction will be used.
ddmax (float or None) The maximum differential displacement value allowed. Values will be kept between +-ddmax by wrapping values with larger absolute values around by adding/subtracting 2*ddmax. Typically, this is set to beb|/2, but can be defect-specific. For instance, fcc a/2<110> dislocations and basal hcp dislocations are typically plotted with ddmax=|b|/4. If set to None, then no wrapping is done.
plotxaxis (str or array-like object, optional) Indicates the Cartesian direction associated with the system’s atomic coordinates to align with the plotting x-axis. Values are either 3D unit vectors, or strings ‘x’, ‘y’, or ‘z’ for the Cartesian axes directions. plotxaxis and plotyaxis must be orthogonal. Default value is ‘x’ = [1, 0, 0].
plotyaxis (str or array-like object, optional) Indicates the Cartesian direction associated with the system’s atomic coordinates to align with the plotting y-axis. Values are either 3D unit vectors, or strings ‘x’, ‘y’, or ‘z’ for the Cartesian axes directions. plotxaxis and plotyaxis must be orthogonal. Default value is ‘y’ = [0, 1, 0].
xlim (tuple, optional) The minimum and maximum coordinates along the plotting x-axis to include in the fit. Values are taken in the specified length_unit. If not given, then the limits are set based on min and max atomic coordinates along the plotting axis.
ylim (tuple, optional) The minimum and maximum coordinates along the plotting y-axis to include in the fit. Values are taken in the specified length_unit. If not given, then the limits are set based on min and max atomic coordinates along the plotting axis.
zlim (tuple, optional) The minimum and maximum coordinates normal to the plotting axes (i.e. plotxaxis X plotyaxis) to include in the fit. Values are taken in the specified length_unit. The optimum zlim should encompass only a single periodic slice. If not given, then the limits are set based on min and max atomic coordinates along the axis.
arrowscale (float, optional) Scaling factor for the magnitude of the differential displacement arrows. Default value is 1: no scaling, vectors are in units of length. For major components, this is often set such that the max differential displacement component after wrapping (see normfactor) is scaled to the distance between the atom pairs in the plot. For minor components, this is often set to a large value simply to make the components visible.
arrowwidth (float, optional) Scaling factor to use for the width of the plotted arrows. Default value is 0.005 = 1/200.
use0z (bool, optional) If False (default), the z coordinates from the reference system will be used for zlim and atomcmap colors. If True, the z coordinates will be used from system0 even if system1 is the reference system.
atomcolor (str or list, optional) Matplotlib color name(s) to use to display the atoms. If str, that color will be assigned to all atypes. If list, must give a color value or None for each atype. Default value (None) will use cmap instead. Note: atomcolor and atomcmap can be used together as long as exactly one color or cmap is given for each unique atype.
atomcmap (str or list, optional) Matplotlib colormap name(s) to use to display the atoms. Atoms will be colored based on their initial positions and scaled using zlim. If str, that cmap will be assigned to all atypes. If list, must give a cmap value or None for each atype. Default value (None) will use ‘hsv’ cmap. Note: atomcolor and atomcmap can be used together as long as exactly one color or cmap is given for each unique atype.
atomsize (float, optional) The circle radius size to use for the plotted atom positions in units of length. Default value is 0.5.
figsize (float or tuple, optional) Specifies the size of the figure to create in inches. If a single value is given, it will be used for the figure’s width, and the height will be scaled based on the xlim and ylim values. Alternatively, both the width and height can be set by passing a tuple of two values, but the plot will not be guaranteed to be “regular” with respect to length dimensions.
Returns
(matplotlib.Figure) The generated figure. This is returned to allow users to further modify it after creation.
5.2.1. Common plotting variations
Plot has a large number of parameters to allow for a wide range of variations and to give users the ability to tinker with the plots to make them look better if needed.
As the differential displacement values are 3D vectors, the plots focus on showing one or two components
‘x’, ‘y’, and ‘z’ are the components along the three Cartesian directions.
Any other Cartesian vector can be specified to see components along it. It is especially useful to be able to plot the component parallel to the Burgers vector for mixed dislocations.
‘projection’ plots the projection of the differential displacement vectors onto a plane perpendicular to the dislocation’s line direction. Unlike the other options, the direction of the arrows correspond to the vector values rather than the positions of the two atoms for which the vector was computed for.
Most of the other method parameters are associated with making the plots look good. The best suggestion is to create one plot, then tinker with the plotting parameters to improve how it looks.
Define common plotting parameters.
[12]:
ddmax = np.linalg.norm(burgers) / 4 # a/2<110> fcc dislocations useb|/4
# Set dict of keyword parameter values (just to make settings same for all plots below)
params = {}
params['plotxaxis'] = 'y'
params['plotyaxis'] = 'z'
params['xlim'] = (-20, 20)
params['ylim'] = (-5, 5)
params['zlim'] = (-0.01, alat*6**0.5 / 2 + 0.01) # Should be one periodic width (plus a small cushion)
params['figsize'] = 10 # Only one value as the other is chosen to make plots "regular"
params['arrowwidth'] = 1/50 # Made bigger to make arrows easier to see
params['arrowscale'] = 2.4 # Typically chosen to make arrows of length ddmax touch the corresponding atom circles
5.2.2. The Cartesian components
[13]:
dd.plot('x', ddmax, **params)
plt.title('x component')
plt.show()
dd.plot('y', ddmax, **params)
plt.title('y component')
plt.show()
dd.plot('z', ddmax, **params)
plt.title('z component')
plt.show()
Same thing, but with all components subplots of the same figure. Since the figure is defined outside the plot function the figure size will need to be manually changed to obtain an optimum appearance.
[14]:
fig = plt.figure(figsize=(10, 9.1))
ax1, ax2, ax3 = fig.subplots(3, 1, sharex=True)
ax1.set_title('x component')
ax2.set_title('y component')
ax3.set_title('z component')
dd.plot('x', ddmax, matplotlib_axes=ax1, **params)
dd.plot('y', ddmax, matplotlib_axes=ax2, **params)
dd.plot('z', ddmax, matplotlib_axes=ax3, **params)
plt.show()
5.2.3. Parallel to Burgers (or any other Cartesian vector)
[15]:
# Set ddmax to be 1/4 the magnitude of the Burgers vector (see above)
ddmax = np.linalg.norm(burgers)/4
dd.plot(burgers, ddmax, **params)
plt.title('Parallel to Burgers')
plt.show()
5.2.4. Radial
[16]:
# Set ddmax to be 1/4 the magnitude of the Burgers vector (see above)
ddmax = np.linalg.norm(burgers)/4
dd.plot('projection', ddmax, **params)
plt.title('In-plane projection')
plt.show()
5.3. Raw differential displacement data
The raw data used to generate the plots can also be retrieved as attributes of the DifferentialDisplacement objects.
ddvectors are the computed differential displacement vectors for each pair of neighbor atoms. Note that these are the full raw values computed during the class initialization and therefore have not been normalized with ddmax.
arrowcenters are the center positions for each of the differential displacement vector arrows. These points also correspond to the midpoints between the pair of neighbor atoms used to compute the differential displacement vector.
arrowuvectors are the unit vectors indicating the direction of the differential displacement vector arrows used by all component options except for ‘projection’.
[17]:
print(dd.ddvectors.shape)
print(dd.ddvectors)
(387340, 3)
[[ 0.00030753 0.00574007 0.0056176 ]
[ 0.00015473 0.00286493 0.00279437]
[-0.0001565 -0.00285422 -0.00276553]
...
[-0.00012418 0.00269864 -0.00278591]
[-0.00012418 0.00269864 -0.00278591]
[-0.00025094 0.00539228 -0.00554786]]
[18]:
print(dd.arrowcenters.shape)
print(dd.arrowcenters)
(387340, 3)
[[ 2.4866357 -112.592277 -109.13996141]
[ 1.24650369 -113.30966101 -109.14137303]
[ 3.72645929 -114.74441346 -109.14415298]
...
[ 5.38094613 108.99060147 -118.48030553]
[ 7.86105735 108.99060147 -118.48030553]
[ 6.62093836 109.70789473 -118.48168651]]
[19]:
print(dd.arrowuvectors.shape)
print(dd.arrowuvectors)
(387340, 3)
[[ 1.07171033e-04 9.99998078e-01 1.95767095e-03]
[-8.65578157e-01 5.00772906e-01 9.75316053e-04]
[ 8.65579630e-01 -5.00770378e-01 -9.65254013e-04]
...
[-8.65627722e-01 5.00687230e-01 -9.72312841e-04]
[ 8.65605989e-01 5.00724801e-01 -9.72385803e-04]
[-8.74616512e-05 9.99998127e-01 -1.93360212e-03]]