Bugfix in LargeSigmaHandler and expand vasp_gam auto checks

src/custodian/vasp/handlers.py

@@ -1365,7 +1365,7 @@ class LargeSigmaHandler(ErrorHandler):
 
     is_monitor: bool = True
 
-    def __init__(self, e_entropy_tol: float = 1e-3, min_sigma: float = 0.01, output_filename: str = "OUTCAR") -> None:
+    def __init__(self, e_entropy_tol: float = 1e-3, min_sigma: float = 0.001, output_filename: str = "OUTCAR") -> None:
         """Initializes the handler with a buffer time."""
         self.e_entropy_tol = e_entropy_tol
         self.min_sigma = min_sigma
@@ -1374,13 +1374,13 @@ def __init__(self, e_entropy_tol: float = 1e-3, min_sigma: float = 0.01, output_
     def check(self, directory="./") -> bool:
         """Check for error."""
         incar = Incar.from_file(os.path.join(directory, "INCAR"))
+        if incar.get("ISMEAR", 1) < 0:
+            # skip check
+            return False
+
         try:
             outcar = load_outcar(os.path.join(directory, self.output_filename))
-        except Exception:
-            # Can't perform check if outcar not valid
-            return False
 
-        if incar.get("ISMEAR", 1) >= 0:
             # get entropy terms, ionic step counts, and number of completed ionic steps
             outcar.read_pattern(
                 {"smearing_entropy": r"entropy T\*S.*= *(\D\d*\.\d*)"},
@@ -1404,7 +1404,10 @@ def check(self, directory="./") -> bool:
             e_step_idx = [step[0] for step in outcar.data.get("electronic_steps", [])]
             smearing_entropy = outcar.data.get("smearing_entropy", [0.0 for _ in e_step_idx])
             for ie_step_idx, ie_step in enumerate(e_step_idx):
-                if ie_step <= completed_ionic_steps:
+                # Because this handler monitors OUTCAR dynamically, it sometimes tries
+                # to retrieve data in OUTCAR before that data is written. To avoid this,
+                # we have two checks for list length here
+                if ie_step <= completed_ionic_steps and ie_step_idx < len(smearing_entropy):
                     entropies_per_atom[ie_step - 1] = smearing_entropy[ie_step_idx]
 
             if len(entropies_per_atom) > 0:
@@ -1413,7 +1416,11 @@ def check(self, directory="./") -> bool:
                 if self.entropy_per_atom > self.e_entropy_tol:
                     return True
 
-        return False
+            return False
+
+        except Exception:
+            # Can't perform check if outcar not valid, or data is missing
+            return False
 
     def correct(self, directory="./"):
         """Perform corrections."""

src/custodian/vasp/jobs.py

@@ -260,8 +260,7 @@ def run(self, directory="./"):
         cmd = list(self.vasp_cmd)
         if self.auto_gamma:
             vi = VaspInput.from_directory(directory)
-            kpts = vi["KPOINTS"]
-            if kpts is not None and kpts.style == Kpoints.supported_modes.Gamma and tuple(kpts.kpts[0]) == (1, 1, 1):
+            if _gamma_point_only_check(vi):
                 if self.gamma_vasp_cmd is not None and which(self.gamma_vasp_cmd[-1]):  # pylint: disable=E1136
                     cmd = self.gamma_vasp_cmd
                 elif which(cmd[-1] + ".gamma"):
@@ -894,12 +893,8 @@ def run(self, directory="./"):
         """
         cmd = list(self.vasp_cmd)
         if self.auto_gamma:
-            kpts = Kpoints.from_file(os.path.join(directory, "KPOINTS"))
-            if kpts.style == Kpoints.supported_modes.Gamma and tuple(kpts.kpts[0]) == (
-                1,
-                1,
-                1,
-            ):
+            vi = VaspInput.from_directory(directory)
+            if _gamma_point_only_check(vi):
                 if self.gamma_vasp_cmd is not None and which(self.gamma_vasp_cmd[-1]):  # pylint: disable=E1136
                     cmd = self.gamma_vasp_cmd
                 elif which(cmd[-1] + ".gamma"):
@@ -987,3 +982,41 @@ def run(self, directory="./") -> None:
 
     def postprocess(self, directory="./") -> None:
         """Dummy postprocess."""
+
+
+def _gamma_point_only_check(vis: VaspInput) -> bool:
+    """
+    Check if only a single k-point is used in this calculation.
+
+    Parameters
+    -----------
+    vis: VaspInput, the VASP input set for the calculation
+
+    Returns:
+    -----------
+    bool: True --> use vasp_gam, False --> use vasp_std
+    """
+    kpts = vis["KPOINTS"]
+    if (
+        kpts is not None
+        and kpts.style == Kpoints.supported_modes.Gamma
+        and tuple(kpts.kpts[0]) == (1, 1, 1)
+        and all(abs(ks) < 1.0e-6 for ks in kpts.kpts_shift)
+    ):
+        return True
+
+    if (kspacing := vis["INCAR"].get("KSPACING")) is not None and vis["INCAR"].get("KGAMMA", True):
+        # Get number of kpoints per axis according to the formula given by VASP:
+        # https://www.vasp.at/wiki/index.php/KSPACING
+        # Note that the VASP definition of the closure relation between reciprocal
+        # lattice vectors b_i and direct lattice vectors a_j is not the conventional
+        # b_i . a_j = 2 pi delta_ij,
+        # and instead places the 2 pi factor in the formula for getting the number
+        # of kpoints per axis.
+        nk = [
+            int(max(1, np.ceil(vis["POSCAR"].structure.lattice.reciprocal_lattice.abc[ik] / kspacing)))
+            for ik in range(3)
+        ]
+        return np.prod(nk) == 1
+
+    return False

tests/vasp/test_handlers.py

@@ -531,7 +531,6 @@ def test_too_large_kspacing(self) -> None:
         handler.check()
         dct = handler.correct()
         assert dct["errors"] == ["dentet"]
-        print(dct["actions"])
         assert dct["actions"] == [{"action": {"_set": {"KSPACING": 1.333333, "KGAMMA": True}}, "dict": "INCAR"}]
 
     def test_nbands_not_sufficient(self) -> None:
@@ -847,6 +846,33 @@ def test_check_correct_large_sigma(self) -> None:
         handler = LargeSigmaHandler(output_filename=zpath("OUTCAR_pass_sigma_check"))
         assert not handler.check()
 
+    def test_no_crash_on_partial_output(self) -> None:
+        from pathlib import Path
+
+        from monty.io import zopen
+        # ensure that the handler doesn't crash when the OUTCAR isn't completely written
+        # this prevents jobs from being killed when the handler itself crashes
+
+        orig_outcar_path = Path(zpath("OUTCAR_pass_sigma_check"))
+        new_outcar_name = str(orig_outcar_path.parent.resolve() / f"temp_{orig_outcar_path.name}")
+        shutil.copy(orig_outcar_path, new_outcar_name)
+
+        # simulate this behavior by manually removing one of the electronic
+        # entropy lines that the handler searches for
+        with zopen(new_outcar_name, "rt") as f:
+            data = f.read().splitlines()
+
+        for rm_idx in range(len(data) - 1, 0, -1):
+            if "T*S" in data[rm_idx]:
+                data.pop(rm_idx)
+                break
+
+        with zopen(new_outcar_name, "wt") as f:
+            f.write("\n".join(data))
+
+        handler = LargeSigmaHandler(output_filename=zpath("OUTCAR_partial_output"))
+        assert not handler.check()
+
 
 class ZpotrfErrorHandlerTest(PymatgenTest):
     def setUp(self) -> None:

tests/vasp/test_jobs.py

@@ -7,9 +7,11 @@
 import pytest
 from monty.os import cd
 from monty.tempfile import ScratchDir
+from pymatgen.core import Structure
 from pymatgen.io.vasp import Incar, Kpoints, Poscar
+from pymatgen.io.vasp.sets import MPRelaxSet
 
-from custodian.vasp.jobs import GenerateVaspInputJob, VaspJob, VaspNEBJob
+from custodian.vasp.jobs import GenerateVaspInputJob, VaspJob, VaspNEBJob, _gamma_point_only_check
 from tests.conftest import TEST_FILES
 
 pymatgen.core.SETTINGS["PMG_VASP_PSP_DIR"] = TEST_FILES
@@ -182,3 +184,33 @@ def test_run(self) -> None:
             assert old_incar["ICHARG"] == 1
             kpoints = Kpoints.from_file("KPOINTS")
             assert str(kpoints.style) == "Reciprocal"
+
+
+class TestAutoGamma:
+    """
+    Test that a VASP job can automatically detect when only 1 k-point at GAMMA is used.
+    """
+
+    def test_gamma_checks(self) -> None:
+        # Isolated atom in PBC
+        structure = Structure(
+            lattice=[[15 + 0.1 * i if i == j else 0.0 for j in range(3)] for i in range(3)],
+            species=["Na"],
+            coords=[[0.5 for _ in range(3)]],
+        )
+
+        vis = MPRelaxSet(structure=structure)
+        assert vis.kpoints.kpts == [(1, 1, 1)]
+        assert vis.kpoints.style.name == "Gamma"
+        assert _gamma_point_only_check(vis.get_input_set())
+
+        # no longer Gamma-centered
+        vis = MPRelaxSet(structure=structure, user_kpoints_settings=Kpoints(kpts_shift=(0.1, 0.0, 0.0)))
+        assert not _gamma_point_only_check(vis.get_input_set())
+
+        # have to increase KSPACING or this will result in a non 1 x 1 x 1 grid
+        vis = MPRelaxSet(structure=structure, user_incar_settings={"KSPACING": 0.5})
+        assert _gamma_point_only_check(vis.get_input_set())
+
+        vis = MPRelaxSet(structure=structure, user_incar_settings={"KSPACING": 0.5, "KGAMMA": False})
+        assert not _gamma_point_only_check(vis.get_input_set())