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,27 @@ 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"]:
+            yield {
+                "title": fchk["title"],
+                "atnums": fchk["Atomic numbers"],
+                "atcorenums": fchk["Nuclear charges"],
+                "energy": fchk.get("Total Energy"),
+                "atcoords": fchk["Current cartesian coordinates"].reshape(-1, 3),
+                "atgradient": fchk.get(
+                    "Cartesian Gradient", np.zeros_like(fchk["Current cartesian coordinates"])
+                ).reshape(-1, 3),
+                "extra": {
+                    "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

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]])