Plot Customisation with doped

All the plotting functions in doped are made to be as customisable as possible, and also return the Matplotlib Figure object, so that you can further customise the plot as you see fit.

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

CdTe Formation Energy Plotting

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

Basic formation energy plot:

plot = CdTe_thermo.plot(limit="Te-rich")  # plot at the Te-rich chemical potential limit

# run this cell to see the possible arguments for this function (or go to the python API documentation)
CdTe_thermo.plot?

dist_tol

In the above plot, we see doped classified Te interstitials into having two separate defect sites. This is dictated by the dist_tol parameter in DefectThermodynamics (= 1.5 Å by default), which groups together defects which have distances between equivalent defect sites within this tolerance. In this case, this is because Te interstitials adopt split-interstitial dimer structures for the +1 and neutral charge states, but a different conventional interstitial type structure for the +2 charge state.
Let’s increase it to 2 Å here to group these Te interstitials together:

CdTe_thermo.dist_tol = 2
plot = CdTe_thermo.plot(limit="Te-rich")  

If we had many inequivalent defects in a system (e.g. in low-symmetry/complex systems such as Sb₂Se₃), we can set dist_tol to a high value to merge together these many inequivalent defects so that our formation energy plot just shows the lowest energy species of each defect type. Let’s quickly look at monoclinic Sb₂O₅ (a promising Sb(V)-based transparent conducting oxide, with several inequivalent sites) as an example case of this:

sb2o5_thermo = loadfn("../tests/data/Sb2O5/sb2o5_thermo.json")  # load our pre-computed Sb2O5 DefectThermodynamics

fig = sb2o5_thermo.plot(limit="O-poor")

Let’s set dist_tol to a large value, which will merge all inequivalent O vacancies/interstitials and Sb interstitials so that only the lowest energy states are shown:

sb2o5_thermo.dist_tol = 10
fig = sb2o5_thermo.plot(limit="O-poor")

all_entries

We can set all_entries = True to show the full formation energy lines for all defect species, not just the lowest energy charge states. This can be useful for analysing metastable defect states, but can also be quite messy when we are plotting many defects at once. For this case, let’s trim down to just the Cd vacancy defects:

from doped.thermodynamics import DefectThermodynamics
v_Cd_thermo = DefectThermodynamics(
    [entry for entry in CdTe_thermo.defect_entries if "v_Cd" in entry.name], 
    chempots=CdTe_thermo.chempots
)  # only Cd vacancy defects

fig = v_Cd_thermo.plot(all_entries=True, limit="Te-rich")

We can instead set all_entries = "faded" to show the full formation energy lines for all defect species, but with all metastable states faded out in grey:

fig = v_Cd_thermo.plot( 
    auto_labels=True,
    xlim=(-0.5, 2.1),
    all_entries="faded",
    limit="Te-rich",
)
ax = fig.gca()
ax.axvline(0.47, ls="--", c="C4", alpha=0.4, zorder=-1)  # add a vertical line at the transition level

<matplotlib.lines.Line2D at 0x282d232b0>

chempots & el_refs

We can also set the chemical potentials to be used in the plot directly in the plotting function (or can of course adjust DefectThermodynamics.chempots/el_refs directly):

# we can always set/overwrite `chempots` in the plotting/tabulation functions like this too:
def_plot = v_Cd_thermo.plot(
    chempots = {"Cd": -0.6255, "Te": -0.625},  # midpoint between Te/Cd-rich  
    el_refs=CdTe_thermo.chempots.get("elemental_refs"),
    auto_labels=True,
)

limit

The limit parameter specifies the chemical potential limit for which to plot formation energies. This can be either:

• None, in which case plots are generated for all limits in chempots.

• “X-rich”/”X-poor” where X is an element in the system, in which case the most X-rich/poor limit will be used (e.g. “Li-rich”).

• A key in the (DefectThermodynamics.)chempots["limits"] dictionary.

The latter two options can only be used if chempots is in the doped format (see chemical potentials tutorial).

# print the chemical potential limits
print(CdTe_thermo.chempots["limits"].keys())

dict_keys(['Cd-CdTe', 'CdTe-Te'])

In this case, we can plot the Cd-rich limit by setting limit="Cd-rich", limit="Cd-CdTe" (the Cd-rich limit here) or by just manually setting the chemical potentials to this limit:

def_plot = CdTe_thermo.plot(limit="Cd-CdTe")

# manual Cd-rich chempots:
def_plot = CdTe_thermo.plot(
    chempots=CdTe_thermo.chempots["limits_wrt_el_refs"]["Cd-CdTe"],
    el_refs=CdTe_thermo.chempots["elemental_refs"]
)

colormap

We can set different (qualitative) colour maps for the plots. This can be a string name of a Matplotlib colormap (see here), or a Colormap object:

plot = CdTe_thermo.plot(limit="Te-rich", colormap="Accent")  

from matplotlib.colors import ListedColormap

custom_cmap = ListedColormap(["#D4447E", "#E9A66C", "#507BAA", "#5FABA2", "#63666A", "#A3A3A3", "#FFD700"])

plot = CdTe_thermo.plot(limit="Te-rich", colormap=custom_cmap)  

style_file

We can adjust the overall style of the plot by using a custom matplotlib style (mplstyle) file:

with open("custom_style.mplstyle", "w") as f:
    f.write("ytick.right : True\nxtick.top : True\nfont.sans-serif : Helvetica")

plot = CdTe_thermo.plot(limit="Te-rich", style_file="custom_style.mplstyle")  

auto_labels

We can include automatic labels for the transition levels: (as you can see, these get very messy when we are plotting many defects, so usually best to use only with a small number of defects)

