From a5224ccea74f9c7ec0f223bd92397aa7829ed56c Mon Sep 17 00:00:00 2001 From: aochuba Date: Fri, 23 Jan 2026 19:11:03 +0530 Subject: [PATCH 1/7] Fix load_many crash on single-atom/single-step FCHK optimizations --- iodata/formats/fchk.py | 48 ++++++++++++++++++++----- iodata/test/data/atom_opt.fchk | 52 ++++++++++++++++++++++++++++ iodata/test/test_fchk_single_atom.py | 33 ++++++++++++++++++ 3 files changed, 125 insertions(+), 8 deletions(-) create mode 100644 iodata/test/data/atom_opt.fchk create mode 100644 iodata/test/test_fchk_single_atom.py diff --git a/iodata/formats/fchk.py b/iodata/formats/fchk.py index f1d333c46..0c055ff67 100644 --- a/iodata/formats/fchk.py +++ b/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 *", @@ -338,14 +346,38 @@ def load_many(lit: LineIterator) -> Iterator[dict]: ) # Determine the type of calculation: IRC or Optimization - if "IRC Number of geometries" in fchk: - prefix = "IRC point" - nsteps = fchk["IRC Number of geometries"] - elif "Optimization Number of geometries" in fchk: - prefix = "Opt point" - nsteps = fchk["Optimization Number of geometries"] - else: - raise LoadError("Cannot find IRC or Optimization trajectory in FCHK file.", lit) + try: + if "IRC Number of geometries" in fchk: + prefix = "IRC point" + nsteps = fchk["IRC Number of geometries"] + elif "Optimization Number of geometries" in fchk: + prefix = "Opt point" + nsteps = fchk["Optimization Number of geometries"] + else: + raise LoadError("Cannot find IRC or Optimization trajectory in FCHK file.", lit) + except LoadError: + # 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 natom = fchk["Atomic numbers"].size for ipoint, nstep in enumerate(nsteps): diff --git a/iodata/test/data/atom_opt.fchk b/iodata/test/data/atom_opt.fchk new file mode 100644 index 000000000..33c1c4909 --- /dev/null +++ b/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 diff --git a/iodata/test/test_fchk_single_atom.py b/iodata/test/test_fchk_single_atom.py new file mode 100644 index 000000000..1bf5b1b2c --- /dev/null +++ b/iodata/test/test_fchk_single_atom.py @@ -0,0 +1,33 @@ +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]]) From 357d963ce52afdd170182cdc1f6ff81492220cfc Mon Sep 17 00:00:00 2001 From: aochuba Date: Fri, 23 Jan 2026 21:46:08 +0530 Subject: [PATCH 2/7] fix lint error E501 and TRY301 --- iodata/formats/fchk.py | 27 ++++++++++++--------------- 1 file changed, 12 insertions(+), 15 deletions(-) diff --git a/iodata/formats/fchk.py b/iodata/formats/fchk.py index 0c055ff67..0d3cdb515 100644 --- a/iodata/formats/fchk.py +++ b/iodata/formats/fchk.py @@ -315,10 +315,10 @@ def load_one(lit: LineIterator) -> dict: - ``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. +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. """ @@ -346,16 +346,13 @@ def load_many(lit: LineIterator) -> Iterator[dict]: ) # Determine the type of calculation: IRC or Optimization - try: - if "IRC Number of geometries" in fchk: - prefix = "IRC point" - nsteps = fchk["IRC Number of geometries"] - elif "Optimization Number of geometries" in fchk: - prefix = "Opt point" - nsteps = fchk["Optimization Number of geometries"] - else: - raise LoadError("Cannot find IRC or Optimization trajectory in FCHK file.", lit) - except LoadError: + if "IRC Number of geometries" in fchk: + prefix = "IRC point" + nsteps = fchk["IRC Number of geometries"] + elif "Optimization Number of geometries" in fchk: + 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. @@ -377,7 +374,7 @@ def load_many(lit: LineIterator) -> Iterator[dict]: }, } return - raise + raise LoadError("Cannot find IRC or Optimization trajectory in FCHK file.", lit) natom = fchk["Atomic numbers"].size for ipoint, nstep in enumerate(nsteps): From c6a857d2e30bb4d20750be4bde7b01103681811f Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Fri, 23 Jan 2026 16:18:53 +0000 Subject: [PATCH 3/7] [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci --- iodata/formats/fchk.py | 6 +++--- iodata/test/test_fchk_single_atom.py | 12 ++++++++---- 2 files changed, 11 insertions(+), 7 deletions(-) diff --git a/iodata/formats/fchk.py b/iodata/formats/fchk.py index 0d3cdb515..c821ac723 100644 --- a/iodata/formats/fchk.py +++ b/iodata/formats/fchk.py @@ -363,9 +363,9 @@ def load_many(lit: LineIterator) -> Iterator[dict]: "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), + "atgradient": fchk.get( + "Cartesian Gradient", np.zeros_like(fchk["Current cartesian coordinates"]) + ).reshape(-1, 3), "extra": { "ipoint": 0, "npoint": 1, diff --git a/iodata/test/test_fchk_single_atom.py b/iodata/test/test_fchk_single_atom.py index 1bf5b1b2c..f5484f111 100644 --- a/iodata/test/test_fchk_single_atom.py +++ b/iodata/test/test_fchk_single_atom.py @@ -1,33 +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]]) From ed694557bea1a8898e6a3584fc826f0091f6162d Mon Sep 17 00:00:00 2001 From: aochuba Date: Sun, 25 Jan 2026 03:01:20 +0530 Subject: [PATCH 4/7] Refactor fchk.py to reduce redundancy in load_many --- iodata/formats/fchk.py | 79 +++++++++++++++++++++++++----------------- 1 file changed, 47 insertions(+), 32 deletions(-) diff --git a/iodata/formats/fchk.py b/iodata/formats/fchk.py index 0d3cdb515..7ec6e49ee 100644 --- a/iodata/formats/fchk.py +++ b/iodata/formats/fchk.py @@ -357,22 +357,17 @@ def load_many(lit: LineIterator) -> Iterator[dict]: # (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, - }, - } + 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) @@ -399,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: From 177736bb8e7b88c50b3ff71edac9921f1734f5ba Mon Sep 17 00:00:00 2001 From: aochuba Date: Mon, 26 Jan 2026 00:46:23 +0530 Subject: [PATCH 5/7] Address PR feedback: move test case and add assertions --- iodata/test/test_fchk.py | 33 +++++++++++++++++++++++++++++++++ 1 file changed, 33 insertions(+) diff --git a/iodata/test/test_fchk.py b/iodata/test/test_fchk.py index 1411ae377..7c0fd5514 100644 --- a/iodata/test/test_fchk.py +++ b/iodata/test/test_fchk.py @@ -816,3 +816,36 @@ def test_dump_load_connectivity_peroxide_roundtrip(tmpdir): assert len(mol2.bonds) == len(mol1.bonds) # Bond order should be preserved assert_equal(mol1.bonds[:, 2], mol2.bonds[:, 2]) + +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 + + # Ensure reaction_coordinate is NOT present (Copilot suggestion) + assert "reaction_coordinate" not in mol.extra + + # 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]]) + + # Check gradient (fallback should be zero if missing) (Copilot suggestion) + assert_allclose(mol.atgradient, [[0.0, 0.0, 0.0]]) + From d352053cb0ca3e34f4d48b74677cc52ee40db591 Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Sun, 25 Jan 2026 19:16:51 +0000 Subject: [PATCH 6/7] [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci --- iodata/test/test_fchk.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/iodata/test/test_fchk.py b/iodata/test/test_fchk.py index 7c0fd5514..c5f2dde6f 100644 --- a/iodata/test/test_fchk.py +++ b/iodata/test/test_fchk.py @@ -817,10 +817,11 @@ def test_dump_load_connectivity_peroxide_roundtrip(tmpdir): # Bond order should be preserved assert_equal(mol1.bonds[:, 2], mol2.bonds[:, 2]) + 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] @@ -848,4 +849,3 @@ def test_load_fchk_single_atom_optimization(): # Check gradient (fallback should be zero if missing) (Copilot suggestion) assert_allclose(mol.atgradient, [[0.0, 0.0, 0.0]]) - From 74763cce2ff66af2ec70542735ec2a3be989abd7 Mon Sep 17 00:00:00 2001 From: Aochuba S Aier Date: Sat, 31 Jan 2026 00:58:39 +0530 Subject: [PATCH 7/7] Delete iodata/test/test_fchk_single_atom.py --- iodata/test/test_fchk_single_atom.py | 37 ---------------------------- 1 file changed, 37 deletions(-) delete mode 100644 iodata/test/test_fchk_single_atom.py diff --git a/iodata/test/test_fchk_single_atom.py b/iodata/test/test_fchk_single_atom.py deleted file mode 100644 index f5484f111..000000000 --- a/iodata/test/test_fchk_single_atom.py +++ /dev/null @@ -1,37 +0,0 @@ -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]])