Sc back (#171)

* working

* working

* working

* working

* working

* pyrho

* pyrho

* data

* pyrho

* up_sample

* up_sample

* up_sample

* up_sample

* up_sample

* update

* update

Author: jmmshn

.github/workflows/testing.yml

@@ -66,26 +66,6 @@ jobs:
           token: ${{ secrets.CODECOV_TOKEN }}
           file: ./coverage.xml
 
-  # figure out if it is possible to upload to testpypi
-  # deploy-test:
-  #   runs-on: ubuntu-latest
-  #   steps:
-  #   - uses: actions/checkout@v3
-  #   - name: Set up Python
-  #     uses: actions/setup-python@v4
-  #     with:
-  #       python-version: '3.x'
-  #   - name: Install dependencies
-  #     run: |
-  #       python -m pip install --upgrade pip
-  #       pip install build
-  #   - name: Build package
-  #     run: python -m build
-  #   - name: Publish package to TestPyPI
-  #     uses: pypa/gh-action-pypi-publish@release/v1
-  #     with:
-  #       password: ${{ secrets.TEST_PYPI_API_TOKEN }}
-  #       repository-url: https://test.pypi.org/legacy/
 
   docs:
     runs-on: ubuntu-latest
@@ -107,6 +87,7 @@ jobs:
 
       - name: Install dependencies
         run: |
+          python -m pip install --upgrade pip
           pip install -e .[strict]
           pip install -e .[docs]
       - name: Build

.pre-commit-config.yaml

@@ -3,15 +3,15 @@ exclude: ^(tests)
 repos:
 - repo: https://github.com/astral-sh/ruff-pre-commit
   # Ruff version.
-  rev: v0.1.8
+  rev: v0.2.2
   hooks:
     # Run the linter.
     - id: ruff
       args: [ --fix ]
     # Run the formatter.
     - id: ruff-format
 - repo: https://github.com/pre-commit/mirrors-mypy
-  rev: v1.5.1
+  rev: v1.8.0
   hooks:
   - id: mypy
     files: ^pymatgen/

docs/source/content/defining-defects.ipynb

@@ -222,8 +222,8 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "#sc_struct_smaller = mg_ga_defect0.get_supercell_structure(max_atoms=100)\n",
-    "#sc_struct_smaller.num_sites"
+    "# sc_struct_smaller = mg_ga_defect0.get_supercell_structure(max_atoms=100)\n",
+    "# sc_struct_smaller.num_sites"
    ]
   },
   {
@@ -234,9 +234,9 @@
     "\n",
     "Intersitial defects are usually hard to define due to a lack of reference points for the site.\n",
     "Extensive symmetry can be done to identifiy highly symmetric sites in the structure for interstitial insertion.\n",
-    "However, the recommended method to create interstitial defects is to use the `ChargeInterstitialGenerator` which analyzes the charge density to identify interstitial sites.  The code snippet to generate the interstitial sites are given below.\n",
+    "However, the recommended method to create interstitial defects is to use the `ChargeInterstitialGenerator` which analyzes the charge density to identify interstitial sites.  The code snippet to generate the interstitial sites is given below.\n",
     "\n",
-    "For more details, check out this [paper](https://www.nature.com/articles/s41524-020-00422-3)\n"
+    "For more details on how the intersitial site identification works, check out this [paper](https://www.nature.com/articles/s41524-020-00422-3)."
    ]
   },
   {
@@ -247,8 +247,11 @@
    },
    "outputs": [],
    "source": [
+    "from pymatgen.analysis.defects.generators import (\n",
+    "    ChargeInterstitialGenerator,\n",
+    "    generate_all_native_defects,\n",
+    ")\n",
     "from pymatgen.io.vasp import Chgcar\n",
-    "from pymatgen.analysis.defects.generators import ChargeInterstitialGenerator, generate_all_native_defects\n",
     "\n",
     "chgcar = Chgcar.from_file(TEST_FILES / \"CHGCAR.Fe3O4.vasp\")\n",
     "cig = ChargeInterstitialGenerator()\n",

pymatgen/analysis/defects/ccd.py

@@ -4,7 +4,7 @@
 import logging
 from dataclasses import dataclass
 from itertools import groupby
-from typing import TYPE_CHECKING, Optional, Sequence, Tuple
+from typing import TYPE_CHECKING
 
 import numpy as np
 from monty.json import MSONable
@@ -20,6 +20,7 @@
 if TYPE_CHECKING:
     from ctypes import Structure
     from pathlib import Path
+    from typing import Optional, Sequence, Tuple
 
     import numpy.typing as npt
     from matplotlib.axes import Axes
@@ -65,9 +66,9 @@ class HarmonicDefect(MSONable):
     charge_state: int
     ispin: int
     vrun: Optional[Vasprun] = None
-    distortions: Optional[list[float]] = None
-    structures: Optional[list[Structure]] = None
-    energies: Optional[list[float]] = None
+    distortions: Optional[Sequence[float]] = None
+    structures: Optional[Sequence[Structure]] = None
+    energies: Optional[Sequence[float]] = None
     defect_band: Optional[Sequence[tuple]] = None
     relaxed_index: Optional[int] = None
     relaxed_bandstructure: Optional[BandStructure] = None
@@ -234,7 +235,7 @@ def _parse_vasprun(vasprun: Vasprun):
             omega=omega,
             charge_state=charge_state,
             ispin=ispin,
-            structures=structures,
+            structures=list(structures),
             distortions=distortions,
             energies=energies,
             defect_band=defect_band,
@@ -313,7 +314,7 @@ def occupation(self, t: npt.ArrayLike | float) -> npt.ArrayLike:
         return 1.0 / (1 - np.exp(-self.omega_eV / KB * t))
 
     def read_wswqs(
-        self, directory: Path, distortions: list[float] | None = None
+        self, directory: Path, distortions: Sequence[float] | None = None
     ) -> None:
         """Read the WSWQ files from a directory.
 

pymatgen/analysis/defects/corrections/freysoldt.py

@@ -412,8 +412,8 @@ def plot_plnr_avg(plot_data, title=None, saved=False, ax=None):
     )
 
     ax.set_xlim(round(x[0]), round(x[-1]))
-    ymin = min(min(v_R), min(dft_diff), min(short_range))
-    ymax = max(max(v_R), max(dft_diff), max(short_range))
+    ymin = min(v_R + dft_diff + short_range)
+    ymax = max(v_R + dft_diff + short_range)
     ax.set_ylim(-0.2 + ymin, 0.2 + ymax)
     ax.set_xlabel(r"distance along axis ($\AA$)", fontsize=15)
     ax.set_ylabel("Potential (V)", fontsize=15)

pymatgen/analysis/defects/generators.py

@@ -21,6 +21,8 @@
 from pymatgen.symmetry.analyzer import SpacegroupAnalyzer
 
 if TYPE_CHECKING:
+    from typing import Sequence
+
     from pymatgen.analysis.defects.core import Defect
     from pymatgen.io.vasp import VolumetricData
 
@@ -86,7 +88,10 @@ def __init__(
         self.angle_tolerance = angle_tolerance
 
     def generate(
-        self, structure: Structure, rm_species: list[str | Species] = None, **kwargs
+        self,
+        structure: Structure,
+        rm_species: list[str | Species] | None = None,
+        **kwargs,
     ) -> Generator[Vacancy, None, None]:
         """Generate a vacancy defects.
 
@@ -254,9 +259,10 @@ def __init__(self, min_dist: float = 0.5) -> None:
     def generate(
         self,
         structure: Structure,
-        insertions: dict[str, list[list[float]]],
-        multiplicities: dict[str, list[int]] | None = None,
-        equivalent_positions: dict[str, list[list[list[float]]]] | None = None,
+        insertions: dict[str, Sequence[Sequence[float]]],
+        multiplicities: dict[str, Sequence[int]] | None = None,
+        equivalent_positions: dict[str, Sequence[Sequence[Sequence[float]]]]
+        | None = None,
         **kwargs,
     ) -> Generator[Interstitial, None, None]:
         """Generate interstitials.
@@ -306,8 +312,8 @@ def generate(
                 )
 
     def _filter_colliding(
-        self, fcoords: list[list[float]], structure: Structure
-    ) -> Generator[tuple[int, list[float]], None, None]:
+        self, fcoords: Sequence[Sequence[float]], structure: Structure
+    ) -> Generator[tuple[int, Sequence[float]], None, None]:
         """Check the sites for collisions.
 
         Args:

pymatgen/analysis/defects/plotting/optics.py

@@ -27,9 +27,9 @@ def plot_optical_transitions(
     kpt_index: int = 0,
     band_window: int = 5,
     user_defect_band: tuple = tuple(),
-    other_defect_bands: list[int] = None,
-    ijdirs: list[tuple] = None,
-    shift_eig: dict[tuple, float] = None,
+    other_defect_bands: list[int] | None = None,
+    ijdirs: list[tuple] | None = None,
+    shift_eig: dict[tuple, float] | None = None,
     x0: float = 0,
     x_width: float = 2,
     ax=None,
@@ -118,9 +118,9 @@ def get_bs_eigenvalues(
     defect: HarmonicDefect,
     kpt_index: int = 0,
     band_window: int = 5,
-    user_defect_band: tuple = None,
-    other_defect_bands: list[int] = None,
-    shift_eig: dict[tuple, float] = None,
+    user_defect_band: tuple | None = None,
+    other_defect_bands: list[int] | None = None,
+    shift_eig: dict[tuple, float] | None = None,
 ) -> dict[tuple, float]:
     """Read the eigenvalues from `HarmonicDefect.relaxed_bandstructure`.
 

pymatgen/analysis/defects/plotting/phases.py

@@ -90,8 +90,8 @@ def plot_chempot_2d(
     )
     ax.add_patch(patch)
 
-    ax.set_xlabel(f"$\Delta\mu_{{{x_element}}}$ (eV)")
-    ax.set_ylabel(f"$\Delta\mu_{{{y_element}}}$ (eV)")
+    ax.set_xlabel(rf"$\Delta\mu_{{{x_element}}}$ (eV)")
+    ax.set_ylabel(rf"$\Delta\mu_{{{y_element}}}$ (eV)")
     ax.set_xlim(x_min - PLOT_PADDING, 0 + PLOT_PADDING)
     ax.set_ylim(y_min - PLOT_PADDING, 0 + PLOT_PADDING)
     if label_lines:
@@ -102,7 +102,7 @@ def _convex_hull_2d(
     points: list[dict],
     x_element: Element,
     y_element: Element,
-    competing_phases: list = None,
+    competing_phases: list | None = None,
 ) -> list[dict]:
     """Compute the convex hull of a set of points in 2D.
 

pymatgen/analysis/defects/supercells.py

@@ -6,14 +6,21 @@
 import math
 from typing import TYPE_CHECKING
 
+import numpy as np
+from monty.dev import deprecated
 from pymatgen.analysis.structure_matcher import ElementComparator, StructureMatcher
+from pymatgen.core import Lattice
+from pymatgen.util.coord_cython import is_coord_subset_pbc, pbc_shortest_vectors
+from pyrho.charge_density import ChargeDensity
 
 # from ase.build import find_optimal_cell_shape, get_deviation_from_optimal_cell_shape
 # from pymatgen.io.ase import AseAtomsAdaptor
 
 if TYPE_CHECKING:
+    import numpy.typing as npt
     from numpy.typing import ArrayLike, NDArray
     from pymatgen.core import Structure
+    from pymatgen.io.vasp.outputs import VolumetricData
 
 __author__ = "Jimmy-Xuan Shen"
 __copyright__ = "Copyright 2022, The Materials Project"
@@ -58,12 +65,13 @@ def get_sc_fromstruct(
     return sc_mat
 
 
-def get_matched_structure_mapping(
+def get_matched_structure_mapping_old(
     uc_struct: Structure, sc_struct: Structure, sm: StructureMatcher | None = None
 ):
     """Get the mapping of the supercell to the unit cell.
 
-    Get the mapping from the supercell defect structure onto the base structure,
+    Get the mapping from the supercell structure onto the base structure,
+    Note: this only works for structures that are exactly matched.
 
     Args:
         uc_struct: host structure, smaller cell
@@ -86,6 +94,38 @@ def get_matched_structure_mapping(
     return sc_m, total_t
 
 
+@deprecated(message="This function was reworked in Feb 2024")
+def get_matched_structure_mapping(
+    uc_struct: Structure, sc_struct: Structure, sm: StructureMatcher | None = None
+):
+    """Get the mapping of the supercell to the unit cell.
+
+    Get the mapping from the supercell structure onto the base structure,
+    Note: this only works for structures that are exactly matched.
+
+    Args:
+        uc_struct: host structure, smaller cell
+        sc_struct: bigger cell
+        sm: StructureMatcher instance
+    Returns:
+        sc_m : supercell matrix to apply to s1 to get s2
+        total_t : translation to apply on s1 * sc_m to get s2
+    """
+    if sm is None:
+        sm = StructureMatcher(
+            primitive_cell=False, comparator=ElementComparator(), attempt_supercell=True
+        )
+    s1, s2 = sm._process_species([sc_struct.copy(), uc_struct.copy()])
+    trans = sm.get_transformation(s1, s2)
+    if trans is None:
+        return None
+    sc, t, mapping = trans
+    temp = s2.copy().make_supercell(sc)
+    ii, jj = 0, mapping[0]
+    vec = np.round(sc_struct[ii].frac_coords - temp[jj].frac_coords)
+    return sc, t + vec
+
+
 def _cubic_cell(
     base_struct: Structure,
     min_atoms: int = 80,
@@ -159,3 +199,66 @@ def _ase_cubic(base_structure, min_atoms: int = 80, max_atoms: int = 240):
     if min_dev[1] is None:
         raise RuntimeError("Could not find a cubic supercell")
     return min_dev[1]
+
+
+def _avg_lat(l1, l2):
+    """Get the average lattice from two lattices."""
+    params = (np.array(l1.parameters) + np.array(l2.parameters)) / 2
+    return Lattice.from_parameters(*params)
+
+
+def _lowest_dist(struct, ref_struct):
+    """For each site, return the lowest distance to any site in the reference structure."""
+    avg_lat = _avg_lat(struct.lattice, ref_struct.lattice)
+    _, d_2 = pbc_shortest_vectors(
+        avg_lat, struct.frac_coords, ref_struct.frac_coords, return_d2=True
+    )
+    return np.min(d_2, axis=1)
+
+
+def get_closest_sc_mat(
+    uc_struct: Structure,
+    sc_struct: Structure,
+    sm: StructureMatcher | None = None,
+    debug: bool = False,
+):
+    """Get the best guess for the supercell matrix that created this defect cell.
+
+    Args:
+        uc_struct: unit cell structure, should be the host structure
+        sc_struct: supercell structure, should be the defect structure
+        sm: StructureMatcher instance, if None, one will be created with default settings
+        debug: bool, if True, return the full list of (distances, lattice, sc_mat) will
+            be returned
+
+    Returns:
+        sc_mat: supercell matrix to apply to s1 to get s2
+        dist: mean distance between the two structures
+    """
+    if sm is None:
+        sm = StructureMatcher(primitive_cell=False, comparator=ElementComparator())
+
+    fu = int(np.round(sc_struct.lattice.volume / uc_struct.lattice.volume))
+    candidate_lattices = tuple(
+        sm._get_lattices(sc_struct.lattice, uc_struct, supercell_size=fu)
+    )
+
+    def _get_mean_dist(lattice, sc_mat):
+        if (
+            np.dot(np.cross(lattice.matrix[0], lattice.matrix[1]), lattice.matrix[2])
+            < 0
+        ):
+            return float("inf")
+        sc2 = uc_struct * sc_mat
+        return np.mean(_lowest_dist(sc2, sc_struct))
+
+    _, best_sc_mat = min(candidate_lattices, key=lambda x: _get_mean_dist(x[0], x[1]))
+    if debug:
+        return sorted(
+            [
+                (_get_mean_dist(lat_, smat_), lat_, smat_)
+                for lat_, smat_ in candidate_lattices
+            ],
+            key=lambda x: x[0],
+        )
+    return best_sc_mat