Repository: materialsproject/custodian
Pull request #362 (2025-02-11 20:12:56)

Cleaner handling of NPAR
## Summary

Anytime that we modify NCORE, we should also unset NPAR if it's present in the INCAR file since NPAR takes precedence. This is done throughout Custodian, but there was one spot missing it. I added it in.


Patch:
diff --git a/src/custodian/vasp/handlers.py b/src/custodian/vasp/handlers.py
index 2b152176..c3ebaaeb 100644
--- a/src/custodian/vasp/handlers.py
+++ b/src/custodian/vasp/handlers.py
@@ -394,6 +394,8 @@ class VaspErrorHandler(ErrorHandler):
                     actions.append({"dict": "INCAR", "action": {"_set": {"POTIM": potim}}})
 
         if self.errors.intersection(["subspacematrix"]):
+            # Sometimes, this error can be due to parallelization issues with running across too many cores
+            # on a small structure. If this is the case, try reducing the number of cores or increasing NCORE or NPAR
             if self.error_count["subspacematrix"] == 0 and vi["INCAR"].get("LREAL", False) is not False:
                 actions.append({"dict": "INCAR", "action": {"_set": {"LREAL": False}}})
             elif self.error_count["subspacematrix"] == 1 and vi["INCAR"].get("PREC", "Normal") != "Accurate":
@@ -570,6 +572,8 @@ class VaspErrorHandler(ErrorHandler):
                 nelect = 1  # dummy value
             if nelect < nprocs:
                 actions.append({"dict": "INCAR", "action": {"_set": {"NCORE": vi["INCAR"].get("NCORE", 1) * 2}}})
+                if "NPAR" in vi["INCAR"]:
+                    actions.append({"dict": "INCAR", "action": {"_unset": {"NPAR": 1}}})
 
         if "grad_not_orth" in self.errors:
             # Often coincides with algo_tet, in which the algo_tet error handler will also resolve grad_not_orth.
@@ -1365,7 +1369,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 +1378,13 @@ class LargeSigmaHandler(ErrorHandler):
     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 +1408,10 @@ class LargeSigmaHandler(ErrorHandler):
             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 +1420,11 @@ class LargeSigmaHandler(ErrorHandler):
                 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."""
diff --git a/src/custodian/vasp/jobs.py b/src/custodian/vasp/jobs.py
index a94957c2..4907ed16 100644
--- a/src/custodian/vasp/jobs.py
+++ b/src/custodian/vasp/jobs.py
@@ -260,8 +260,7 @@ class VaspJob(Job):
         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 @@ class VaspNEBJob(VaspJob):
         """
         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 @@ class GenerateVaspInputJob(Job):
 
     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


Test patch:
diff --git a/tests/vasp/test_handlers.py b/tests/vasp/test_handlers.py
index 6147862f..2904ddb8 100644
--- a/tests/vasp/test_handlers.py
+++ b/tests/vasp/test_handlers.py
@@ -531,7 +531,6 @@ class VaspErrorHandlerTest(PymatgenTest):
         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 @@ class LargeSigmaHandlerTest(PymatgenTest):
         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:
diff --git a/tests/vasp/test_jobs.py b/tests/vasp/test_jobs.py
index bcbebce7..7e80caa8 100644
--- a/tests/vasp/test_jobs.py
+++ b/tests/vasp/test_jobs.py
@@ -7,9 +7,11 @@ import pymatgen
 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 @@ class TestGenerateVaspInputJob:
             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())