Advanced Defect Analysis

Here we describe some more targeted analysis you can do for your defect calculations parsed with doped (including comparing the relaxed configurations for different initial interstitial positions, structure & bond length analysis of defects, and plotting/analysis of the defect charge corrections), which may be useful for in certain cases.

Tip

You can run this notebook interactively through Google Colab or Binder – just click the links! If running on Colab, then you’ll need to run !pip install doped in a cell to install the package, and !git clone https://github.com/SMTG-Bham/doped to download the example data (and update paths in the code cells accordingly).

Defect-Induced Site Displacements (Strain)

Using the DefectEntry.plot_site_displacements() method, we can analyse the displacements of atoms around the defect site (i.e. defect-induced local strain) during geometry relaxation.

%matplotlib inline
from monty.serialization import loadfn
CdTe_defects_thermo = loadfn("CdTe/CdTe_thermo_wout_meta.json.gz")  # load our DefectThermodynamics object

Absolute Displacements

Let’s look at the displacements of atoms around the Cd vacancy in CdTe, and how this changes with charge state:

import matplotlib.pyplot as plt
from doped.utils.plotting import format_defect_name

v_Cd_entries = [entry for entry in CdTe_defects_thermo.defect_entries.values() if "v_Cd" in entry.name]
for defect_entry in v_Cd_entries:
    fig = defect_entry.plot_site_displacements(separated_by_direction=False)
    fig.suptitle(format_defect_name(defect_entry.name, include_site_info_in_name=False), fontsize=18)
plt.show()

If we want to get an interactive plot, which can be useful for analysis, we can set use_plotly = True to return plotly versions of these plots:

for defect_entry in v_Cd_entries:  # returns interactive plots when run in a notebook!
    fig = defect_entry.plot_site_displacements(separated_by_direction=False, use_plotly=True); fig.show()

Separated by Direction

for defect_entry in v_Cd_entries:
    fig = defect_entry.plot_site_displacements(separated_by_direction=True)
    fig.suptitle(format_defect_name(defect_entry.name, include_site_info_in_name=False),
                 y=1.2, fontsize=22)
plt.show()

Here we see that \(V_{Cd}^{-2}\) has isotropic (symmetric) displacements of atoms around the vacancy site in the x/y/z directions, which makes sense as it adopts a tetrahedral (Td) geometry (as shown in the get_symmetries_and_degeneracies output below and discussed in detail in this paper). As expected, we see an exponential tail-off in the site displacement magnitudes as we move away from the defect, and it is the Te nearest neighbours which are most strongly perturbed.

For \(V_{Cd}^{-1}\), we have a Jahn-Teller-distorted \(C_{3v}\) geometry where a neighbouring Te atom displaces along the [111] direction away from the vacancy site (while the other Te atoms displace towards the vacancy but by a smaller degree), while for \(V_{Cd}^{0}\), we have a \(C_{2v}\) geometry where two neighbouring Te displace significantly towards the vacancy and each other (around 1 Å) forming a dimer bond, while the other two Te move a smaller distance towards the vacancy (around 0.2 Å) as seen in the displacement plots and illustrated below.

from doped.thermodynamics import DefectThermodynamics
v_Cd_thermo = DefectThermodynamics(v_Cd_entries)
v_Cd_thermo.get_symmetries_and_degeneracies()

		Site_Symm	Defect_Symm	g_Orient	g_Spin	g_Total	Mult
Defect	q
v_Cd	0	Td	C2v	6.0	1	6.0	1.0
	-1	Td	C3v	4.0	2	8.0	1.0
	-2	Td	Td	1.0	1	1.0	1.0

The high symmetry of \(V_{Cd}^{-2}\) is evident from the displacement plots above, where it looks like there are much fewer atoms in the plots, however this is just because we have many symmetry-equivalent atoms in this case and so we end up with many overlapping points (and so much less distinct points). Then for \(C_{3v}\) \(V_{Cd}^{-1}\) we have more distinct sites appearing, and then more so for the lower-symmetry \(C_{2v}\) \(V_{Cd}^{0}\) structure.

