Fix load_many crash on single-atom/single-step FCHK optimizations

iodata/formats/fchk.py

@@ -313,6 +313,12 @@ def load_one(lit: LineIterator) -> dict:
 - ``istep`` is the counter within one "point"
 - ``nstep`` is the total number of geometries within in a "point".
 - ``reaction_coordinate`` is only present in case of an IRC calculation.
+
+Load Error Recovery:
+In some cases (e.g. single-step optimizations), the trajectory block (IRC/Optimization Number of
+geometries) may be missing. In these cases, ``load_many`` will attempt to recover by constructing a
+single-frame trajectory from the current geometry and properties (Energy, Gradient) if the file is
+identified as an optimization or scan run.
 """
 
 
@@ -331,6 +337,8 @@ def load_many(lit: LineIterator) -> Iterator[dict]:
             "Atomic numbers",
             "Current cartesian coordinates",
             "Nuclear charges",
+            "Total Energy",
+            "Cartesian Gradient",
             "IRC *",
             "Optimization *",
             "Opt point *",
@@ -345,6 +353,22 @@ def load_many(lit: LineIterator) -> Iterator[dict]:
         prefix = "Opt point"
         nsteps = fchk["Optimization Number of geometries"]
     else:
+        # If the trajectory is missing, we might have a single-frame optimization
+        # (e.g. because it converged in one step). In that case, we return the
+        # current geometry as the only frame.
+        if fchk.get("command") in ["FOpt", "Scan", "IRC"]:
+            atcoords = fchk["Current cartesian coordinates"].reshape(-1, 3)
+            atgradient = fchk.get(
+                "Cartesian Gradient", np.zeros_like(fchk["Current cartesian coordinates"])
+            ).reshape(-1, 3)
+            yield _create_frame(
+                fchk,
+                fchk.get("Total Energy"),
+                atcoords,
+                atgradient,
+                {"ipoint": 0, "npoint": 1, "istep": 0, "nstep": 1},
+            )
+            return
         raise LoadError("Cannot find IRC or Optimization trajectory in FCHK file.", lit)
 
     natom = fchk["Atomic numbers"].size
@@ -370,23 +394,43 @@ def load_many(lit: LineIterator) -> Iterator[dict]:
             )
         reference_energy = fchk.get("Optimization Reference Energy", 0.0)
         for istep, (energy, recor, atcoords, gradients) in enumerate(trajectory):
-            data = {
-                "title": fchk["title"],
-                "atnums": fchk["Atomic numbers"],
-                "atcorenums": fchk["Nuclear charges"],
-                "energy": energy + reference_energy,
-                "atcoords": atcoords,
-                "atgradient": gradients,
-                "extra": {
-                    "ipoint": ipoint,
-                    "npoint": len(nsteps),
-                    "istep": istep,
-                    "nstep": len(trajectory),
-                },
+            extra = {
+                "ipoint": ipoint,
+                "npoint": len(nsteps),
+                "istep": istep,
+                "nstep": len(trajectory),
             }
-            if prefix == "IRC point":
-                data["extra"]["reaction_coordinate"] = recor
-            yield data
+            yield _create_frame(
+                fchk,
+                energy + reference_energy,
+                atcoords,
+                gradients,
+                extra,
+                recor if prefix == "IRC point" else None,
+            )
+
+
+def _create_frame(
+    fchk: dict,
+    energy: float | None,
+    atcoords: NDArray[float],
+    atgradient: NDArray[float],
+    extra: dict,
+    reaction_coordinate: float | None = None,
+) -> dict:
+    """Create a frame dictionary from FCHK data."""
+    data = {
+        "title": fchk["title"],
+        "atnums": fchk["Atomic numbers"],
+        "atcorenums": fchk["Nuclear charges"],
+        "energy": energy,
+        "atcoords": atcoords,
+        "atgradient": atgradient,
+        "extra": extra,
+    }
+    if reaction_coordinate is not None:
+        data["extra"]["reaction_coordinate"] = reaction_coordinate
+    return data
 
 
 def _load_fchk_low(lit: LineIterator, label_patterns: list[str] | None = None) -> dict:

iodata/test/data/atom_opt.fchk

@@ -0,0 +1,52 @@
+h_sto3g
+FOpt      UHF                                                         STO-3G
+Number of atoms                            I                1
+Charge                                     I                0
+Multiplicity                               I                2
+Number of electrons                        I                1
+Number of alpha electrons                  I                1
+Number of beta electrons                   I                0
+Number of basis functions                  I                1
+Number of contracted shells                I                1
+Highest angular momentum                   I                0
+Largest degree of contraction              I                3
+Number of primitive shells                 I                3
+Virial Ratio                               R      1.613897736040718E+00
+SCF Energy                                 R     -4.665818503844346E-01
+Total Energy                               R     -4.665818503844346E-01
+Atomic numbers                             I   N=           1
+           1
+Nuclear charges                            R   N=           1
+  1.00000000E+00
+Current cartesian coordinates              R   N=           3
+  0.00000000E+00  0.00000000E+00  0.00000000E+00
+Integer atomic weights                     I   N=           1
+           1
+Real atomic weights                        R   N=           1
+  1.00782504E+00
+Shell types                                I   N=           1
+           0
+Number of primitives per shell             I   N=           1
+           3
+Shell to atom map                          I   N=           1
+           1
+Primitive exponents                        R   N=           3
+  3.42525091E+00  6.23913730E-01  1.68855404E-01
+Contraction coefficients                   R   N=           3
+  1.54328967E-01  5.35328142E-01  4.44634542E-01
+Coordinates of each shell                  R   N=           3
+  0.00000000E+00  0.00000000E+00  0.00000000E+00
+Alpha Orbital Energies                     R   N=           1
+ -4.66581850E-01
+Beta Orbital Energies                      R   N=           1
+  3.08024094E-01
+Alpha MO coefficients                      R   N=           1
+  1.00000000E+00
+Beta MO coefficients                       R   N=           1
+  1.00000000E+00
+Total SCF Density                          R   N=           1
+  1.00000000E+00
+Spin SCF Density                           R   N=           1
+  1.00000000E+00
+Dipole Moment                              R   N=           3
+  0.00000000E+00  0.00000000E+00  0.00000000E+00

iodata/test/test_fchk_single_atom.py

@@ -0,0 +1,37 @@
+from importlib.resources import as_file, files
+
+from numpy.testing import assert_allclose
+
+from ..api import load_many
+
+
+def load_fchk_trj_helper(fn_fchk):
+    """Load a trajectory from a testing fchk file with iodata.iodata.load_many."""
+    with as_file(files("iodata.test.data").joinpath(fn_fchk)) as fn:
+        return list(load_many(fn))
+
+
+def test_load_fchk_single_atom_optimization():
+    """Test loading a single-atom optimization (missing trajectory block)."""
+    trj = load_fchk_trj_helper("atom_opt.fchk")
+    # Should load exactly one frame
+    assert len(trj) == 1
+    mol = trj[0]
+
+    # Check basic properties
+    assert mol.natom == 1
+    assert mol.atnums[0] == 1
+    assert mol.atcorenums[0] == 1.0
+
+    # Check fallback values
+    assert mol.extra["ipoint"] == 0
+    assert mol.extra["npoint"] == 1
+    assert mol.extra["istep"] == 0
+    assert mol.extra["nstep"] == 1
+
+    # Check that properties were correctly loaded from the single frame
+    assert mol.energy is not None
+    assert_allclose(mol.energy, -4.665818503844346e-01)
+
+    # Check coordinates
+    assert_allclose(mol.atcoords, [[0.0, 0.0, 0.0]])