fix(vasp): skip malformed OUTCAR force frames (#1030)

Author: njzjz-bot

dpdata/formats/vasp/outcar.py

@@ -126,6 +126,10 @@ def get_outcar_block(fp, ml=False):
     return blk
 
 
+class _IncompleteForceTableError(ValueError):
+    """Indicate that a ``TOTAL-FORCE`` table ended before all atom rows."""
+
+
 def check_outputs(coord, cell, force):
     if len(force) == 0:
         raise ValueError("cannot find forces in OUTCAR block")
@@ -167,9 +171,21 @@ def _get_frames_lower(fp, fname, begin=0, step=1, ml=False, convergence_check=Tr
     rec_failed = []
     while len(blk) > 0:
         if cc >= begin and (cc - begin) % step == 0:
-            coord, cell, energy, force, virial, is_converge = analyze_block(
-                blk, ntot, nelm, ml
-            )
+            try:
+                coord, cell, energy, force, virial, is_converge = analyze_block(
+                    blk, ntot, nelm, ml
+                )
+            except _IncompleteForceTableError as exc:
+                # An interrupted VASP process can leave one damaged ionic
+                # step between otherwise usable frames.  Its atom arrays are
+                # incomplete, so skip precisely that step instead of aborting
+                # the complete trajectory.
+                warnings.warn(
+                    f"incomplete labels in frame {cc + 1}; it is ignored: {exc}"
+                )
+                blk = get_outcar_block(fp, ml)
+                cc += 1
+                continue
             if energy is None:
                 break
             if nwrite == 0:
@@ -267,9 +283,33 @@ def analyze_block(lines, ntot, nelm, ml=False):
             virial[0][2] = tmp_v[5]
             virial[2][0] = tmp_v[5]
         elif "TOTAL-FORCE" in ii and (("ML" in ii) == ml):
+            block_coord = []
+            block_force = []
             for jj in range(idx + 2, idx + 2 + ntot):
-                tmp_l = lines[jj]
-                info = [float(ss) for ss in tmp_l.split()]
-                coord.append(info[:3])
-                force.append(info[3:6])
+                if jj >= len(lines):
+                    # VASP output can end while a force table is being
+                    # written.  Do not index beyond the block or retain the
+                    # preceding subset of atoms as a complete frame.
+                    raise _IncompleteForceTableError(
+                        f"expected {ntot} atom rows, found {len(block_force)}"
+                    )
+                fields = lines[jj].split()
+                if len(fields) < 6:
+                    raise _IncompleteForceTableError(
+                        f"expected 6 numeric columns, got {len(fields)} in {lines[jj]!r}"
+                    )
+                try:
+                    info = [float(ss) for ss in fields[:6]]
+                except ValueError as exc:
+                    # A non-numeric record with at least six fields, such as
+                    # the following energy record, means the table ended
+                    # before ``ntot`` atom rows.  Parse into temporary lists
+                    # so none of the incomplete table leaks into the frame.
+                    raise _IncompleteForceTableError(
+                        f"non-numeric atom row {lines[jj]!r}"
+                    ) from exc
+                block_coord.append(info[:3])
+                block_force.append(info[3:6])
+            coord = block_coord
+            force = block_force
     return coord, cell, energy, force, virial, is_converge

tests/test_vasp_outcar.py

@@ -1,11 +1,14 @@
 from __future__ import annotations
 
+import io
 import unittest
+import warnings
 
 import numpy as np
 from comp_sys import CompLabeledSys, IsPBC
 from context import dpdata
 
+from dpdata.formats.vasp.outcar import _get_frames_lower
 from dpdata.utils import uniq_atom_names
 
 
@@ -21,6 +24,83 @@ def setUp(self):
         self.v_places = 4
 
 
+class TestVaspOUTCARIncompleteForceTable(unittest.TestCase):
+    """Regression tests for force tables cut short before all atom rows."""
+
+    @staticmethod
+    def _block(force_rows, *, include_header=False, energy=-1.0):
+        lines = []
+        if include_header:
+            lines.extend(
+                [
+                    " TITEL  = PAW_PBE H 15Jun2001",
+                    " NELM = 60; maximum number of electronic SC steps",
+                    " ions per type = 2",
+                ]
+            )
+        lines.extend(
+            [
+                " VOLUME and BASIS-vectors are now :",
+                " filler",
+                " filler",
+                " filler",
+                " filler",
+                " 1.0 0.0 0.0",
+                " 0.0 1.0 0.0",
+                " 0.0 0.0 1.0",
+                " POSITION                                       TOTAL-FORCE (eV/Angst)",
+                " -----------------------------------------------------------------------------------",
+                *force_rows,
+                f" free  energy   TOTEN  =       {energy:.6f} eV",
+            ]
+        )
+        return "\n".join(lines) + "\n"
+
+    def test_non_numeric_record_before_all_atoms_skips_incomplete_frame(self):
+        # The first frame is valid, but the second table reaches its energy
+        # record after only one of the two expected atoms. Older tests used
+        # complete tables, so converting that non-numeric record was never
+        # exercised.
+        valid_rows = ["0 0 0 1 2 3", "1 1 1 4 5 6"]
+        incomplete_rows = ["2 2 2 7 8 9"]
+        contents = self._block(valid_rows, include_header=True) + self._block(
+            incomplete_rows, energy=-2.0
+        )
+
+        with warnings.catch_warnings(record=True) as caught:
+            warnings.simplefilter("always")
+            frames = _get_frames_lower(io.StringIO(contents), "OUTCAR")
+
+        self.assertEqual(frames[3].shape, (1, 3, 3))
+        self.assertEqual(frames[4].shape, (1, 2, 3))
+        self.assertEqual(frames[6].shape, (1, 2, 3))
+        self.assertTrue(
+            any("incomplete labels in frame 2" in str(item.message) for item in caught)
+        )
+
+    def test_truncated_table_does_not_index_past_block(self):
+        # A file truncated immediately after its first atom previously raised
+        # IndexError while the parser blindly indexed all ``ntot`` rows.
+        valid_rows = ["0 0 0 1 2 3", "1 1 1 4 5 6"]
+        truncated_block = "\n".join(
+            [
+                " POSITION                                       TOTAL-FORCE (eV/Angst)",
+                " -----------------------------------------------------------------------------------",
+                "0 0 0 1 2 3",
+            ]
+        )
+        contents = self._block(valid_rows, include_header=True) + truncated_block
+
+        with warnings.catch_warnings(record=True) as caught:
+            warnings.simplefilter("always")
+            frames = _get_frames_lower(io.StringIO(contents), "OUTCAR")
+
+        self.assertEqual(frames[4].shape, (1, 2, 3))
+        self.assertTrue(
+            any("expected 2 atom rows, found 1" in str(item.message) for item in caught)
+        )
+
+
 class TestVaspOUTCARTypeMap(unittest.TestCase, CompLabeledSys, IsPBC):
     def setUp(self):
         sys0 = dpdata.LabeledSystem("poscars/OUTCAR.ch4.unconverged", fmt="vasp/outcar")