Modelling the Mobility of Dislocations with matscipy

The dislocation modelling toolkit provided in matscipy.dislocation includes functions to assist in characterising the mobility of dislocations in various systems. Currently, only motion in the glide direction is supported.

Dislocation Glide

Dislocation Glide is where the dislocation line migrates in the glide direction. There are two main mechanisms for attempting to describe dislocation glide dynamics.

The first mechanism is the (double) kink mechanism. Starting from a perfectly straight dislocation line, a small segment nucleates out in the glide directionforming a pair of dislocation kinks. These kinks can then migrate along the dislocation line.

A second, simpler mechanism is to consider the entire dislocation line moving at once - this is much less physical than the double kink mechanism, but can be a good first approximation. The 2D nature of this mechanism can also make it much simpler to study, rather than having to deal with the fully 3D kinks.

Glide in Dislocation Cylinders

As an example of modelling dislocation glide in cylindrical cells, let’s look at the Diamond Screw dislocation in Si, using the Holland and Marder potential:

from matscipy.dislocation import get_elastic_constants, DiamondGlide60Degree
# the calculator to provide forces and energies from the potential
from matscipy.calculators.manybody.explicit_forms.stillinger_weber import StillingerWeber,\
                                                                Holland_Marder_PRL_80_746_Si
from matscipy.calculators.manybody import Manybody
calc = Manybody(**StillingerWeber(Holland_Marder_PRL_80_746_Si))
import numpy as np

# the function accepts any ASE type of calculator
alat, C11, C12, C44 = get_elastic_constants(calculator=calc, symbol="Si", verbose=False)

Si_disloc = DiamondGlide60Degree(alat, C11, C12, C44, symbol="Si")

Si_disloc_bulk, Si_disloc_dislo = Si_disloc.build_cylinder(radius=20)
Si_disloc.view_cyl(Si_disloc_dislo, mode="dxa")

