✨ PwOutput: Parse forces from stdout as fallback

src/qe_tools/outputs/parsers/pw.py

@@ -75,6 +75,22 @@ def parse(content):
     r"\s*wk\s*=\s*([\-\d.E+]+)"
 )
 
+# Captures the `Forces acting on atoms (cartesian axes, Ry/au):` header plus the
+# following contiguous block of `atom N type M   force =  Fx  Fy  Fz` rows. QE
+# follows this block with `The non-local contrib. to forces`, `The ionic
+# contribution to forces`, and `The local contribution to forces` blocks that
+# reuse the row format -- anchoring on the `Forces acting on atoms` header (and
+# stopping at the first non-row line) ensures only the total forces are captured.
+_FORCES_BLOCK_RE = re.compile(
+    r"Forces acting on atoms\s*\(cartesian axes,\s*Ry/au\):\s*\n\s*\n"
+    r"(?P<rows>(?:\s*atom\s+\d+\s+type\s+\d+\s+force\s*=\s*"
+    r"[\-\d.E+]+\s+[\-\d.E+]+\s+[\-\d.E+]+\s*\n)+)"
+)
+_FORCES_ROW_RE = re.compile(
+    r"atom\s+\d+\s+type\s+\d+\s+force\s*=\s*"
+    r"([\-\d.E+]+)\s+([\-\d.E+]+)\s+([\-\d.E+]+)"
+)
+
 
 class PwStdoutParser(BaseStdoutParser):
     """
@@ -100,6 +116,14 @@ def parse(content: str) -> dict:
         if energy_matches:
             parsed_data["total_energy"] = float(energy_matches[-1])
 
+        # Forces: take the last block, since relax/md runs print one per ionic step.
+        force_blocks = list(_FORCES_BLOCK_RE.finditer(content))
+        if force_blocks:
+            parsed_data["forces"] = [
+                [float(v) for v in row.groups()]
+                for row in _FORCES_ROW_RE.finditer(force_blocks[-1].group("rows"))
+            ]
+
         # k-points: take the last block, since vc-relax reprints after relaxation.
         # The stdout values are in 2pi/alat; convert to 1/Å using `alat` (also from
         # stdout, in bohr) so downstream Specs match the XML-derived units.

src/qe_tools/outputs/pw.py

@@ -197,15 +197,29 @@ class _PwMapping:
     forces: Annotated[
         list,
         Spec(
-            (
-                "xml.output.forces",
-                lambda forces: [
-                    [
-                        value * CONSTANTS.hartree_to_ev / CONSTANTS.bohr_to_ang
-                        for value in forces["$"][atom_index * 3 : (atom_index + 1) * 3]
-                    ]
-                    for atom_index in range(forces["@dims"][1])
-                ],
+            Coalesce(
+                (
+                    "xml.output.forces",
+                    lambda forces: [
+                        [
+                            value * CONSTANTS.hartree_to_ev / CONSTANTS.bohr_to_ang
+                            for value in forces["$"][
+                                atom_index * 3 : (atom_index + 1) * 3
+                            ]
+                        ]
+                        for atom_index in range(forces["@dims"][1])
+                    ],
+                ),
+                (
+                    "stdout.forces",
+                    lambda forces: [
+                        [
+                            value * CONSTANTS.ry_to_ev / CONSTANTS.bohr_to_ang
+                            for value in row
+                        ]
+                        for row in forces
+                    ],
+                ),
             )
         ),
         Unit("eV/angstrom"),

tests/outputs/test_pw_output.py

@@ -100,6 +100,35 @@ def test_total_energy_nscf_clobber():
     )
 
 
+def test_forces_stdout_xml_agree():
+    """XML and stdout `forces` agree within stdout's printed precision.
+
+    QE prints forces in `cartesian axes, Ry/au` to 8 decimal digits (~1e-8 Ry/bohr).
+    `xml.output.forces` is in Hartree/bohr; multiplying the stdout values by 2
+    (1 Ha = 2 Ry) puts both on the same scale.
+    """
+    pw_directory = Path(__file__).parent / "fixtures" / "pw" / "default_xml_220603"
+
+    pw_out = PwOutput.from_dir(pw_directory)
+    xml_forces = pw_out.raw_outputs["xml"]["output"]["forces"]["$"]
+    stdout_forces = np.array(pw_out.raw_outputs["stdout"]["forces"]).flatten()
+
+    np.testing.assert_allclose(2 * np.asarray(xml_forces), stdout_forces, atol=1e-8)
+
+
+def test_forces_stdout_fallback_without_xml():
+    """`PwOutput.from_files(stdout=...)` exposes `forces` via the stdout fallback."""
+    stdout_file = (
+        Path(__file__).parent / "fixtures" / "pw" / "default_xml_220603" / "pw.out"
+    )
+
+    pw_out = PwOutput.from_files(stdout=stdout_file)
+
+    assert "xml" not in pw_out.raw_outputs
+    assert pw_out.get_output("forces") is not None
+    assert len(pw_out.get_output("forces")) == 2
+
+
 def test_k_points_stdout_xml_agree():
     """Stdout-derived k-points agree with the XML values within stdout's printed precision.
 

