fix(vasp): warn when POTCAR and POSCAR species order disagree

dpdata/formats/vasp/outcar.py

@@ -28,6 +28,99 @@ def atom_name_from_potcar_string(instr: str) -> str:
         return instr
 
 
+# VASP echoes the POSCAR title into a fixed-width field, so a long formula is
+# cut off; anything reaching the full width may have lost its last species.
+_POSCAR_TITLE_WIDTH = 40
+_FORMULA_TOKEN = re.compile(r"^([A-Z][a-z]?)(\d*)$")
+
+
+def composition_from_poscar_title(title: str) -> dict[str, int] | None:
+    """Read a composition from a POSCAR title when every count is explicit.
+
+    The title is free-form text, so this returns ``None`` unless every token
+    reads as an element symbol with a count -- ``Li3 F39 K3`` yields
+    ``{"Li": 3, "F": 39, "K": 3}``, while both ``H C`` and ``POSCAR file
+    written by OVITO`` yield ``None``. Requiring counts avoids interpreting a
+    comment that merely lists elements as a declaration of their order. A
+    title filling the whole field is treated as truncated and its final token
+    is dropped.
+
+    Parameters
+    ----------
+    title : str
+        the text following ``POSCAR =`` in the OUTCAR
+
+    Returns
+    -------
+    Optional[dict[str, int]]
+        element counts, or None if the title is not an explicit composition
+    """
+    from dpdata.periodic_table import ELEMENTS
+
+    truncated = len(title.rstrip()) >= _POSCAR_TITLE_WIDTH
+    tokens = title.split()
+    if truncated:
+        tokens = tokens[:-1]
+    if len(tokens) < 2:
+        # A single symbol cannot disagree on ordering, and a one-word title is
+        # far more likely to be prose than a formula.
+        return None
+    composition = {}
+    for token in tokens:
+        matched = _FORMULA_TOKEN.match(token)
+        if matched is None or matched.group(1) not in ELEMENTS or not matched.group(2):
+            return None
+        name = matched.group(1)
+        composition[name] = composition.get(name, 0) + int(matched.group(2))
+    return composition
+
+
+def check_potcar_poscar_order(
+    atom_names: list[str], atom_numbs: list[int], poscar_title: str | None
+) -> None:
+    """Warn when POTCAR-paired counts contradict a POSCAR-title composition.
+
+    VASP pairs the ``ions per type`` counts, which come from the POSCAR, with
+    the species order of the POTCAR. When the two files disagree, VASP neither
+    reorders nor complains, so counts can silently land on the wrong elements
+    and dpdata faithfully reports the mislabeled system. The POSCAR title is
+    only a comment, however, so compare compositions rather than token order.
+    """
+    if poscar_title is None:
+        return
+    title_composition = composition_from_poscar_title(poscar_title)
+    if title_composition is None:
+        return
+
+    potcar_composition = {}
+    for name, count in zip(atom_names, atom_numbs):
+        potcar_composition[name] = potcar_composition.get(name, 0) + count
+
+    truncated = len(poscar_title.rstrip()) >= _POSCAR_TITLE_WIDTH
+    if truncated:
+        matches = all(
+            potcar_composition.get(name) == count
+            for name, count in title_composition.items()
+        )
+    else:
+        matches = title_composition == potcar_composition
+    if matches:
+        return
+
+    def format_composition(composition: dict[str, int]) -> str:
+        return " ".join(f"{name}{count}" for name, count in composition.items())
+
+    warnings.warn(
+        "the composition produced by pairing the POTCAR species with 'ions per "
+        f"type' in this OUTCAR ({format_composition(potcar_composition)}) does "
+        "not match the explicit composition in the POSCAR title "
+        f"({format_composition(title_composition)}). The atom names reported "
+        "here are what VASP actually computed; if the title describes the "
+        "intended structure, the POTCAR may have been concatenated in a "
+        "different order and the calculation used the wrong potentials."
+    )
+
+
 def system_info(
     lines: list[str],
     type_idx_zero: bool = False,
@@ -59,6 +152,7 @@ def system_info(
     atom_numbs = None
     nelm = None
     nwrite = None
+    poscar_title = None
     for ii in lines:
         if "TITEL" in ii:
             # get atom names from POTCAR info, tested only for PAW_PBE ...
@@ -79,6 +173,12 @@ def system_info(
             m = re.search(r"NWRITE\s*=\s*(\d+)", ii)
             if m:
                 nwrite = int(m.group(1))
+        if poscar_title is None:
+            # the POSCAR title echoed among the start parameters, e.g.
+            # POSCAR =  Li3 F39 K3 Mg3 Ca3 Na3 Al10 O6
+            m = re.match(r"\s*POSCAR\s*=\s*(.*)$", ii)
+            if m:
+                poscar_title = m.group(1)
         if "ions per type" in ii:
             atom_numbs_ = [int(s) for s in ii.split()[4:]]
             if atom_numbs is None:
@@ -103,6 +203,7 @@ def system_info(
             f"Please try to convert data from vasprun.xml instead."
         )
     atom_names = atom_names[: len(atom_numbs)]
+    check_potcar_poscar_order(atom_names, atom_numbs, poscar_title)
     atom_types = []
     for idx, ii in enumerate(atom_numbs):
         for jj in range(ii):

tests/test_vasp_outcar_species_order.py

@@ -0,0 +1,178 @@
+from __future__ import annotations
+
+import os
+import shutil
+import tempfile
+import unittest
+import warnings
+
+from context import dpdata
+
+from dpdata.formats.vasp.outcar import (
+    check_potcar_poscar_order,
+    composition_from_poscar_title,
+)
+
+SOURCE = os.path.join("poscars", "OUTCAR.ch4.1step")
+
+
+class TestPoscarTitleFormula(unittest.TestCase):
+    """The POSCAR title is free-form text and only sometimes a formula."""
+
+    def test_formula_with_counts(self):
+        self.assertEqual(
+            composition_from_poscar_title(" Li3 F39 K3 Mg3 Ca3 Na3 Al10 O6"),
+            {"Li": 3, "F": 39, "K": 3, "Mg": 3, "Ca": 3, "Na": 3, "Al": 10, "O": 6},
+        )
+
+    def test_elements_without_counts_are_ambiguous(self):
+        self.assertIsNone(composition_from_poscar_title("  H  C  "))
+
+    def test_prose_is_not_a_formula(self):
+        for title in (
+            " POSCAR file written by OVITO",
+            " File generated by python recompute scrip",
+            " generated by SeeK-path",
+        ):
+            with self.subTest(title=title):
+                self.assertIsNone(composition_from_poscar_title(title))
+
+    def test_unknown_symbol_is_rejected(self):
+        self.assertIsNone(composition_from_poscar_title(" Xx2 Yy3"))
+
+    def test_single_token_is_rejected(self):
+        # One species cannot disagree on ordering, and a one-word title is far
+        # more likely to be prose.
+        self.assertIsNone(composition_from_poscar_title(" Si"))
+
+    def test_truncated_title_drops_the_cut_token(self):
+        # VASP pads the title to 40 characters, cutting the last species.
+        title = " Cr3 Mn5 Fe3 Co1 Ni4 Cu4 Zn2 Y5 Zr4 Nb5 M"
+        self.assertEqual(
+            composition_from_poscar_title(title),
+            {
+                "Cr": 3,
+                "Mn": 5,
+                "Fe": 3,
+                "Co": 1,
+                "Ni": 4,
+                "Cu": 4,
+                "Zn": 2,
+                "Y": 5,
+                "Zr": 4,
+                "Nb": 5,
+            },
+        )
+
+
+class TestPotcarPoscarOrderCheck(unittest.TestCase):
+    def _warnings(self, atom_names, atom_numbs, title):
+        with warnings.catch_warnings(record=True) as caught:
+            warnings.simplefilter("always")
+            check_potcar_poscar_order(atom_names, atom_numbs, title)
+        return [str(ii.message) for ii in caught]
+
+    def test_matching_order_is_silent(self):
+        self.assertEqual(self._warnings(["H", "C"], [2, 1], "  H2  C1"), [])
+
+    def test_missing_title_is_silent(self):
+        self.assertEqual(self._warnings(["H", "C"], [2, 1], None), [])
+
+    def test_prose_title_is_silent(self):
+        self.assertEqual(
+            self._warnings(["O", "H"], [1, 2], " POSCAR file written by OVITO"),
+            [],
+        )
+
+    def test_comment_order_does_not_imply_species_order(self):
+        # Pymatgen commonly writes a reduced formula as the free-form comment
+        # even when the actual POSCAR species line uses another order.
+        self.assertEqual(self._warnings(["O", "H"], [1, 2], " H2 O1"), [])
+
+    def test_conflicting_composition_warns(self):
+        messages = self._warnings(
+            ["Li", "O", "F"],
+            [3, 39, 3],
+            " Li3 F39 K3",
+        )
+        self.assertEqual(len(messages), 1)
+        self.assertIn("does not match the explicit composition", messages[0])
+        self.assertIn("Li3 O39 F3", messages[0])
+        self.assertIn("Li3 F39 K3", messages[0])
+
+    def test_truncated_title_compares_only_shared_species(self):
+        # Only the species the title still carries may be compared.
+        title = " Cr3 Mn5 Fe3 Co1 Ni4 Cu4 Zn2 Y5 Zr4 Nb5 M"
+        self.assertEqual(
+            self._warnings(
+                ["Cr", "Mn", "Fe", "Co", "Ni", "Cu", "Zn", "Y", "Zr", "Nb", "Mo"],
+                [3, 5, 3, 1, 4, 4, 2, 5, 4, 5, 7],
+                title,
+            ),
+            [],
+        )
+
+
+class TestOutcarOrderMismatch(unittest.TestCase):
+    """End-to-end: an OUTCAR whose POTCAR order contradicts its POSCAR title."""
+
+    def setUp(self):
+        with open(SOURCE) as fp:
+            self.lines = fp.read().split("\n")
+        self.tmp_dir = tempfile.mkdtemp()
+
+    def tearDown(self):
+        shutil.rmtree(self.tmp_dir, ignore_errors=True)
+
+    def _load(self, lines):
+        path = os.path.join(self.tmp_dir, "OUTCAR")
+        with open(path, "w") as fp:
+            fp.write("\n".join(lines))
+        with warnings.catch_warnings(record=True) as caught:
+            warnings.simplefilter("always")
+            system = dpdata.LabeledSystem(path, fmt="vasp/outcar")
+        messages = [
+            str(ii.message)
+            for ii in caught
+            if "composition produced by pairing" in str(ii.message)
+        ]
+        return system, messages
+
+    def test_consistent_outcar_is_silent(self):
+        system, messages = self._load(self.lines)
+        self.assertEqual(system.data["atom_names"], ["H", "C"])
+        self.assertEqual(messages, [])
+
+    def test_swapped_potcar_order_warns(self):
+        # Swap the two TITEL records so the POTCAR order becomes C, H while
+        # the POSCAR title still says H C.
+        titel = [i for i, ll in enumerate(self.lines) if "TITEL" in ll]
+        self.assertEqual(len(titel), 2)
+        swapped = list(self.lines)
+        swapped[titel[0]], swapped[titel[1]] = (
+            self.lines[titel[1]],
+            self.lines[titel[0]],
+        )
+        title = next(i for i, line in enumerate(swapped) if "POSCAR =" in line)
+        swapped[title] = " POSCAR =   H4 C1"
+        system, messages = self._load(swapped)
+        # dpdata keeps reporting what VASP computed, but says the inputs
+        # disagree.
+        self.assertEqual(system.data["atom_names"], ["C", "H"])
+        self.assertEqual(len(messages), 1)
+        self.assertIn("does not match the explicit composition", messages[0])
+
+    def test_poscar_separator_padding_is_not_part_of_title(self):
+        # VASP commonly emits multiple spaces after ``POSCAR =``. They are
+        # separator padding, not part of the fixed-width title field used to
+        # decide whether the last formula token was truncated.
+        padded = list(self.lines)
+        title = next(i for i, line in enumerate(padded) if "POSCAR =" in line)
+        padded[title] = " POSCAR =   H11111111111111111111111111111111111 O1"
+        _, messages = self._load(padded)
+        self.assertEqual(len(messages), 1)
+        self.assertIn("does not match the explicit composition", messages[0])
+
+
+if __name__ == "__main__":
+    unittest.main()