/usr/lib/python3/dist-packages/scipy/__init__.py:146: UserWarning: A NumPy version >=1.17.3 and <1.25.0 is required for this version of SciPy (detected version 1.26.4
  warnings.warn(f"A NumPy version >={np_minversion} and <{np_maxversion}"

/usr/local/lib/python3.10/dist-packages/matscipy/dislocation.py:485: FutureWarning: Import StrainFilter from ase.filters
  sf = StrainFilter(unit_cell)

To generate the initial and final structures for a full glide event, we can use the build_glide_configurations function. The structures contain straight dislocation lines, separated by the dislocation glide distance. These structures will be similar to those produced by a call to build_cylinder, except extra bulk is added (creating a pill shape, not a circle) in order to make the initial and final configurations symmetric.

We can also then use the ase.neb tools to smoothly interpolate an approximate glide path, which allows us to generate structures for the simpler dislocation glide mechanism discussed above.

bulk, glide_ini, glide_fin = Si_disloc.build_glide_configurations(radius=20)

from ase.neb import NEB

nims = 5

glide_neb = NEB([glide_ini.copy() for i in range(nims-1)] + [glide_fin.copy()])

glide_neb.interpolate(method="idpp", apply_constraint=True)

/tmp/ipykernel_8330/2158173092.py:7: FutureWarning: Please import NEB from ase.mep, not ase.neb.
  glide_neb = NEB([glide_ini.copy() for i in range(nims-1)] + [glide_fin.copy()])

To visualise the glide structures, we will combine ase’s plot_atoms to convert a structure to a matplotlib plot, and then use FuncAnimation to animate the glide structures:

def animate_glide(images, diamond=True, radii=None):
    from ase.visualize.plot import plot_atoms
    import matplotlib.pyplot as plt
    from matplotlib.animation import FuncAnimation
    from matscipy.utils import get_structure_types
    from visualisation import show_HTML

    fig, ax = plt.subplots(figsize=(10, 10))
    ax.axis("off")

    # Add extra reps of start and end points for clarity
    anim_images = [images[0]] * 3 + images + [images[-1]] * 3

    def plot_frame(framedata):
        ax.clear()
        # Plot an individual frame of the animation 
        framenum, atoms = framedata

        # get CNA colours to enhance plot
        atom_labels, struct_names, colors = get_structure_types(atoms, 
                                                                diamond_structure=diamond)
        atom_colors = [colors[atom_label] for atom_label in atom_labels]

        plot_atoms(atoms, ax=ax, colors=atom_colors, radii=radii)


    animation = FuncAnimation(fig, plot_frame, frames=enumerate(anim_images),
                                save_count=len(anim_images),
                                init_func=lambda: None,
                                interval=200)
    
    # Need to convert animation to HTML in order for it to be visible on the docs
    return show_HTML(animation)

animate_glide(glide_neb.images, radii=0.3)

Once Loop Reflect

This is also possible in the dissociated case:

bulk, glide_ini, glide_fin = Si_disloc.build_glide_configurations(radius=20, partial_distance=5)

nims = 5

glide_neb = NEB([glide_ini.copy() for i in range(nims-1)] + [glide_fin.copy()])

glide_neb.interpolate(method="idpp", apply_constraint=True)

animate_glide(glide_neb.images, radii=0.3)

/tmp/ipykernel_8330/2427532584.py:5: FutureWarning: Please import NEB from ase.mep, not ase.neb.
  glide_neb = NEB([glide_ini.copy() for i in range(nims-1)] + [glide_fin.copy()])

Once Loop Reflect

Glide in Dislocation Quadrupoles

As quadrupoles are extremely useful for modelling dislocations using plane-wave DFT, it can be convenient to be able to set up initial guesses for complex processes such as dislocation glide.

We can use the function build_glide_quadrupoles to construct a set of images for this system, which can optionally model the glide of either the “left” (\(+\mathbf{b}\)) or “right” (\(-\mathbf{b}\)) dislocation cores, or both at once.

from matscipy.dislocation import Quadrupole

Si_quad = Quadrupole(DiamondGlide60Degree, alat, C11, C12, C44, symbol="Si")

num_images = 5

glide_quads = Si_quad.build_glide_quadrupoles(nims=num_images, 
                                            glide_left=True, # Allow left dislocation glide
                                            glide_right=True, # Allow right dislocation glide
                                            glide_separation=6,
                                            verbose=False)

animate_glide(glide_quads)

Once Loop Reflect

num_images = 5

glide_quads = Si_quad.build_glide_quadrupoles(nims=num_images, 
                                            glide_left=False, # Prevent left dislocation glide
                                            glide_right=True, # Allow right dislocation glide
                                            partial_distance=2, # Dissociated Glide
                                            glide_separation=8,
                                            verbose=False)

animate_glide(glide_quads)

Once Loop Reflect

Dislocation Kink

Dislocation kink is often the preferred mechanism for migration of the dislocation line. It involves a short segment of the dislocation line hopping by one glide vector, which then provides a low barrier for the rest of the dislocation line to migrate.

Here the space of structures to explore is greater, due to the 3D nature of dislocation kink. Kink nucleation is also possible on segments of an already kinked-out dislocation line, which can lead to a full network of dislocation kinks.

Kink maps

Dislocation kink in matscipy is controlled by a kink map, which controls the position of the dislocation line as a function of the vertical coordinate. The kink map is a list of integers, where the values give a line position in units of glide distances, and each element corresponds to an extra cell in z (replication of disloc.unit_cell).

A kink map of [0, 0, 1, 1] means that the dislocation line initially starts at position zero for two repetitions, and then moves out by one glide distance for two repetitions.

Both dislocation cylinders and quadrupoles have support for this kind of kink map, but the periodic boundaries are treated differently.

In dislocation cylinders (using build_kink_cyl), the periodicity in z is enforced by requiring that the dislocation line returns back to the initial position across the periodic boundary. The example kink map of [0, 0, 1, 1] will therefore have a single kink out in the center of the cell, and a single kink back at the periodic boundary (the dislocation line will go like 0, 0, 1, 1, 0, 0, 1, 1, ...).

In quadrupoles (using build_kink_quadrupole_network), the cell is modified such that the dislocation position at the top of the cell becomes the new dislocation line position at the bottom of the cell. This means that for quadrupoles an extra kink will not be created, and that the example map of [0, 0, 1, 1] will create only one kink in the center of the cell.

Since the kink map is just a list of integers, it can be more convenient to exploit list addition, and specify kink maps in a form similar to [0] * 5 + [1] * 5, which would be identical to the input of [0, 0, 0, 0, 0, 1, 1, 1, 1, 1].

Kink in cylinders

Si_bulk, Si_kink = Si_disloc.build_kink_cylinder(
                    radius=20,
                    kink_map= [0] * 5 + [1] * 5
                )

Si_disloc.view_cyl(Si_kink)

Kink in quadrupoles

Like with dislocation cylinders, we can build networks of kinks in dislocation quadrupoles:

Si_quad_bulk, Si_quad_kink = Si_quad.build_kink_quadrupole(
    glide_separation=8,
    kink_map=[0]*5 + [1]*5,
    verbose=False
)

Si_quad.view_quad(Si_quad_kink, mode="dxa")

Warning

For the cell shift to always produce the correct crystalstructure, some atoms need to be deleted for some values of kink_map[-1]. This means that some kink maps may not conserve stoichiometry for some multispecies systems.

There is also another routine build_minimal_kink_quadrupole for building only the smallest possible kink structures. This is where the kink happens in a single burgers vector cell. Here, the n_kink parameter controls the number of glide distances covered by the kink, and the direction of the kink: n_kink=2 builds a compressed version of the kink map [0, 2], and n_kink=-1 constructs a compressed version of [0, -1].

Si_quad_bulk, Si_quad_kink = Si_quad.build_minimal_kink_quadrupole(
    glide_separation=8,
    n_kink=1,
    verbose=False
)

Si_quad.view_quad(Si_quad_kink, mode="dxa")

Improving the kink structures

So far, we have only used the Continuum Linear Elastic (CLE) solutions when building kink structures. As the kink structures are built by interpolating between glide structures, we can get a better approximation of the kink if we relax these glide structures before building the kink.

To do this, we can replace a call to build_kink_cylinder with a call to build_kink_glide_structs to build the required glide structures, and then build_kink_from_glide_cyls to actually construct the kink. For quadrupoles, the equivalent is replacing build_kink_quadrupole with build_kink_quadrupole_glide_structs and build_kink_quadrupole_from_glide_structs.

from matscipy.calculators.manybody.explicit_forms.stillinger_weber import StillingerWeber,\
                                                                Holland_Marder_PRL_80_746_Si
from matscipy.calculators.manybody import Manybody
from ase.optimize.precon import PreconLBFGS


calc = Manybody(**StillingerWeber(Holland_Marder_PRL_80_746_Si))

kink_map = [0] * 5 + [1] * 5

ref_bulk, glide_structs, struct_map = Si_disloc.build_kink_glide_structs(kink_map, radius=40)

# glide_structs has a length of 2, as only two unique values in the kink map

for struct in glide_structs:
    struct.calc = calc
    opt = PreconLBFGS(struct, logfile=None)
    opt.run(1e-1)

Si_bulk, Si_kink = Si_disloc.build_kink_from_glide_cyls(ref_bulk, glide_structs, struct_map)

Si_disloc.view_cyl(Si_kink)

Kink in Dissociated Dislocations

For dissociated dislocations, kinks can nucleate at either core independently, with the energetics of the process largely determined by the burgers vectors of both partial dislocations, and by their separation. Kink can be modelled in such systems by using a 2D kink map, where the 2nd axis is the number of dislocation cores.

As these structures feature complex 3D geometry, we will hide bulk-coordinated atoms so that the dislocation cores are easily visible.

kink_map = np.array(
    #   L, R ordering (L = self.left_dislocation, R = self.right_dislocation)
      [[0, 5],] * 2 # Start dissociated
    + [[0, 6],] * 5 # Kink out the right partial by one
    + [[1, 6],] * 5 # Kink out the left partial by one
    + [[0, 5],] * 2 # Wrap back to initial position
) 
# Kink map gets wrapped at boundary anyway, but manually specifying the wrapping
# in this way results in an easier-to-read visualisation

print(kink_map.shape)

Si_bulk, Si_kink = Si_disloc.build_kink_cylinder(
                    radius=20,
                    kink_map=kink_map
                )

Si_disloc.view_cyl(Si_kink, hide_bulk=True)

(14, 2)

Bulk-terminated Loops of Dissociated Dislocations

To push these tools to their limit, let’s look at a dislocation loop enclosed in bulk, where the dislocation is able to dissociate. We can use the Quadrupole class to easily create a dislocation model with four total dislocation cores: the two partial dislocations forming the left side of the dislocation loop, and the two partials forming the right side of the loop.

To do this, we specify 4 kink positions per vertical segment, and we construct the loop in vertical slices.

# Quadrupole & CubicCrystalDissociatedDislocation have self.left_dislocation and 
# self.right_dislocation attributes
# Combining them in this way results in an ordering of:
# LL, LR, RL, RR, where LL is self.left_dislocation.left_dislocation, 
# RL is self.right_dislocation.left_dislocation, etc 


kink_map = np.array(
    #   LL, LR, RL, RR ordering
      [[ 0,  0,  0,  0]] * 7 # Start with perfect bulk (as dislocation cores overlapping)
    + [[-1, -1,  1,  1]] * 5 # Kink out both perfect dislocations
    + [[-2, -2,  2,  2]] * 4 
    + [[-3, -3,  3,  3]] * 4
    + [[-4, -3,  3,  4]] * 4 # Start dissociation of both perfect dislocations
    + [[-5, -3,  3,  5]] * 4
    + [[-6, -3,  3,  6]] * 8 # Full extent of the dislocation loop
    + [[-5, -3,  3,  5]] * 4 # Invert the process, moving cores closer together
    + [[-4, -3,  3,  4]] * 4
    + [[-3, -3,  3,  3]] * 4
    + [[-2, -2,  2,  2]] * 4 
    + [[-1, -1,  1,  1]] * 5
    + [[ 0,  0,  0,  0]] * 7 # End with perfect bulk
)

Si_bulk, Si_kink = Si_quad.build_kink_cylinder(
                    radius=50,
                    kink_map=kink_map
                )

Si_disloc.view_cyl(Si_kink, hide_bulk=True)