Merge pull request #83 from materialsproject/kumaga

Kumagai EFNV correction

.github/workflows/testing.yml

@@ -82,6 +82,7 @@ jobs:
       - name: Install dependencies
         run: |
           pip install -e .[strict]
+          pip install -e .[docs]
       - name: Build
         run: jupyter-book build --path-output docs docs/source
         # Develop branch only

.pre-commit-config.yaml

@@ -72,7 +72,7 @@ repos:
     language: python
     types: [text]
     args: [
-      --ignore-words-list, 'titel,statics,ba,nd,te,mater,commun',
+      --ignore-words-list, 'titel,statics,ba,nd,te,mater,commun,vise',
       --skip, "*.ipynb"
       ]
 - repo: https://github.com/kynan/nbstripout

pymatgen/analysis/defects/corrections/freysoldt.py

@@ -3,7 +3,6 @@
 from __future__ import annotations
 
 import logging
-from collections import namedtuple
 from typing import Optional
 
 import matplotlib.pyplot as plt
@@ -15,6 +14,7 @@
 from scipy import stats
 
 from pymatgen.analysis.defects.utils import (
+    CorrectionResult,
     QModel,
     ang_to_bohr,
     converge,
@@ -36,16 +36,6 @@
 Rewritten to be functional instead of object oriented.
 """
 
-"""
-Named tuple for storing result of correction result.
-
-Metadata contains plotting data for the planar average electrostatic potential.
-key 0, 1, 2 correspond to the x, y, z axes respectively.
-"""
-FreysoldtSummary = namedtuple(
-    "FreysoldtSummary", ["electrostatic", "potential_alignment", "metadata"]
-)
-
 
 def get_freysoldt_correction(
     q: int,
@@ -58,17 +48,9 @@ def get_freysoldt_correction(
     mad_tol: float = 1e-4,
     q_model: Optional[QModel] = None,
     step: float = 1e-4,
-) -> FreysoldtSummary:
+) -> CorrectionResult:
     """Gets the Freysoldt correction for a defect entry.
 
-    Get the Freysoldt correction for a defect. The result is given
-    as a dictionary with the following keys:
-    - "freysoldt_electrostatic": the electrostatic portion of the correction
-    - "freysoldt_potential_alignment": the potential alignment portion of the correction
-
-    The second return value is a dictionary used to plot the planar average
-    electrostatic potential (used for debugging).
-
     Args:
         q:
             Charge state of defect
@@ -90,8 +72,11 @@ def get_freysoldt_correction(
             Step size for numerical integration.
 
     Returns:
-        dict: Freysoldt correction values as a dictionary
-        dict: Plot data for planar average electrostatic potential.
+        CorrectionResult: Correction summary object. The metadata contains
+            plotting data for the planar average electrostatic potential.
+            ```
+            plot_plnr_avg(result.metadata[0], title="Lattice Direction 1")
+            ```
     """
     # dielectric has to be a float
     if isinstance(dielectric, (int, float)):
@@ -163,10 +148,9 @@ def get_freysoldt_correction(
         plot_data[axis] = md
 
     pot_corr = np.mean(list(pot_corrs.values()))
-
-    return FreysoldtSummary(
-        electrostatic=es_corr,
-        potential_alignment=pot_corr / (-q) if q else 0,
+    pot_align = pot_corr / (-q) if q else 0
+    return CorrectionResult(
+        correction_energy=es_corr + pot_align,
         metadata=plot_data,
     )
 

pymatgen/analysis/defects/corrections/kumagai.py

@@ -0,0 +1,112 @@
+"""Extended FNV correction by Kumagai/Oba.
+
+Works as a light wrapper around pydefect.
+https://github.com/kumagai-group/pydefect
+"""
+
+from __future__ import annotations
+
+import logging
+import math
+from pathlib import Path
+
+from pymatgen.analysis.defects.utils import CorrectionResult, get_zfile
+
+# check that pydefect is installed
+try:
+    pass
+except ImportError:
+    raise ImportError("pydefect is not installed. Please install it first.")
+
+# Disable messages from pydefect import
+from vise import user_settings
+
+user_settings.logger.setLevel(logging.CRITICAL)
+
+from pydefect.analyzer.calc_results import CalcResults
+from pydefect.cli.vasp.make_efnv_correction import make_efnv_correction
+from pymatgen.core import Structure
+from pymatgen.io.vasp import Outcar, Vasprun
+
+__author__ = "Jimmy-Xuan Shen"
+__copyright__ = "Copyright 2022, The Materials Project"
+__maintainer__ = "Jimmy-Xuan Shen"
+__email__ = "[email protected]"
+
+_logger = logging.getLogger(__name__)
+# suppress pydefect INFO messages
+logging.getLogger("pydefect").setLevel(logging.WARNING)
+
+
+def get_structure_with_pot(directory: Path) -> Structure:
+    """Reads vasprun.xml and OUTCAR files in a directory.
+
+    Args:
+        directory: Directory containing vasprun.xml and OUTCAR files.
+
+    Returns:
+        Structure with "potential" site property.
+    """
+    d_ = Path(directory)
+    f_vasprun = get_zfile(d_, "vasprun.xml")
+    f_outcar = get_zfile(d_, "OUTCAR")
+    vasprun = Vasprun(f_vasprun)
+    outcar = Outcar(f_outcar)
+
+    calc = CalcResults(
+        structure=vasprun.final_structure,
+        energy=outcar.final_energy,
+        magnetization=outcar.total_mag or 0.0,
+        potentials=[-p for p in outcar.electrostatic_potential],
+        electronic_conv=vasprun.converged_electronic,
+        ionic_conv=vasprun.converged_ionic,
+    )
+
+    struct = calc.structure.copy(site_properties={"potential": calc.potentials})
+    return struct
+
+
+def get_efnv_correction(
+    charge: int,
+    defect_structure: Structure,
+    bulk_structure: Structure,
+    dielectric_tensor: list[list[float]],
+    **kwargs,
+) -> CorrectionResult:
+    """Returns the Kumagai/Oba EFNV correction for a given defect.
+
+    Args:
+        charge: Charge of the defect.
+        defect_structure: Defect structure.
+        bulk_structure: Bulk structure.
+        dielectric_tensor: Dielectric tensor.
+        **kwargs: Keyword arguments to pass to `make_efnv_correction`.
+    """
+    # ensure that the structures have the "potential" site property
+    bulk_potentials = [site.properties["potential"] for site in defect_structure]
+    defect_potentials = [site.properties["potential"] for site in bulk_structure]
+
+    defect_calc_results = CalcResults(
+        structure=defect_structure,
+        energy=math.inf,
+        magnetization=math.inf,
+        potentials=defect_potentials,
+    )
+    bulk_calc_results = CalcResults(
+        structure=bulk_structure,
+        energy=math.inf,
+        magnetization=math.inf,
+        potentials=bulk_potentials,
+    )
+
+    efnv_corr = make_efnv_correction(
+        charge=charge,
+        calc_results=defect_calc_results,
+        perfect_calc_results=bulk_calc_results,
+        dielectric_tensor=dielectric_tensor,
+        **kwargs,
+    )
+
+    return CorrectionResult(
+        correction_energy=efnv_corr.correction_energy, metadata={"efnv_corr": efnv_corr}
+    )

pymatgen/analysis/defects/thermo.py

@@ -21,12 +21,9 @@
 from scipy.spatial import ConvexHull
 
 from pymatgen.analysis.defects.core import Defect
-from pymatgen.analysis.defects.corrections.freysoldt import (
-    FreysoldtSummary,
-    get_freysoldt_correction,
-)
+from pymatgen.analysis.defects.corrections.freysoldt import get_freysoldt_correction
 from pymatgen.analysis.defects.finder import DefectSiteFinder
-from pymatgen.analysis.defects.utils import get_zfile
+from pymatgen.analysis.defects.utils import CorrectionResult, get_zfile
 
 __author__ = "Jimmy-Xuan Shen, Danny Broberg, Shyam Dwaraknath"
 __copyright__ = "Copyright 2022, The Materials Project"
@@ -72,7 +69,7 @@ def __post_init__(self) -> None:
         """Post-initialization."""
         self.charge_state = int(self.charge_state)
         self.corrections: dict = {} if self.corrections is None else self.corrections
-        self.correction_type: str = "Freysoldt"
+        self.correction_type: str = "freysoldt"
         self.correction_metadata: dict = (
             {} if self.correction_metadata is None else self.correction_metadata
         )
@@ -85,7 +82,7 @@ def get_freysoldt_correction(
         defect_struct: Optional[Structure] = None,
         bulk_struct: Optional[Structure] = None,
         **kwargs,
-    ) -> FreysoldtSummary:
+    ) -> CorrectionResult:
         """Calculate the Freysoldt correction.
 
         Updates the corrections dictionary with the Freysoldt correction
@@ -144,18 +141,16 @@ def get_freysoldt_correction(
         )
         self.corrections.update(
             {
-                "electrostatic": frey_corr.electrostatic,
-                "potential_alignment": frey_corr.potential_alignment,
+                "freysoldt": frey_corr.correction_energy,
             }
         )
-        self.correction_metadata.update(frey_corr.metadata.copy())
-        self.correction_type = "Freysoldt"
+        self.correction_metadata.update({"freysoldt": frey_corr.metadata.copy()})
         return frey_corr
 
     @property
     def corrected_energy(self) -> float:
         """The energy of the defect entry with all corrections applied."""
-        return self.sc_entry.energy + sum(self.corrections.values())
+        return self.sc_entry.energy + self.corrections["freysoldt"]
 
 
 @dataclass

pymatgen/analysis/defects/utils.py

@@ -9,6 +9,7 @@
 import operator
 from collections import defaultdict
 from copy import deepcopy
+from dataclasses import dataclass
 from functools import cached_property
 from pathlib import Path
 from typing import Any, Callable, Generator
@@ -969,3 +970,17 @@ def get_symmetry_labeled_structures(
     # Label the groups by structure matching
     site_labels = generic_groupby(inserted_structs, comp=sm.fit)
     return [*zip(sites.tolist(), site_labels)]
+
+
+@dataclass
+class CorrectionResult(MSONable):
+    """A summary of the corrections applied to a structure.
+
+    Attributes:
+        electrostatic: The electrostatic correction.
+        potential_alignment: The potential alignment correction.
+        metadata: A dictionary of metadata for plotting and intermediate analysis.
+    """
+
+    correction_energy: float
+    metadata: dict[Any, Any]

pyproject.toml

@@ -30,15 +30,16 @@ dev = ["pre-commit>=2.12.1"]
 docs = [
   "jupyter-book>=0.13.1",
 ]
+pydefect = ["pydefect>=0.6.2"]
 
 strict = [
   "pymatgen==2022.11.7",
   "dscribe==1.2.2",
   "scikit-image==0.19.3",
   "pytest==7.2.0",
   "pytest-cov==4.0.0",
-  "jupyter-book==0.13.1",
   "pre-commit==2.21.0",
+  "pydefect==0.6.2",
 ]
 
 tests = ["pytest==7.2.0", "pytest-cov==4.0.0"]

tests/test_corrections.py

@@ -4,6 +4,10 @@
     get_freysoldt_correction,
     plot_plnr_avg,
 )
+from pymatgen.analysis.defects.corrections.kumagai import (
+    get_efnv_correction,
+    get_structure_with_pot,
+)
 
 
 def test_freysoldt(data_Mg_Ga):
@@ -17,8 +21,7 @@ def test_freysoldt(data_Mg_Ga):
         bulk_locpot=bulk_locpot,
         defect_frac_coords=[0.5, 0.5, 0.5],
     )
-    assert freysoldt_summary.electrostatic == pytest.approx(0, abs=1e-4)
-    assert freysoldt_summary.potential_alignment == pytest.approx(0, abs=1e-4)
+    assert freysoldt_summary.correction_energy == pytest.approx(0, abs=1e-4)
 
     # simple check that the plotter works
     plot_plnr_avg(freysoldt_summary.metadata[0])
@@ -41,3 +44,20 @@ def test_freysoldt(data_Mg_Ga):
         bulk_locpot=[*map(bulk_locpot.get_axis_grid, [0, 1, 2])],
         defect_frac_coords=[0.5, 0.5, 0.5],
     )
+
+
+def test_kumagai(test_dir):
+
+    sb = get_structure_with_pot(test_dir / "Mg_Ga" / "bulk_sc")
+    sd0 = get_structure_with_pot(test_dir / "Mg_Ga" / "q=0")
+    sd1 = get_structure_with_pot(test_dir / "Mg_Ga" / "q=1")
+
+    res0 = get_efnv_correction(
+        0, sd0, sb, dielectric_tensor=[[1, 0, 0], [0, 1, 0], [0, 0, 1]]
+    )
+    assert res0.correction_energy == pytest.approx(0, abs=1e-4)
+
+    res1 = get_efnv_correction(
+        1, sd1, sb, dielectric_tensor=[[1, 0, 0], [0, 1, 0], [0, 0, 1]]
+    )
+    assert res1.correction_energy > 0