Source code for pymatgen.analysis.defects.corrections
# coding: utf-8
# Copyright (c) Pymatgen Development Team.
# Distributed under the terms of the MIT License.
"""
Implementation of defect correction methods.
"""
import logging
import numpy as np
import scipy
from scipy import stats
from pymatgen.analysis.defects.core import DefectCorrection
from pymatgen.analysis.defects.utils import ang_to_bohr, hart_to_ev, eV_to_k, \
generate_reciprocal_vectors_squared, QModel, converge, tune_for_gamma, \
generate_R_and_G_vecs, kumagai_to_V
import matplotlib.pyplot as plt
__author__ = "Danny Broberg, Shyam Dwaraknath"
__copyright__ = "Copyright 2018, The Materials Project"
__version__ = "1.0"
__maintainer__ = "Shyam Dwaraknath"
__email__ = "shyamd@lbl.gov"
__status__ = "Development"
__date__ = "Mar 15, 2018"
logger = logging.getLogger(__name__)
[docs]class FreysoldtCorrection(DefectCorrection):
"""
A class for FreysoldtCorrection class. Largely adapated from PyCDT code
If this correction is used, please reference Freysoldt's original paper.
doi: 10.1103/PhysRevLett.102.016402
"""
def __init__(self, dielectric_const, q_model=None, energy_cutoff=520, madetol=0.0001, axis=None):
"""
Initializes the FreysoldtCorrection class
Args:
dielectric_const (float or 3x3 matrix): Dielectric constant for the structure
q_model (QModel): instantiated QModel object or None.
Uses default parameters to instantiate QModel if None supplied
energy_cutoff (int): Maximum energy in eV in reciprocal space to perform
integration for potential correction.
madeltol(float): Convergence criteria for the Madelung energy for potential correction
axis (int): Axis to calculate correction.
If axis is None, then averages over all three axes is performed.
"""
self.q_model = QModel() if not q_model else q_model
self.energy_cutoff = energy_cutoff
self.madetol = madetol
self.dielectric_const = dielectric_const
if isinstance(dielectric_const, int) or \
isinstance(dielectric_const, float):
self.dielectric = float(dielectric_const)
else:
self.dielectric = float(np.mean(np.diag(dielectric_const)))
self.axis = axis
self.metadata = {"pot_plot_data": {}, "pot_corr_uncertainty_md": {}}
[docs] def get_correction(self, entry):
"""
Gets the Freysoldt correction for a defect entry
Args:
entry (DefectEntry): defect entry to compute Freysoldt correction on.
Requires following keys to exist in DefectEntry.parameters dict:
axis_grid (3 x NGX where NGX is the length of the NGX grid
in the x,y and z axis directions. Same length as planar
average lists):
A list of 3 numpy arrays which contain the cartesian axis
values (in angstroms) that correspond to each planar avg
potential supplied.
bulk_planar_averages (3 x NGX where NGX is the length of
the NGX grid in the x,y and z axis directions.):
A list of 3 numpy arrays which contain the planar averaged
electrostatic potential for the bulk supercell.
defect_planar_averages (3 x NGX where NGX is the length of
the NGX grid in the x,y and z axis directions.):
A list of 3 numpy arrays which contain the planar averaged
electrostatic potential for the defective supercell.
initial_defect_structure (Structure) structure corresponding to
initial defect supercell structure (uses Lattice for charge correction)
defect_frac_sc_coords (3 x 1 array) Fractional co-ordinates of
defect location in supercell structure
Returns:
FreysoldtCorrection values as a dictionary
"""
if self.axis is None:
list_axis_grid = np.array(entry.parameters["axis_grid"])
list_bulk_plnr_avg_esp = np.array(entry.parameters["bulk_planar_averages"])
list_defect_plnr_avg_esp = np.array(entry.parameters["defect_planar_averages"])
list_axes = range(len(list_axis_grid))
else:
list_axes = np.array(self.axis)
list_axis_grid, list_bulk_plnr_avg_esp, list_defect_plnr_avg_esp = [], [], []
for ax in list_axes:
list_axis_grid.append(np.array(entry.parameters["axis_grid"][ax]))
list_bulk_plnr_avg_esp.append(np.array(entry.parameters["bulk_planar_averages"][ax]))
list_defect_plnr_avg_esp.append(np.array(entry.parameters["defect_planar_averages"][ax]))
lattice = entry.parameters["initial_defect_structure"].lattice.copy()
defect_frac_coords = entry.parameters["defect_frac_sc_coords"]
q = entry.defect.charge
es_corr = self.perform_es_corr(lattice, entry.charge)
pot_corr_tracker = []
for x, pureavg, defavg, axis in zip(list_axis_grid, list_bulk_plnr_avg_esp, list_defect_plnr_avg_esp,
list_axes):
tmp_pot_corr = self.perform_pot_corr(
x, pureavg, defavg, lattice, entry.charge, defect_frac_coords,
axis, widthsample=1.0)
pot_corr_tracker.append(tmp_pot_corr)
pot_corr = np.mean(pot_corr_tracker)
entry.parameters["freysoldt_meta"] = dict(self.metadata)
entry.parameters["potalign"] = pot_corr / (-q) if q else 0.
return {"freysoldt_electrostatic": es_corr, "freysoldt_potential_alignment": pot_corr}
[docs] def perform_es_corr(self, lattice, q, step=1e-4):
"""
Peform Electrostatic Freysoldt Correction
Args:
lattice: Pymatgen lattice object
q (int): Charge of defect
step (float): step size for numerical integration
Return:
Electrostatic Point Charge contribution to Freysoldt Correction (float)
"""
logger.info("Running Freysoldt 2011 PC calculation (should be " "equivalent to sxdefectalign)")
logger.debug("defect lattice constants are (in angstroms)" + str(lattice.abc))
[a1, a2, a3] = ang_to_bohr * np.array(lattice.get_cartesian_coords(1))
logging.debug("In atomic units, lat consts are (in bohr):" + str([a1, a2, a3]))
vol = np.dot(a1, np.cross(a2, a3)) # vol in bohr^3
def e_iso(encut):
gcut = eV_to_k(encut) # gcut is in units of 1/A
return scipy.integrate.quad(lambda g: self.q_model.rho_rec(g * g) ** 2, step, gcut)[0] * (q ** 2) / np.pi
def e_per(encut):
eper = 0
for g2 in generate_reciprocal_vectors_squared(a1, a2, a3, encut):
eper += (self.q_model.rho_rec(g2) ** 2) / g2
eper *= (q ** 2) * 2 * round(np.pi, 6) / vol
eper += (q ** 2) * 4 * round(np.pi, 6) * self.q_model.rho_rec_limit0 / vol
return eper
eiso = converge(e_iso, 5, self.madetol, self.energy_cutoff)
logger.debug("Eisolated : %f", round(eiso, 5))
eper = converge(e_per, 5, self.madetol, self.energy_cutoff)
logger.info("Eperiodic : %f hartree", round(eper, 5))
logger.info("difference (periodic-iso) is %f hartree", round(eper - eiso, 6))
logger.info("difference in (eV) is %f", round((eper - eiso) * hart_to_ev, 4))
es_corr = round((eiso - eper) / self.dielectric * hart_to_ev, 6)
logger.info("Defect Correction without alignment %f (eV): ", es_corr)
return es_corr
[docs] def perform_pot_corr(self,
axis_grid,
pureavg,
defavg,
lattice,
q,
defect_frac_position,
axis,
widthsample=1.0):
"""
For performing planar averaging potential alignment
Args:
axis_grid (1 x NGX where NGX is the length of the NGX grid
in the axis direction. Same length as pureavg list):
A numpy array which contain the cartesian axis
values (in angstroms) that correspond to each planar avg
potential supplied.
pureavg (1 x NGX where NGX is the length of the NGX grid in
the axis direction.):
A numpy array for the planar averaged
electrostatic potential of the bulk supercell.
defavg (1 x NGX where NGX is the length of the NGX grid in
the axis direction.):
A numpy array for the planar averaged
electrostatic potential of the defect supercell.
lattice: Pymatgen Lattice object of the defect supercell
q (float or int): charge of the defect
defect_frac_position: Fracitional Coordinates of the defect in the supercell
axis (int): axis for performing the freysoldt correction on
widthsample (float): width (in Angstroms) of the region in between defects
where the potential alignment correction is averaged. Default is 1 Angstrom.
Returns:
Potential Alignment contribution to Freysoldt Correction (float)
"""
logging.debug("run Freysoldt potential alignment method for axis " + str(axis))
nx = len(axis_grid)
# shift these planar averages to have defect at origin
axfracval = defect_frac_position[axis]
axbulkval = axfracval * lattice.abc[axis]
if axbulkval < 0:
axbulkval += lattice.abc[axis]
elif axbulkval > lattice.abc[axis]:
axbulkval -= lattice.abc[axis]
if axbulkval:
for i in range(nx):
if axbulkval < axis_grid[i]:
break
rollind = len(axis_grid) - i
pureavg = np.roll(pureavg, rollind)
defavg = np.roll(defavg, rollind)
# if not self._silence:
logger.debug("calculating lr part along planar avg axis")
reci_latt = lattice.reciprocal_lattice
dg = reci_latt.abc[axis]
dg /= ang_to_bohr # convert to bohr to do calculation in atomic units
# Build background charge potential with defect at origin
v_G = np.empty(len(axis_grid), np.dtype("c16"))
v_G[0] = 4 * np.pi * -q / self.dielectric * self.q_model.rho_rec_limit0
g = np.roll(np.arange(-nx / 2, nx / 2, 1, dtype=int), int(nx / 2)) * dg
g2 = np.multiply(g, g)[1:]
v_G[1:] = 4 * np.pi / (self.dielectric * g2) * -q * self.q_model.rho_rec(g2)
v_G[nx // 2] = 0 if not (nx % 2) else v_G[nx // 2]
# Get the real space potential by peforming a fft and grabbing the imaginary portion
v_R = np.fft.fft(v_G)
if abs(np.imag(v_R).max()) > self.madetol:
raise Exception("imaginary part found to be %s", repr(np.imag(v_R).max()))
v_R /= (lattice.volume * ang_to_bohr ** 3)
v_R = np.real(v_R) * hart_to_ev
# get correction
short = (np.array(defavg) - np.array(pureavg) - np.array(v_R))
checkdis = int((widthsample / 2) / (axis_grid[1] - axis_grid[0]))
mid = int(len(short) / 2)
tmppot = [short[i] for i in range(mid - checkdis, mid + checkdis + 1)]
logger.debug("shifted defect position on axis (%s) to origin", repr(axbulkval))
logger.debug("means sampling region is (%f,%f)", axis_grid[mid - checkdis], axis_grid[mid + checkdis])
C = -np.mean(tmppot)
logger.debug("C = %f", C)
final_shift = [short[j] + C for j in range(len(v_R))]
v_R = [elmnt - C for elmnt in v_R]
logger.info("C value is averaged to be %f eV ", C)
logger.info("Potentital alignment energy correction (-q*delta V): %f (eV)", -q * C)
self.pot_corr = -q * C
# log plotting data:
self.metadata["pot_plot_data"][axis] = {
"Vr": v_R,
"x": axis_grid,
"dft_diff": np.array(defavg) - np.array(pureavg),
"final_shift": final_shift,
"check": [mid - checkdis, mid + checkdis + 1]
}
# log uncertainty:
self.metadata["pot_corr_uncertainty_md"][axis] = {"stats": stats.describe(tmppot)._asdict(), "potcorr": -q * C}
return self.pot_corr
[docs] def plot(self, axis, title=None, saved=False):
"""
Plots the planar average electrostatic potential against the Long range and
short range models from Freysoldt. Must run perform_pot_corr or get_correction
(to load metadata) before this can be used.
Args:
axis (int): axis to plot
title (str): Title to be given to plot. Default is no title.
saved (bool): Whether to save file or not. If False then returns plot
object. If True then saves plot as str(title) + "FreyplnravgPlot.pdf"
"""
if not self.metadata["pot_plot_data"]:
raise ValueError("Cannot plot potential alignment before running correction!")
x = self.metadata['pot_plot_data'][axis]['x']
v_R = self.metadata['pot_plot_data'][axis]['Vr']
dft_diff = self.metadata['pot_plot_data'][axis]['dft_diff']
final_shift = self.metadata['pot_plot_data'][axis]['final_shift']
check = self.metadata['pot_plot_data'][axis]['check']
plt.figure()
plt.clf()
plt.plot(x, v_R, c="green", zorder=1, label="long range from model")
plt.plot(x, dft_diff, c="red", label="DFT locpot diff")
plt.plot(x, final_shift, c="blue", label="short range (aligned)")
tmpx = [x[i] for i in range(check[0], check[1])]
plt.fill_between(tmpx, -100, 100, facecolor="red", alpha=0.15, label="sampling region")
plt.xlim(round(x[0]), round(x[-1]))
ymin = min(min(v_R), min(dft_diff), min(final_shift))
ymax = max(max(v_R), max(dft_diff), max(final_shift))
plt.ylim(-0.2 + ymin, 0.2 + ymax)
plt.xlabel(r"distance along axis ($\AA$)", fontsize=15)
plt.ylabel("Potential (V)", fontsize=15)
plt.legend(loc=9)
plt.axhline(y=0, linewidth=0.2, color="black")
plt.title(str(title) + " defect potential", fontsize=18)
plt.xlim(0, max(x))
if saved:
plt.savefig(str(title) + "FreyplnravgPlot.pdf")
return
else:
return plt
[docs]class KumagaiCorrection(DefectCorrection):
"""
A class for KumagaiCorrection class. Largely adapated from PyCDT code
If this correction is used, please reference Kumagai and Oba's original paper
(doi: 10.1103/PhysRevB.89.195205) as well as Freysoldt's original
paper (doi: 10.1103/PhysRevLett.102.016402)
NOTE that equations 8 and 9 from Kumagai et al. reference are divided by (4 pi) to get SI units
"""
def __init__(self,
dielectric_tensor,
sampling_radius=None,
gamma=None):
"""
Initializes the Kumagai Correction
Args:
dielectric_tensor (float or 3x3 matrix): Dielectric constant for the structure
optional data that can be tuned:
sampling_radius (float): radius (in Angstrom) which sites must be outside
of to be included in the correction. Publication by Kumagai advises to
use Wigner-Seitz radius of defect supercell, so this is default value.
gamma (float): convergence parameter for gamma function.
Code will automatically determine this if set to None.
"""
self.metadata = {"gamma": gamma, "sampling_radius": sampling_radius, "potalign": None}
if isinstance(dielectric_tensor, int) or \
isinstance(dielectric_tensor, float):
self.dielectric = np.identity(3) * dielectric_tensor
else:
self.dielectric = np.array(dielectric_tensor)
[docs] def get_correction(self, entry):
"""
Gets the Kumagai correction for a defect entry
Args:
entry (DefectEntry): defect entry to compute Kumagai correction on.
Requires following parameters in the DefectEntry to exist:
bulk_atomic_site_averages (list): list of bulk structure"s atomic site averaged ESPs * charge,
in same order as indices of bulk structure
note this is list given by VASP's OUTCAR (so it is multiplied by a test charge of -1)
defect_atomic_site_averages (list): list of defect structure"s atomic site averaged ESPs * charge,
in same order as indices of defect structure
note this is list given by VASP's OUTCAR (so it is multiplied by a test charge of -1)
site_matching_indices (list): list of corresponding site index values for
bulk and defect site structures EXCLUDING the defect site itself
(ex. [[bulk structure site index, defect structure"s corresponding site index], ... ]
initial_defect_structure (Structure): Pymatgen Structure object representing un-relaxed defect
structure
defect_frac_sc_coords (array): Defect Position in fractional coordinates of the supercell
given in bulk_structure
Returns:
KumagaiCorrection values as a dictionary
"""
bulk_atomic_site_averages = entry.parameters["bulk_atomic_site_averages"]
defect_atomic_site_averages = entry.parameters["defect_atomic_site_averages"]
site_matching_indices = entry.parameters["site_matching_indices"]
defect_sc_structure = entry.parameters["initial_defect_structure"]
defect_frac_sc_coords = entry.parameters["defect_frac_sc_coords"]
lattice = defect_sc_structure.lattice
volume = lattice.volume
q = entry.defect.charge
if not self.metadata["gamma"]:
self.metadata["gamma"] = tune_for_gamma(lattice, self.dielectric)
prec_set = [25, 28]
g_vecs, recip_summation, r_vecs, real_summation = generate_R_and_G_vecs(self.metadata["gamma"],
prec_set, lattice, self.dielectric)
pot_shift = self.get_potential_shift(self.metadata["gamma"], volume)
si = self.get_self_interaction(self.metadata["gamma"])
es_corr = [(real_summation[ind] + recip_summation[ind] + pot_shift + si) for ind in range(2)]
# increase precision if correction is not converged yet
# TODO: allow for larger prec_set to be tried if this fails
if abs(es_corr[0] - es_corr[1]) > 0.0001:
logger.debug("Es_corr summation not converged! ({} vs. {})\nTrying a larger prec_set...".format(es_corr[0],
es_corr[1]))
prec_set = [30, 35]
g_vecs, recip_summation, r_vecs, real_summation = generate_R_and_G_vecs(self.metadata["gamma"],
prec_set, lattice, self.dielectric)
es_corr = [(real_summation[ind] + recip_summation[ind] + pot_shift + si) for ind in range(2)]
if abs(es_corr[0] - es_corr[1]) < 0.0001:
raise ValueError("Correction still not converged after trying prec_sets up to 35... serious error.")
es_corr = es_corr[0] * -(q ** 2.) * kumagai_to_V / 2. # [eV]
# if no sampling radius specified for pot align, then assuming Wigner-Seitz radius:
if not self.metadata["sampling_radius"]:
wz = lattice.get_wigner_seitz_cell()
dist = []
for facet in wz:
midpt = np.mean(np.array(facet), axis=0)
dist.append(np.linalg.norm(midpt))
self.metadata["sampling_radius"] = min(dist)
# assemble site_list based on matching indices
# [[defect_site object, Vqb for site], .. repeat for all non defective sites]
site_list = []
for bs_ind, ds_ind in site_matching_indices:
Vqb = -(defect_atomic_site_averages[int(ds_ind)] - bulk_atomic_site_averages[int(bs_ind)])
site_list.append([defect_sc_structure[int(ds_ind)], Vqb])
pot_corr = self.perform_pot_corr(defect_sc_structure, defect_frac_sc_coords, site_list,
self.metadata["sampling_radius"], q, r_vecs[0],
g_vecs[0], self.metadata["gamma"])
entry.parameters["kumagai_meta"] = dict(self.metadata)
entry.parameters["potalign"] = pot_corr / (-q) if q else 0.
return {"kumagai_electrostatic": es_corr, "kumagai_potential_alignment": pot_corr}
[docs] def perform_es_corr(self, gamma, prec, lattice, charge):
"""
Peform Electrostatic Kumagai Correction
Args:
gamma (float): Ewald parameter
prec (int): Precision parameter for reciprical/real lattice vector generation
lattice: Pymatgen Lattice object corresponding to defect supercell
charge (int): Defect charge
Return:
Electrostatic Point Charge contribution to Kumagai Correction (float)
"""
volume = lattice.volume
g_vecs, recip_summation, r_vecs, real_summation = generate_R_and_G_vecs(gamma, [prec], lattice, self.dielectric)
recip_summation = recip_summation[0]
real_summation = real_summation[0]
es_corr = (recip_summation + real_summation + self.get_potential_shift(gamma, volume) +
self.get_self_interaction(gamma))
es_corr *= -(charge ** 2.) * kumagai_to_V / 2. # [eV]
return es_corr
[docs] def perform_pot_corr(self, defect_structure, defect_frac_coords, site_list,
sampling_radius, q, r_vecs, g_vecs, gamma):
"""
For performing potential alignment in manner described by Kumagai et al.
Args:
defect_structure: Pymatgen Structure object corrsponding to the defect supercell
defect_frac_coords (array): Defect Position in fractional coordinates of the supercell
given in bulk_structure
site_list: list of corresponding site index values for
bulk and defect site structures EXCLUDING the defect site itself
(ex. [[bulk structure site index, defect structure"s corresponding site index], ... ]
sampling_radius (float): radius (in Angstrom) which sites must be outside
of to be included in the correction. Publication by Kumagai advises to
use Wigner-Seitz radius of defect supercell, so this is default value.
q (int): Defect charge
r_vecs: List of real lattice vectors to use in summation
g_vecs: List of reciprocal lattice vectors to use in summation
gamma (float): Ewald parameter
Return:
Potential alignment contribution to Kumagai Correction (float)
"""
volume = defect_structure.lattice.volume
potential_shift = self.get_potential_shift(gamma, volume)
pot_dict = {} # keys will be site index in the defect structure
for_correction = [] # region to sample for correction
# for each atom, do the following:
# (a) get relative_vector from defect_site to site in defect_supercell structure
# (b) recalculate the recip and real summation values based on this r_vec
# (c) get information needed for pot align
for site, Vqb in site_list:
dist, jimage = site.distance_and_image_from_frac_coords(defect_frac_coords)
vec_defect_to_site = defect_structure.lattice.get_cartesian_coords(site.frac_coords -
jimage - defect_frac_coords)
dist_to_defect = np.linalg.norm(vec_defect_to_site)
if abs(dist_to_defect - dist) > 0.001:
raise ValueError("Error in computing vector to defect")
relative_real_vectors = [r_vec - vec_defect_to_site for r_vec in r_vecs[:]]
real_sum = self.get_real_summation(gamma, relative_real_vectors)
recip_sum = self.get_recip_summation(gamma, g_vecs, volume, r=vec_defect_to_site[:])
Vpc = (real_sum + recip_sum + potential_shift) * kumagai_to_V * q
defect_struct_index = defect_structure.index(site)
pot_dict[defect_struct_index] = {
"Vpc": Vpc,
"Vqb": Vqb,
"dist_to_defect": dist_to_defect
}
logger.debug("For atom {}\n\tbulk/defect DFT potential difference = "
"{}".format(defect_struct_index, Vqb))
logger.debug("\tanisotropic model charge: {}".format(Vpc))
logger.debug("\t\treciprocal part: {}".format(recip_sum * kumagai_to_V * q))
logger.debug("\t\treal part: {}".format(real_sum * kumagai_to_V * q))
logger.debug("\t\tself interaction part: {}".format(potential_shift * kumagai_to_V * q))
logger.debug("\trelative_vector to defect: {}".format(vec_defect_to_site))
if dist_to_defect > sampling_radius:
logger.debug("\tdistance to defect is {} which is outside minimum sampling "
"radius {}".format(dist_to_defect,
sampling_radius))
for_correction.append(Vqb - Vpc)
else:
logger.debug("\tdistance to defect is {} which is inside minimum sampling "
"radius {} (so will not include for correction)"
"".format(dist_to_defect, sampling_radius))
if len(for_correction):
pot_alignment = np.mean(for_correction)
else:
logger.info("No atoms sampled for_correction radius!"
" Assigning potential alignment value of 0.")
pot_alignment = 0.
self.metadata["potalign"] = pot_alignment
pot_corr = -q * pot_alignment
# log uncertainty stats:
self.metadata["pot_corr_uncertainty_md"] = {
"stats": stats.describe(for_correction)._asdict(),
"number_sampled": len(for_correction)
}
self.metadata["pot_plot_data"] = pot_dict
logger.info("Kumagai potential alignment (site averaging): %f", pot_alignment)
logger.info("Kumagai potential alignment correction energy: %f eV", pot_corr)
return pot_corr
[docs] def get_real_summation(self, gamma, real_vectors):
"""
Get real summation term from list of real-space vectors
"""
real_part = 0
invepsilon = np.linalg.inv(self.dielectric)
rd_epsilon = np.sqrt(np.linalg.det(self.dielectric))
for r_vec in real_vectors:
if np.linalg.norm(r_vec) > 1e-8:
loc_res = np.sqrt(np.dot(r_vec, np.dot(invepsilon, r_vec)))
nmr = scipy.special.erfc(gamma * loc_res)
real_part += nmr / loc_res
real_part /= (4 * np.pi * rd_epsilon)
return real_part
[docs] def get_recip_summation(self, gamma, recip_vectors, volume, r=[0., 0., 0.]):
"""
Get Reciprocal summation term from list of reciprocal-space vectors
"""
recip_part = 0
for g_vec in recip_vectors:
# dont need to avoid G=0, because it will not be
# in recip list (if generate_R_and_G_vecs is used)
Gdotdiel = np.dot(g_vec, np.dot(self.dielectric, g_vec))
summand = np.exp(-Gdotdiel / (4 * (gamma ** 2))) * np.cos(np.dot(g_vec, r)) / Gdotdiel
recip_part += summand
recip_part /= volume
return recip_part
[docs] def get_self_interaction(self, gamma):
"""
Args:
gamma ():
Returns:
Self-interaction energy of defect.
"""
determ = np.linalg.det(self.dielectric)
return - gamma / (2. * np.pi * np.sqrt(np.pi * determ))
[docs] def get_potential_shift(self, gamma, volume):
"""
Args:
gamma (float): Gamma
volume (float): Volume.
Returns:
Potential shift for defect.
"""
return - 0.25 / (volume * gamma ** 2.)
[docs] def plot(self, title=None, saved=False):
"""
Plots the AtomicSite electrostatic potential against the Long range and short range models
from Kumagai and Oba (doi: 10.1103/PhysRevB.89.195205)
"""
if "pot_plot_data" not in self.metadata.keys():
raise ValueError("Cannot plot potential alignment before running correction!")
sampling_radius = self.metadata["sampling_radius"]
site_dict = self.metadata["pot_plot_data"]
potalign = self.metadata["potalign"]
plt.figure()
plt.clf()
distances, sample_region = [], []
Vqb_list, Vpc_list, diff_list = [], [], []
for site_ind, site_dict in site_dict.items():
dist = site_dict["dist_to_defect"]
distances.append(dist)
Vqb = site_dict["Vqb"]
Vpc = site_dict["Vpc"]
Vqb_list.append(Vqb)
Vpc_list.append(Vpc)
diff_list.append(Vqb - Vpc)
if dist > sampling_radius:
sample_region.append(Vqb - Vpc)
plt.plot(distances, Vqb_list,
color='r', marker='^', linestyle='None',
label='$V_{q/b}$')
plt.plot(distances, Vpc_list,
color='g', marker='o', linestyle='None',
label='$V_{pc}$')
plt.plot(distances, diff_list, color='b', marker='x', linestyle='None',
label='$V_{q/b}$ - $V_{pc}$')
x = np.arange(sampling_radius, max(distances) * 1.05, 0.01)
y_max = max(max(Vqb_list), max(Vpc_list), max(diff_list)) + .1
y_min = min(min(Vqb_list), min(Vpc_list), min(diff_list)) - .1
plt.fill_between(x, y_min, y_max, facecolor='red',
alpha=0.15, label='sampling region')
plt.axhline(y=potalign, linewidth=0.5, color='red',
label='pot. align. / -q')
plt.legend(loc=0)
plt.axhline(y=0, linewidth=0.2, color='black')
plt.ylim([y_min, y_max])
plt.xlim([0, max(distances) * 1.1])
plt.xlabel(r'Distance from defect ($\AA$)', fontsize=20)
plt.ylabel('Potential (V)', fontsize=20)
plt.title(str(title) + " atomic site potential plot", fontsize=20)
if saved:
plt.savefig(str(title) + "KumagaiESPavgPlot.pdf")
else:
return plt
[docs]class BandFillingCorrection(DefectCorrection):
"""
A class for BandFillingCorrection class. Largely adapted from PyCDT code
"""
def __init__(self, resolution=0.01):
"""
Initializes the Bandfilling correction
Args:
resolution (float): energy resolution to maintain for gap states
"""
self.resolution = resolution
self.metadata = {
"num_hole_vbm": None,
"num_elec_cbm": None,
"potalign": None
}
[docs] def get_correction(self, entry):
"""
Gets the BandFilling correction for a defect entry
Args:
entry (DefectEntry): defect entry to compute BandFilling correction on.
Requires following parameters in the DefectEntry to exist:
eigenvalues
dictionary of defect eigenvalues, as stored in a Vasprun object
kpoint_weights (list of floats)
kpoint weights corresponding to the dictionary of eigenvalues,
as stored in a Vasprun object
potalign (float)
potential alignment for the defect calculation
Only applies to non-zero charge,
When using potential alignment correction (freysoldt or kumagai),
need to divide by -q
cbm (float)
CBM of bulk calculation (or band structure calculation of bulk);
calculated on same level of theory as the defect
(ex. GGA defects -> requires GGA cbm)
vbm (float)
VBM of bulk calculation (or band structure calculation of bulk);
calculated on same level of theory as the defect
(ex. GGA defects -> requires GGA vbm)
Returns:
Bandfilling Correction value as a dictionary
"""
eigenvalues = entry.parameters["eigenvalues"]
kpoint_weights = entry.parameters["kpoint_weights"]
potalign = entry.parameters["potalign"]
vbm = entry.parameters["vbm"]
cbm = entry.parameters["cbm"]
bf_corr = self.perform_bandfill_corr(eigenvalues, kpoint_weights, potalign, vbm, cbm)
entry.parameters["bandfilling_meta"] = dict(self.metadata)
return {"bandfilling_correction": bf_corr}
[docs] def perform_bandfill_corr(self, eigenvalues, kpoint_weights, potalign, vbm, cbm):
"""
This calculates the band filling correction based on excess of electrons/holes in CB/VB...
Note that the total free holes and electrons may also be used for a "shallow donor/acceptor"
correction with specified band shifts:
+num_elec_cbm * Delta E_CBM (or -num_hole_vbm * Delta E_VBM)
"""
bf_corr = 0.
self.metadata["potalign"] = potalign
self.metadata["num_hole_vbm"] = 0.
self.metadata["num_elec_cbm"] = 0.
core_occupation_value = list(eigenvalues.values())[0][0][0][1] # get occupation of a core eigenvalue
if len(eigenvalues.keys()) == 1:
# needed because occupation of non-spin calcs is sometimes still 1... should be 2
spinfctr = 2. if core_occupation_value == 1. else 1.
elif len(eigenvalues.keys()) == 2:
spinfctr = 1.
else:
raise ValueError("Eigenvalue keys greater than 2")
# for tracking mid gap states...
shifted_cbm = potalign + cbm # shift cbm with potential alignment
shifted_vbm = potalign + vbm # shift vbm with potential alignment
for spinset in eigenvalues.values():
for kptset, weight in zip(spinset, kpoint_weights):
for eig, occu in kptset: # eig is eigenvalue and occu is occupation
if (occu and (eig > shifted_cbm - self.resolution)): # donor MB correction
bf_corr += weight * spinfctr * occu * (eig - shifted_cbm) # "move the electrons down"
self.metadata["num_elec_cbm"] += weight * spinfctr * occu
elif (occu != core_occupation_value) and (
eig <= shifted_vbm + self.resolution): # acceptor MB correction
bf_corr += weight * spinfctr * (core_occupation_value - occu) * (
shifted_vbm - eig) # "move the holes up"
self.metadata["num_hole_vbm"] += weight * spinfctr * (core_occupation_value - occu)
bf_corr *= -1 # need to take negative of this shift for energetic correction
return bf_corr
[docs]class BandEdgeShiftingCorrection(DefectCorrection):
"""
A class for BandEdgeShiftingCorrection class. Largely adapted from PyCDT code
"""
def __init__(self):
"""
Initializes the BandEdgeShiftingCorrection class
"""
self.metadata = {
"vbmshift": 0.,
"cbmshift": 0.,
}
[docs] def get_correction(self, entry):
"""
Gets the BandEdge correction for a defect entry
Args:
entry (DefectEntry): defect entry to compute BandFilling correction on.
Requires some parameters in the DefectEntry to properly function:
hybrid_cbm (float)
CBM of HYBRID bulk calculation one wishes to shift to
hybrid_vbm (float)
VBM of HYBRID bulk calculation one wishes to shift to
cbm (float)
CBM of bulk calculation (or band structure calculation of bulk);
calculated on same level of theory as the defect
(ex. GGA defects -> requires GGA cbm)
vbm (float)
VBM of bulk calculation (or band structure calculation of bulk);
calculated on same level of theory as the defect
(ex. GGA defects -> requires GGA vbm)
Returns:
BandfillingCorrection value as a dictionary
"""
hybrid_cbm = entry.parameters["hybrid_cbm"]
hybrid_vbm = entry.parameters["hybrid_vbm"]
vbm = entry.parameters["vbm"]
cbm = entry.parameters["cbm"]
self.metadata["vbmshift"] = hybrid_vbm - vbm # note vbmshift has UPWARD as positive convention
self.metadata["cbmshift"] = hybrid_cbm - cbm # note cbmshift has UPWARD as positive convention
charge = entry.charge
bandedgeshifting_correction = charge * self.metadata["vbmshift"]
entry.parameters["bandshift_meta"] = dict(self.metadata)
return {
"bandedgeshifting_correction": bandedgeshifting_correction
}