Patch locpot (#182)

* sc locpot

* ccd bug

* volumetric data

Author: jmmshn

pymatgen/analysis/defects/core.py

@@ -164,7 +164,12 @@ def get_charge_states(self, padding: int = 1) -> list[int]:
         if isinstance(self.oxi_state, int) or self.oxi_state.is_integer():
             oxi_state = int(self.oxi_state)
         else:  # pragma: no cover
-            raise ValueError("Oxidation state must be an integer")
+            sign = -1 if self.oxi_state < 0 else 1
+            oxi_state = sign * int(np.ceil(abs(self.oxi_state)))
+            _logger.warn(
+                "Non-integer oxidation state detected."
+                f"Rounding to integer with larger absolute value: {self.oxi_state} -> {oxi_state}"
+            )
 
         if oxi_state >= 0:
             charges = [*range(-padding, oxi_state + padding + 1)]

pymatgen/analysis/defects/corrections/freysoldt.py

@@ -18,7 +18,7 @@
     generate_reciprocal_vectors_squared,
     hart_to_ev,
 )
-from pymatgen.io.vasp.outputs import Locpot
+from pymatgen.io.vasp.outputs import Locpot, VolumetricData
 from scipy import stats
 
 if TYPE_CHECKING:
@@ -104,7 +104,7 @@ def get_freysoldt_correction(
 
     q_model = QModel() if q_model is None else q_model
 
-    if isinstance(defect_locpot, Locpot):
+    if isinstance(defect_locpot, VolumetricData):
         list_axis_grid = [*map(defect_locpot.get_axis_grid, [0, 1, 2])]
         list_defect_plnr_avg_esp = [
             *map(defect_locpot.get_average_along_axis, [0, 1, 2])
@@ -130,7 +130,7 @@ def get_freysoldt_correction(
         raise ValueError("defect_locpot must be of type Locpot or dict")
 
     # TODO this can be done with regridding later
-    if isinstance(bulk_locpot, Locpot):
+    if isinstance(bulk_locpot, VolumetricData):
         list_bulk_plnr_avg_esp = [*map(bulk_locpot.get_average_along_axis, [0, 1, 2])]
     elif isinstance(bulk_locpot, dict):
         bulk_locpot_ = {int(k): v for k, v in bulk_locpot.items()}

pymatgen/analysis/defects/supercells.py

@@ -161,11 +161,15 @@ def _cubic_cell(
     try:
         cst.apply_transformation(base_struct)
     except BaseException:
-        return _ase_cubic(base_struct, min_atoms, max_atoms)
+        return _ase_cubic(
+            base_struct, min_atoms=min_atoms, max_atoms=max_atoms, min_length=min_length
+        )
     return cst.transformation_matrix
 
 
-def _ase_cubic(base_structure, min_atoms: int = 80, max_atoms: int = 240):
+def _ase_cubic(
+    base_structure, min_atoms: int = 80, max_atoms: int = 240, min_length=10.0
+):
     """Generate the best supercell from a unit cell.
 
     Use ASE's find_optimal_cell_shape function to find the best supercell.
@@ -174,6 +178,7 @@ def _ase_cubic(base_structure, min_atoms: int = 80, max_atoms: int = 240):
         base_structure: structure of the unit cell
         max_atoms: Maximum number of atoms allowed in the supercell.
         min_atoms: Minimum number of atoms allowed in the supercell.
+        min_length: Minimum length of the smallest supercell lattice vector.
 
     Returns:
         3x3 matrix: supercell matrix
@@ -194,6 +199,10 @@ def _ase_cubic(base_structure, min_atoms: int = 80, max_atoms: int = 240):
             ase_atoms.cell, target_size=size, target_shape="sc"
         )
         sc_cell = aaa.get_atoms(base_structure * sc).cell
+        lattice_lens = np.linalg.norm(sc_cell, axis=1)
+        _logger.warn(f"{lattice_lens}, {min_length}, {min_dev}")
+        if min(lattice_lens) < min_length:
+            continue
         deviation = get_deviation_from_optimal_cell_shape(sc_cell, target_shape="sc")
         min_dev = min(min_dev, (deviation, sc))
     if min_dev[1] is None:

pymatgen/analysis/defects/thermo.py

@@ -23,7 +23,7 @@
 from pymatgen.core import Composition, Element
 from pymatgen.electronic_structure.dos import FermiDos
 from pymatgen.entries.computed_entries import ComputedEntry
-from pymatgen.io.vasp import Locpot, Vasprun
+from pymatgen.io.vasp import Chgcar, Locpot, Vasprun, VolumetricData
 from pyrho.charge_density import get_volumetric_like_sc
 from scipy.constants import value as _cd
 from scipy.optimize import bisect
@@ -148,12 +148,12 @@ def get_freysoldt_correction(
         else:  # pragma: no cover
             defect_fpos = self.sc_defect_frac_coords
 
-        if isinstance(defect_locpot, Locpot):
+        if isinstance(defect_locpot, VolumetricData):
             defect_gn = defect_locpot.dim
         elif isinstance(defect_locpot, dict):
             defect_gn = tuple(map(len, (defect_locpot for k in ["0", "1", "2"])))
 
-        if isinstance(bulk_locpot, Locpot):
+        if isinstance(bulk_locpot, VolumetricData):
             bulk_sc_locpot = get_sc_locpot(
                 uc_locpot=bulk_locpot,
                 defect_struct=defect_struct,