plot = CdTe_thermo.plot(limit="Te-rich", auto_labels=True)  

chempot_table

We can control whether or not to show the chemical potentials above the plot:

plot = CdTe_thermo.plot(limit="Te-rich", chempot_table=False)  

xlim & ylim

We can adjust the axis limits:

plot = CdTe_thermo.plot(limit="Te-rich", xlim=(-0.5, 2.1), ylim=(-0.5, 4.0))

fermi_level

We can show a vertical line at the predicted Fermi level position (see the doped thermodynamics tutorial for details on how we might go about doing this):

plot = CdTe_thermo.plot(limit="Te-rich", fermi_level=0.4)

filename

We can set filename to save the plot to file:

plot = CdTe_thermo.plot(limit="Te-rich", filename="CdTe_defects_plot.png")  

Finite-Size Charge Correction Plots

Both the Freysoldt (FNV) and Kumagai (eFNV) charge correction plots in doped are also quite customisable, and as always return the Matplotlib Figure object to allow further customisation.

eFNV (Kumagai) Correction

As shown in the doped docs Tips page, we can set defect_region_radius and/or excluded_indices with the eFNV correction (get_kumagai_correction()) to control which sites away from the defect site are used to determine the potential alignment – which we may want to do with e.g. low-dimensional materials.

te_cd_entry = [entry for entry in CdTe_thermo.defect_entries if entry.name == "Te_Cd_+1"][0]

correction, plot = te_cd_entry.get_kumagai_correction(plot=True)

Calculated Kumagai (eFNV) correction is 0.238 eV

As with the defect formation energy plots, because this function also returns the Matplotlib Figure object, we can further customise the plot as we see fit. To illustrate, here we remove the legend from the plot and shade the region where the potential is positive in yellow:

correction, plot = te_cd_entry.get_kumagai_correction(plot=True)
ax = plot.gca()  # get axis object
ax.get_legend().remove()  # remove legend
ax.axhspan(0, 100, alpha=0.2, color="yellow")

Calculated Kumagai (eFNV) correction is 0.238 eV

<matplotlib.patches.Polygon at 0x17da912d0>

style_file

As with the defect formation energy plots, we can adjust the overall style of the plot by using a custom matplotlib style (mplstyle) file:

with open("custom_style.mplstyle", "w") as f:
    f.write("ytick.right : True\nxtick.top : True\nfont.sans-serif : Helvetica")

correction, plot = te_cd_entry.get_kumagai_correction(plot=True, style_file="custom_style.mplstyle")  

Calculated Kumagai (eFNV) correction is 0.238 eV

FNV (Freysoldt) Correction

Remember, the Freysoldt (FNV) correction is only valid for isotropic systems!

v_cd_entry = [entry for entry in CdTe_thermo.defect_entries if entry.name == "v_Cd_-2"][0]

correction, plot = v_cd_entry.get_freysoldt_correction(plot=True)

Calculated Freysoldt (FNV) correction is 0.738 eV

For the FNV correction, it calculates the electrostatic potential alignment along the a, b and c axes and then computes the average. We can focus on just one axis at a time by setting axis:

correction_a, plot_a = v_cd_entry.get_freysoldt_correction(plot=True, axis=0)

Calculated Freysoldt (FNV) correction is 0.738 eV

Again, because this function also returns the Matplotlib Figure object, we can further customise the plot as we see fit. To illustrate, here we remove the legends from each plot and shade the region where the potential is positive in yellow:

correction, plot = v_cd_entry.get_freysoldt_correction(plot=True)
for ax in plot.get_axes():
    ax.get_legend().remove()  # remove legend
    ax.artists[0].remove()  # shade the region where the potential is positive
    ax.axhspan(0, 100, alpha=0.2, color="yellow")

Calculated Freysoldt (FNV) correction is 0.738 eV

style_file

As with the defect formation energy plots, we can adjust the overall style of the plot by using a custom matplotlib style (mplstyle) file:

with open("custom_style.mplstyle", "w") as f:
    f.write("ytick.right : True\nxtick.top : True\nfont.sans-serif : Helvetica")

correction, plot = v_cd_entry.get_freysoldt_correction(plot=True, style_file="custom_style.mplstyle")  

Calculated Freysoldt (FNV) correction is 0.738 eV

Site Displacement / Strain Plots

As shown in more detail in the advanced analysis tutorial, the DefectEntry.plot_site_displacements() method and functions in doped.utils.displacements can be quite useful for analysing the lattice response (site displacements) and local strain for a relaxed defect structure.

Again, these functions return the Matplotlib Figure object to allow further customisation (more examples shown in the advanced analysis tutorial):

te_cd_entry = [entry for entry in CdTe_thermo.defect_entries if entry.name == "Te_Cd_+1"][0]
plot = te_cd_entry.plot_site_displacements()
ax = plot.gca(); ax.set_xlim(2, 11); ax.set_ylim(0, 0.22);  # adjust x/y-lims

Eigenvalue Plots

As shown in more detail in the advanced analysis tutorial, the DefectEntry.get_eigenvalue_analysis() method and functions in doped.utils.eigenvalues can be quite useful for analysing the electronic structure of our defect calculations.

Again, these functions return the Matplotlib Figure object to allow further customisation:

te_cd_entry = [entry for entry in CdTe_thermo.defect_entries if entry.name == "Te_Cd_+1"][0]
bes, fig = te_cd_entry.get_eigenvalue_analysis()
# adjust ylim and add line around defect level:
ax = fig.axes[0]; ax.set_ylim(-0.2, 1.8), ax.axhline(1.1);