Relative Displacements

Instead of plotting the absolute displacements in the x, y, z directions, we can also plot the displacements of atoms relative to the defect site (i.e. displacement of the atom along the line connecting itself to the defect). A negative displacement indicates that the atom moves towards the defect (compressive strain) while a positive displacement indicates that the atom moves away from the defect (tensile strain):

for defect_entry in v_Cd_entries:
    fig = defect_entry.plot_site_displacements(relative_to_defect=True)
    fig.suptitle(format_defect_name(defect_entry.name, include_site_info_in_name=False),
                 fontsize=18)
plt.show()

The above plots clearly show the different reconstructions of each charge state:

• The formation of the neutral vacancy introduces two holes which localise in the Te-Te dimer. These Te atoms forming the dimer displace significantly towards the defect, while the two other NN Te atoms also move towards the vacancy, but to a smaller degree.

• The singly charged vacancy introduces a single hole which localises in one of the Te NNs. This Te atom displaces away from the vacancy in the (-1,-1,-1) direction, while the other 3 Te atoms displace towards the vacancy.

• The doubly charged vacancy keeps the \(T_d\) symmetry, with the 4 Te NNs displacing towards the vacancy by 0.2 Å (from the original Cd-Te bond distance of 2.83 Å to a distance of 2.61 Å) to allow for greater hybridization between dangling bonds.

Note

The data for the atomic site displacements in the relaxed defect supercell is stored in the DefectEntry.calculation_metadata["site_displacements"] attributes, which has a dictionary of site displacement vectors (relative to the unperturbed bulk positions) and their (relaxed) distances to the defect site, ordered by the atom type in the bulk formula, then by distance from the defect site.

Because of the high symmetry and thus overlapping displacement datapoints in the \(V_{Cd}^{-2}\) plot, we can also directly access the site displacement data (stored in DefectEntry.calculation_metadata["site_displacements"], or by directly using calc_site_displacements from doped.utils.displacements) to confirm the number of sites with a given displacement, or to perform other analyses:

import numpy as np

defect_entry = CdTe_defects_thermo.defect_entries["v_Cd_-2"]

# this is a dict of 'displacements' vectors and (original) 'distances' lists, 
# sorted by atom type, then by distance from the defect site:  
displacements_dict = defect_entry.calculation_metadata["site_displacements"]
Te_indices = [i for i, site in enumerate(defect_entry.defect_supercell) if site.specie.symbol == "Te"]
Te_displacements = np.array(displacements_dict["displacements"])[Te_indices]
for i, disp in enumerate(Te_displacements[:5]):
    print(f"Te atom {i} displacement vector: {disp}, magnitude: {np.linalg.norm(disp):.2f} Å")

Te atom 0 displacement vector: [ 0.127 -0.127 -0.127], magnitude: 0.22 Å
Te atom 1 displacement vector: [-0.127 -0.127  0.127], magnitude: 0.22 Å
Te atom 2 displacement vector: [-0.127  0.127 -0.127], magnitude: 0.22 Å
Te atom 3 displacement vector: [0.127 0.127 0.127], magnitude: 0.22 Å
Te atom 4 displacement vector: [ 0.00363  0.02292 -0.00363], magnitude: 0.02 Å

Nice, we can see that it is 4 Te neighbours that each displace by the same magnitude along the {111} directions by 0.22 Å toward \(V_{Cd}^{-2}\), which thus retains the \(T_d\) symmetry of this site.

Displacements Along Specific Directions

We can also analyse the displacements along a specific direction. For instance, for the \(V_{Cd}^{-1}\) defect, we can plot the displacements of atoms along the (-1, -1, -1) direction (the vector along which one of the Te NNs displaces away from the vacancy):

defect_entry = CdTe_defects_thermo.defect_entries["v_Cd_-1"]
fig = defect_entry.plot_site_displacements(vector_to_project_on=[-1,-1,-1])
_ = fig.suptitle(format_defect_name(defect_entry.name, include_site_info_in_name=False),
                y=1.2, fontsize=18)
plt.show()

