Major Fixes to Selective Dynamics and Oxidation (#127)

* fixed ylim

* limit number of interstitials

* selective dyanmics bug fix

- fixed bug in centering the selective dyanamics sphere
- allow perturbation of the sites

* remove none

* fixed frac/cart

* fixed tests

* modified oxi-state for substitution

* modified oxi-state for substitution

* max insertions

* Update test_core.py

coverage

* [pre-commit.ci] auto fixes from pre-commit.com hooks

for more information, see https://pre-commit.ci

* fixed test

---------

Co-authored-by: pre-commit-ci[bot] <66853113+pre-commit-ci[bot]@users.noreply.github.com>

docs/source/content/freysoldt-correction.ipynb

@@ -1,6 +1,7 @@
 {
  "cells": [
   {
+   "attachments": {},
    "cell_type": "markdown",
    "metadata": {},
    "source": [
@@ -44,6 +45,7 @@
    ]
   },
   {
+   "attachments": {},
    "cell_type": "markdown",
    "metadata": {},
    "source": [
@@ -83,6 +85,7 @@
    ]
   },
   {
+   "attachments": {},
    "cell_type": "markdown",
    "metadata": {},
    "source": [
@@ -189,6 +192,7 @@
    ]
   },
   {
+   "attachments": {},
    "cell_type": "markdown",
    "metadata": {},
    "source": [

pymatgen/analysis/defects/core.py

@@ -164,6 +164,7 @@ def get_supercell_structure(
         min_length: float = 10.0,
         force_diagonal: bool = False,
         relax_radius: float | str | None = None,
+        perturb: float | None = None,
     ) -> Structure:
         """Generate the supercell for a defect.
 
@@ -175,6 +176,8 @@ def get_supercell_structure(
             min_length: Minimum length of the smallest supercell lattice vector.
             force_diagonal: If True, return a transformation with a diagonal transformation matrix.
             relax_radius: Relax the supercell atoms to a sphere of this radius around the defect site.
+            perturb: The amount to perturb the sites in the supercell. Only perturb the sites with
+                selective dynamics set to True. So this setting only works with `relax_radius`.
 
         Returns:
             Structure: The supercell structure.
@@ -206,9 +209,12 @@ def get_supercell_structure(
 
         _set_selective_dynamics(
             structure=sc_defect_struct,
-            site_pos=sc_pos,
+            site_pos=sc_site.coords,
             relax_radius=relax_radius,
         )
+        if perturb is not None:
+            _perturb_dynamic_sites(sc_defect_struct, distance=perturb)
+
         return sc_defect_struct
 
     @property
@@ -356,16 +362,11 @@ def name(self) -> str:
     def defect_structure(self) -> Structure:
         """Returns the defect structure."""
         struct: Structure = self.structure.copy()
-        rm_oxi = struct.sites[self.defect_site_index].specie.oxi_state
         struct.remove_sites([self.defect_site_index])
-        sub_states = self.site.specie.icsd_oxidation_states
-        if len(sub_states) == 0:
-            sub_states = self.site.specie.oxidation_states
-        sub_oxi = min(sub_states, key=lambda x: abs(x - rm_oxi))
-        sub_specie = Species(self.site.specie.symbol, sub_oxi)
+        insert_el = get_element(self.site.specie)
         struct.insert(
             self.defect_site_index,
-            species=sub_specie,
+            species=insert_el,
             coords=np.mod(self.site.frac_coords, 1),
         )
         return struct
@@ -406,9 +407,15 @@ def _guess_oxi_state(self) -> float:
         Returns:
             float: The oxidation state of the defect.
         """
-        orig_site = self.defect_site
-        sub_site = self.defect_structure[self.defect_site_index]
-        return sub_site.specie.oxi_state - orig_site.specie.oxi_state
+        rm_oxi = self.structure[self.defect_site_index].specie.oxi_state
+        sub_states = self.site.specie.common_oxidation_states
+        if len(sub_states) == 0:
+            raise ValueError(
+                f"No common oxidation states found for {self.site.specie}."
+                "Please specify the oxidation state manually."
+            )
+        sub_oxi = sub_states[0]
+        return sub_oxi - rm_oxi
 
     def __repr__(self) -> str:
         """Representation of a substitutional defect."""
@@ -579,6 +586,7 @@ def get_supercell_structure(
         min_length: float = 10.0,
         force_diagonal: bool = False,
         relax_radius: float | str | None = None,
+        perturb: float | None = None,
     ) -> Structure:
         """Generate the supercell for a defect.
 
@@ -590,6 +598,9 @@ def get_supercell_structure(
             min_atoms: Minimum number of atoms allowed in the supercell.
             min_length: Minimum length of the smallest supercell lattice vector.
             force_diagonal: If True, return a transformation with a diagonal transformation matrix.
+            relax_radius: Relax the supercell atoms to a sphere of this radius around the defect site.
+            perturb: The amount to perturb the sites in the supercell. Only perturb the sites with
+                selective dynamics set to True. So this setting only works with `relax_radius`.
 
         Returns:
             Structure: The supercell structure.
@@ -619,9 +630,13 @@ def get_supercell_structure(
 
         _set_selective_dynamics(
             structure=sc_defect_struct,
-            site_pos=sc_pos,
+            site_pos=sc_site.coords,
             relax_radius=relax_radius,
         )
+
+        if perturb is not None:
+            _perturb_dynamic_sites(sc_defect_struct, distance=perturb)
+
         return sc_defect_struct
 
 
@@ -723,10 +738,9 @@ def _set_selective_dynamics(structure, site_pos, relax_radius):
         relax_radius = min(get_plane_spacing(structure.lattice.matrix)) / 2.0
     if not isinstance(relax_radius, float):
         raise ValueError("relax_radius must be a float or 'auto' or None")
-
     structure.get_sites_in_sphere(site_pos, relax_radius)
     relax_sites = structure.get_sites_in_sphere(
-        [0, 0, 0], relax_radius, include_index=True
+        site_pos, relax_radius, include_index=True
     )
     relax_indices = [site.index for site in relax_sites]
     relax_mask = [[False, False, False]] * len(structure)
@@ -735,6 +749,56 @@ def _set_selective_dynamics(structure, site_pos, relax_radius):
     structure.add_site_property("selective_dynamics", relax_mask)
 
 
+def perturb_sites(
+    structure,
+    distance: float,
+    min_distance: float | None = None,
+    site_indices: list | None = None,
+) -> None:
+    """
+    Performs a random perturbation of the sites in a structure to break
+    symmetries.
+
+    Args:
+        structure (Structure): Input structure.
+        distance (float): Distance in angstroms by which to perturb each
+            site.
+        min_distance (None, int, or float): if None, all displacements will
+            be equal amplitude. If int or float, perturb each site a
+            distance drawn from the uniform distribution between
+            'min_distance' and 'distance'.
+        site_indices (list): List of site indices on which to perform the
+            perturbation. If None, all sites will be perturbed.
+
+    """
+
+    def get_rand_vec():
+        # deals with zero vectors.
+        vector = np.random.randn(3)
+        vnorm = np.linalg.norm(vector)
+        dist = distance
+        if isinstance(min_distance, (float, int)):
+            dist = np.random.uniform(min_distance, dist)
+        return vector / vnorm * dist if vnorm != 0 else get_rand_vec()
+
+    if site_indices is None:
+        site_indices_ = list(range(len(structure._sites)))
+    else:
+        site_indices_ = site_indices
+
+    for i in site_indices_:
+        structure.translate_sites([i], get_rand_vec(), frac_coords=False)
+
+
+def _perturb_dynamic_sites(structure, distance):
+    free_indices = [
+        i
+        for i, site in enumerate(structure)
+        if site.properties["selective_dynamics"][0]
+    ]
+    perturb_sites(structure=structure, distance=distance, site_indices=free_indices)
+
+
 # TODO: matching defect complexes might be done with some kind of CoM site to fix the periodicity
 # Get this by taking the periodic average of all the provided sites.
 # class DefectComplex(DummySpecies):

pymatgen/analysis/defects/finder.py

@@ -48,18 +48,26 @@ def __init__(self, symprec: float = 0.01, angle_tolerance: float = 5.0):
         self.angle_tolerance = angle_tolerance
 
     def get_defect_fpos(
-        self, defect_structure: Structure, base_structure: Structure
+        self,
+        defect_structure: Structure,
+        base_structure: Structure,
+        remove_oxi: bool = True,
     ) -> ArrayLike:
         """Get the position of a defect in the pristine structure.
 
         Args:
             defect_structure: Relaxed structure containing the defect
             base_structure: Structure for the pristine cell
+            remove_oxi: Whether to remove oxidation states from the structures
 
         Returns:
             ArrayLike: Position of the defect in the pristine structure
             (in fractional coordinates)
         """
+        if remove_oxi:
+            defect_structure.remove_oxidation_states()
+            base_structure.remove_oxidation_states()
+
         if self._is_impurity(defect_structure, base_structure):
             return self.get_impurity_position(defect_structure, base_structure)
         else:

pymatgen/analysis/defects/generators.py

@@ -353,13 +353,15 @@ def __init__(
         min_dist: float = 1.0,
         avg_radius: float = 0.4,
         max_avg_charge: float = 0.9,
+        max_insertions: int | None = None,
     ) -> None:
         self.clustering_tol = clustering_tol
         self.ltol = ltol
         self.stol = stol
         self.angle_tol = angle_tol
         self.avg_radius = avg_radius
         self.max_avg_charge = max_avg_charge
+        self.max_insertions = max_insertions
         super().__init__(min_dist=min_dist)
 
     def generate(self, chgcar: Chgcar, insert_species: set[str] | list[str], **kwargs) -> Generator[Interstitial, None, None]:  # type: ignore[override]
@@ -375,6 +377,8 @@ def generate(self, chgcar: Chgcar, insert_species: set[str] | list[str], **kwarg
         cand_sites_and_mul = [*self._get_candidate_sites(chgcar)]
         for species in insert_species:
             cand_sites = [cand_site for cand_site, mul in cand_sites_and_mul]
+            if self.max_insertions is not None:
+                cand_sites = cand_sites[: self.max_insertions]
             multiplicity = [mul for cand_site, mul in cand_sites_and_mul]
             yield from super().generate(
                 chgcar.structure,
@@ -404,8 +408,20 @@ def generate_all_native_defects(
     sub_generator: SubstitutionGenerator | None = None,
     vac_generator: VacancyGenerator | None = None,
     int_generator: ChargeInterstitialGenerator | None = None,
+    max_insertions: int | None = None,
 ):
-    """Generate all native defects."""
+    """Generate all native defects.
+
+    Convenience function to generate all native defects for a host structure or chgcar object.
+
+    Args:
+        host: The host structure or chgcar object.
+        sub_generator: The substitution generator, if None, a default generator is used.
+        vac_generator: The vacancy generator, if None, a default generator is used.
+        int_generator: The interstitial generator, if None, a default
+            ChargeInterstitialGenerator is used.
+        max_insertions: The maximum number of interstitials to attempt to insert.
+    """
     if isinstance(host, Chgcar):
         struct = host.structure
         chgcar = host
@@ -426,7 +442,9 @@ def generate_all_native_defects(
         yield from sub_generator.generate(struct, {sp2: sp1})
     # generate interstitials if a chgcar is provided
     if chgcar is not None:
-        int_generator = int_generator or ChargeInterstitialGenerator()
+        int_generator = int_generator or ChargeInterstitialGenerator(
+            max_insertions=max_insertions
+        )
         yield from int_generator.generate(chgcar, insert_species=species)
 
 

pymatgen/analysis/defects/thermo.py

@@ -1007,19 +1007,20 @@ def plot_formation_energy_diagrams(
     formation_energy_diagrams = list(filter(filterfunction, formation_energy_diagrams))
 
     band_gap = formation_energy_diagrams[0].band_gap
-    if not xlim and not band_gap:
+    if not xlim and band_gap is None:
         raise ValueError("Must specify xlim or set band_gap attribute")
 
     if axis is None:
         _, axis = plt.subplots()
-    axis.plot([0, 0], [0, 1], color=band_edge_color, linestyle="--", linewidth=1)
-    if not xlim and band_gap:
+    axis.axvline(band_gap, color=band_edge_color, linestyle="--", linewidth=1)
+    axis.axvline(0, color=band_edge_color, linestyle="--", linewidth=1)
+    if not xlim:
         xmin, xmax = np.subtract(-0.2, alignment), np.subtract(
             band_gap + 0.2, alignment
         )
     else:
         xmin, xmax = xlim
-    ymin, ymax = 0.0, 1.0
+    ymin, ymax = np.inf, -np.inf
     legends_txt = []
     artists = []
     fontwidth = 12
@@ -1051,9 +1052,12 @@ def plot_formation_energy_diagrams(
                 lowerlines, np.add(xmin, alignment), np.add(xmax, alignment)
             )
         )
+        trans_y = trans[:, 1]
+        ymin = min(ymin, min(trans_y))
+        ymax = max(ymax, max(trans_y))
         axis.plot(
             np.subtract(trans[:, 0], alignment),
-            trans[:, 1],
+            trans_y,
             color=colors[fid],
             ls=linestyle,
             lw=linewidth,

tests/test_core.py

@@ -1,6 +1,5 @@
 import numpy as np
 from pymatgen.core.periodic_table import Element, Specie
-from pymatgen.io.vasp import Poscar
 
 from pymatgen.analysis.defects.core import (
     Adsorbate,
@@ -10,6 +9,7 @@
     Substitution,
     Vacancy,
 )
+from pymatgen.analysis.defects.finder import DefectSiteFinder
 
 
 def test_vacancy(gan_struct):
@@ -47,10 +47,27 @@ def test_substitution(gan_struct):
     assert sub.element_changes == {Element("N"): -1, Element("O"): 1}
 
     # test supercell with locking
-    sc_locked = sub.get_supercell_structure(relax_radius=1.0)
-    poscar_locked = Poscar(sc_locked)
-    poscar_string = poscar_locked.get_string()
-    assert len(poscar_string.split("\n")) == 138
+    sc_locked = sub.get_supercell_structure(relax_radius=5.0)
+    free_sites = [
+        i
+        for i, site in enumerate(sc_locked)
+        if site.properties["selective_dynamics"][0]
+    ]
+
+    finder = DefectSiteFinder()
+    fpos = finder.get_defect_fpos(sc_locked, sub.structure)
+    cpos = sc_locked.lattice.get_cartesian_coords(fpos)
+    free_sites_ref = sc_locked.get_sites_in_sphere(cpos, 5.0, include_index=True)
+    free_sites_ref = [site.index for site in free_sites_ref]
+    free_sites_union = set(free_sites_ref) | set(free_sites)
+    free_sites_intersection = set(free_sites_ref) & set(free_sites)
+    assert len(free_sites_intersection) / len(free_sites_union) == 1.0
+
+    # test perturbation
+    sc_locked = sub.get_supercell_structure(relax_radius=5.0, perturb=0.0)
+    free_sites_ref2 = sc_locked.get_sites_in_sphere(cpos, 5.0, include_index=True)
+    free_sites_ref2 = [site.index for site in free_sites_ref2]
+    assert set(free_sites_ref2) == set(free_sites_ref)
 
     # test for user defined charge
     dd = sub.as_dict()