diff --git a/iodata/formats/fchk.py b/iodata/formats/fchk.py index f1d333c46..7ec6e49ee 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 *", @@ -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: 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.py b/iodata/test/test_fchk.py index 1411ae377..c5f2dde6f 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]])