which shows how one of the Te NNs displaces in the (-1,-1,-1) direction away from the vacancy, while the other 3 Te atoms displace mostly perpendicular to this direction (overlaped points in the right plot at a distance of 2.5 Å). This is also shown when orienting the defect environment along the (-1,-1,-1) direction, where we see that the Te atom behind the vacancy (in shaded black) has displaced away from the vacancy (from an original distance of 2.8 to 3.0 Å).

Similarly, for the neutral case, we can plot the displacements along the vector connecting the Te dimer ((1,1,0) vector). This will show two Te atoms significantly displacing along this vector, but in opposite directions (moving towards each other to form the dimer bond):

defect_entry = CdTe_defects_thermo.defect_entries["v_Cd_0"]  # Neutral state
vector_to_project_on = [1,1,0]  # Vector connecting Te dimer
fig = defect_entry.plot_site_displacements(vector_to_project_on=vector_to_project_on)
_ = fig.suptitle(format_defect_name(defect_entry.name, include_site_info_in_name=False),
                y=1.2, fontsize=18)
plt.show()

which indeed shows one Te moving by 1 Å along the (1,1,0) direction while the other Te moves by ~0.9 Å along the (-1,-1,0) direction, reducing the original Te-Te bond distance from 4.63 Å to 2.75 Å. This can also be visualised by orienting the defect environment along the (1,1,0) direction:

Substitution Example

For illustration purposes, here we plot the site displacements for \(Te_{Cd}^{+1}\):

Te_Cd_1 = CdTe_defects_thermo.defect_entries["Te_Cd_+1"]
fig = Te_Cd_1.plot_site_displacements()
plt.show()

We see that the antisite Te atom displaces significantly from the original Cd lattice site (on which it is substituting, at \(x\) = 0 in the plot below), and then three of the neighbouring Te displace by ~0.2 Å while the fourth Te neighbour only displaces by ~0.03 Å, resulting in the \(C_{3v}\) symmetry of this defect:

Te_Cd_1.calculation_metadata["relaxed point symmetry"]

'C3v'

Interstitial Example

Cd_i_2 = CdTe_defects_thermo.defect_entries["Cd_i_Td_Cd2.83_+2"]
fig = Cd_i_2.plot_site_displacements(relative_to_defect=True)
plt.show()

Here for \(Cd_{i}^{+2}\) (coordinated by Cd), we see that the coordinating \(Cd^{+2}\) ions displace away from the positively-charged interstitial, while the 2nd-neighbour shell of \(Te^{-2}\) anions displace toward the interstitial site, as expected. The high (\(T_d\)) symmetry of this fully-ionised interstitial (and thus atomic displacments) means only a few plotted points are visible due to overlapping.

Displacement Ellipsoid Analysis

from doped.utils.displacements import calc_displacements_ellipsoid, plot_displacements_ellipsoid
ellipsoid_plot, anisotropy_plot = plot_displacements_ellipsoid(
    CdTe_defects_thermo.defect_entries["v_Cd_0"], plot_anisotropy=True, quantile=0.8  # neutral V_Cd
)
plt.show()

The first image shows an ellipsoid with atoms inside that have a displacement greater than a certain threshold. The black dots represent atoms and the dark lines represent crystal lattices. The threshold can be adjusted by the quantile of the displacement norm. The shape of this ellipsoid shows the anisotropy of the displacement around the defect.

The left side of the second figure shows the box plot of displacement norm. The right side of the second figure plots the ratio of the radius of the longest ellipsoid to the radius of the other two. Hereafter, this graph is called an anisotropy plot. In the anisotropy plot, when a point is in the upper right region, it indicates that the displacement around the defect is isotropic. When a point is in the lower region, it indicates that the displacement around the defect has a two-dimensional spread.

Again, we can set use_plotly = True to get an interactive plot of the displacement ellipsoid:

ellipsoid_plot, anisotropy_plot = plot_displacements_ellipsoid(
    CdTe_defects_thermo.defect_entries["v_Cd_0"], use_plotly=True, plot_anisotropy=True, quantile=0.8
)
ellipsoid_plot.show()
anisotropy_plot.show()