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:
Identification of the position of the defect in the simulation cell.
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}")
---------------------------------------------------------------------------
ImportError Traceback (most recent call last)
Cell In[2], line 1
----> 1 finder = DefectSiteFinder()
2 bulk_locpot = Locpot.from_file(DEFECT_DIR / "bulk/LOCPOT.gz")
3 freysoldt_results = dict()
File ~/work/pymatgen-analysis-defects/pymatgen-analysis-defects/pymatgen/analysis/defects/finder.py:61, in DefectSiteFinder.__init__(self, symprec, angle_tolerance)
59 if SOAP is None:
60 msg = "dscribe is required to use DefectSiteFinder. Install with ``pip install dscribe``."
---> 61 raise ImportError(
62 msg,
63 )
64 self.symprec = symprec
65 self.angle_tolerance = angle_tolerance
ImportError: dscribe is required to use DefectSiteFinder. Install with ``pip install dscribe``.
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}")