tests/outputs/test_pw_output/test_default_xml_211101_.yml

@@ -101,6 +101,13 @@ base_outputs:
 raw_outputs:
   stdout:
     code_version: '7.0'
+    forces:
+    - - 0.0
+      - 0.0
+      - -0.00666975
+    - - 0.0
+      - 0.0
+      - 0.00666975
     k_points_cartesian:
     - - 0.0
       - 0.0

tests/outputs/test_pw_output/test_default_xml_220603_.yml

@@ -117,6 +117,13 @@ base_outputs:
 raw_outputs:
   stdout:
     code_version: '7.1'
+    forces:
+    - - 0.0
+      - 0.0
+      - -1.4e-06
+    - - 0.0
+      - -0.0
+      - 1.4e-06
     k_points_cartesian:
     - - 0.0
       - 0.0

tests/outputs/test_pw_output/test_default_xml_230310_.yml

@@ -168,6 +168,19 @@ base_outputs:
 raw_outputs:
   stdout:
     code_version: '7.2'
+    forces:
+    - - 0.0
+      - 0.0
+      - 0.0
+    - - 0.0
+      - 0.0
+      - -0.02705325
+    - - 0.0
+      - 0.0
+      - 0.02705325
+    - - 0.0
+      - 0.0
+      - -0.0
     highest_occupied_level: 9.0922
     k_points_cartesian:
     - - 0.0

tests/outputs/test_pw_output/test_default_xml_240411_.yml

@@ -153,6 +153,13 @@ base_outputs:
 raw_outputs:
   stdout:
     code_version: 7.3.1
+    forces:
+    - - 0.0
+      - 0.0
+      - -2.3e-07
+    - - -0.0
+      - -0.0
+      - 2.3e-07
     k_points_cartesian:
     - - 0.0
       - 0.0

tests/outputs/test_pw_output/test_default_xml_250521_.yml

@@ -152,6 +152,31 @@ base_outputs:
 raw_outputs:
   stdout:
     code_version: '7.5'
+    forces:
+    - - -0.06858972
+      - 0.01765138
+      - 0.01964435
+    - - 0.0388041
+      - 0.03512021
+      - 0.02191125
+    - - 0.03004469
+      - 0.01294412
+      - -0.01093989
+    - - -0.00017394
+      - -0.04589712
+      - -0.02011802
+    - - -0.02178025
+      - -0.02664816
+      - -0.01644639
+    - - 0.01784844
+      - 0.02551909
+      - 0.00949881
+    - - -0.00075029
+      - 0.00122727
+      - -0.00636038
+    - - 0.00459697
+      - -0.01991679
+      - 0.00281027
     k_points_cartesian:
     - - 0.0
       - 0.0