Freysoldt Correction

The canonical tool for performing Freysoldt correction is the sxdefectalign tool shipped as part of the SPHInX DFT code. The standalone download for the tool can be found here and the source code for the SPHInX package can be found here.

Here, we demonstrate a python implementation that removes the remaining parts of the sxdefectalign that requires user inputs:

1. Identification of the position of the defect in the simulation cell.

2. Identification of the plateau for potential alignment.

This document serves as both a validation of the python code against sxdefectalign as well as a tutorials for users of sxdefectalign to migrate their workflow over to the pure python implementation.

To get the Freysoldt corrections energies we only need to use the LOCPOT files from the DFT calculations and the charge state for the model charge. We will use our python code to obtain the defect position in the simulation cell and potential alignment and use those as inputs in the sxdefectalign tool to compare the results against our pure python implementation.

The raw data from the DFT calculations as well as the script used to generate the correction data can be found in the /tests/test_files/v_N_GaN.

from pathlib import Path

import numpy as np
from matplotlib import pyplot as plt
from monty.io import reverse_readfile
from pymatgen.analysis.defects.corrections.freysoldt import (
    Locpot,
    get_freysoldt_correction,
    plot_plnr_avg,
)
from pymatgen.analysis.defects.finder import DefectSiteFinder

TEST_FILES = Path("../../../tests/test_files")
DEFECT_DIR = TEST_FILES / "v_N_GaN"

First, we perform the analysis in python and print out the defect positions center and the potential alignment values align. These are used as the --center and --align parameters of the sxdefectalign tool.

The full result of correction is also computed here and saved for our later comparison with sxdefectalign results.

finder = DefectSiteFinder()
bulk_locpot = Locpot.from_file(DEFECT_DIR / "bulk/LOCPOT.gz")
freysoldt_results = dict()
for q in range(-1, 3):
    folder = DEFECT_DIR / f"q={q}"
    defect_locpot = Locpot.from_file(folder / "LOCPOT.gz")
    center = finder.get_defect_fpos(
        defect_structure=defect_locpot.structure, base_structure=bulk_locpot.structure
    )
    center -= np.round(center)
    print(f"Computing Freysoldt Correction for: q={q}")
    print(f"center={center}")
    freysoldt_results[q] = get_freysoldt_correction(
        q=q,
        dielectric=5,
        bulk_locpot=bulk_locpot,
        defect_locpot=defect_locpot,
        defect_frac_coords=center,
    )
    align = np.mean(list(freysoldt_results[q].metadata["alignments"].values()))
    print(f"alignment correction: {align}")

Computing Freysoldt Correction for: q=-1
center=[0.16666602 0.33333386 0.16692487]
alignment correction: 0.14016700517571093

Computing Freysoldt Correction for: q=0
center=[0.16666664 0.33333324 0.1657773 ]
alignment correction: -0.08398513278652599

Computing Freysoldt Correction for: q=1
center=[0.16666658 0.33333331 0.16008467]
alignment correction: -0.32682064580418513

Computing Freysoldt Correction for: q=2
center=[-0.01832584  0.39998153  0.22418941]
alignment correction: -0.6271082552052011

The correction results from sxdefectalign was generated by running the command at the top of the corr.out and corr_align.out files in the /tests/test_files/v_N_GaN/q=* directories. THe parameters of sxdefectalign such as the --center and --align arguments are obtained above and the results are stored in the corr.out (not accounting for alignment) and corr_align.out (with potential alignment) files. To verify the relative alignment values in the plots below, we only kept the vline-eV-a#.dat files from the unaligned called.

Now we can plot the planar averaged potential values from both methods. The documentation for sxdefectalign suggests running xmgrace -nxy vline-eV-a#.dat for each direction and visually identify the plateaus The plot_sxd_output below simulates the xmgrace -nxy call and we draw a vertical horizontal line at the automatically determined plateau value to see if they match visually.

def plot_sxd_output(q, direction, C, ax):
    """Read the vline-eV-a#.dat files computed by the `sxdefectalign` program.

    Args:
        q (int): the charge (electron is negative)
        direction (int): the direction of the line (0, 1, or 2)
        C (float): the plateau value for the short-range correction.
        ax (matplotlib.axes.Axes): the axes to plot on
    """
    v_line_file = DEFECT_DIR / f"q={q}/vline-eV-a{direction}.dat"
    # split the file into blocks based on the location of the string "&"
    with open(v_line_file, "r") as f:
        blocks = f.read().split("&")

    # convert each block into a list of lines and convert to a numpy array
    blocks = [[line.split("\t") for line in block.splitlines()] for block in blocks]
    blocks[1].pop(0)  # remove the first line of the second block
    blocks = [np.array(block, dtype=float) for block in blocks]
    colors_gen = iter(["black", "red", "green"])
    legends = iter(
        ["long range from model", "DFT Locpot diff", "short range (not aligned)"]
    )
    for b in blocks:
        for col in range(1, b.shape[1]):
            ax.plot(b[:, 0], b[:, col], c=next(colors_gen), label=next(legends))

    ax.axhline(C, color="k", linestyle="--")
    ax.text(4, C, f"C={C:0.3f}", va="bottom")
    ax.legend()


def comparison_plots(q):
    """"""
    fig, axs = plt.subplots(3, 2, figsize=(15, 15), sharey=True)
    for direction in range(3):
        plot_sxd_output(
            q,
            direction,
            freysoldt_results[q].metadata["alignments"][direction],
            axs[direction, 1],
        )
        plot_plnr_avg(
            freysoldt_results[q].metadata["plot_data"][direction], ax=axs[direction, 0]
        )
        axs[direction, 0].set_title(f"q={q} direction={direction} (Python)")
        axs[direction, 1].set_title(f"q={q} direction={direction} (SXDefectAlign)")
        axs[direction, 0].set_xlabel("")
        axs[direction, 1].set_xlabel("")

The results from sxdefectalign are shown in the right column. The peak in the LOCPOT difference corresponds to the position of the defect. We find that in the region furthest away from that peak, the plateau value computed by Python matched the plateau value visually. Note that in the case of the Python-only plots the defect position is set to the origin and the different alignment terms are set to zero far from the origin.

We recommend checking the results of your Freysoldt correction by running something like

plot_plnr_avg(
    result.metadata["plot_data"][0],
)

This immediately lets you know that, the defect positions were identified correctly if the LOCPOT diff is peaked at the origin and the curves are relatively flat in the sampling region.

Note that in the example here, the simulation cell is far too small for standard defect calculations but the defect finder is still able to find the defect position automatically. Expect it to behave better in larger simulation cells.

comparison_plots(2)

Result from larger simulation cell The LOCPOT data is too large to be included in the repository. For a 144-atom cell, the positions of the defect is visually obvious.

After accounting for the electrostatic and potential alignment corrections we can see how the Python correction compares with the results of sxdefectalign

def get_sxd_result(q):
    """Read the output of sxdefectalign"""
    fname = str(DEFECT_DIR / f"q={q}/corr_align.out")
    gen = reverse_readfile(fname)
    next(gen)
    line = next(gen)
    return float(line.split()[3])

for q in range(-1, 3):
    python_res = freysoldt_results[q].correction_energy
    sxd_res = get_sxd_result(q)
    print(f"python: {python_res:0.4f} / sxdefectalign: {sxd_res:0.4f}")

python: 0.3664 / sxdefectalign: 0.3666
python: 0.0000 / sxdefectalign: 0.0000
python: 0.1798 / sxdefectalign: 0.1799
python: 0.7721 / sxdefectalign: 0.7728