load_many from fchk fails if the Gaussian job is an optimization of a single atom

[mbocus]

Hi all,

I have stumbled upon an issue when loading a formatted checkpoint file of Gaussian (*.fchk). I have a script that checks the job type and, if it is an optimization, uses the load_many function to parse the output. If the optimization of a species with more than one atom finishes in one step (i.e., the structure is already optimized), all works as intended. For example:

for mol in load_many("water.fchk"):
    print(mol.atcoords)

returns

[[-3.13242678e-34 -7.31649715e-35  2.40261291e-01]
 [-2.33989326e-32  1.43244075e+00 -9.61045165e-01]
 [-1.75423398e-16 -1.43244075e+00 -9.61045165e-01]]

On the other hand, in the limiting case in which a single atom is "optimized" (which is obviously meaningless, but something that can happen in automated workflows), the reading fails with the following error:

Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "/home/mbocus/miniconda3/envs/iodata_env/lib/python3.11/site-packages/iodata/api.py", line 219, in load_many
    for data in format_module.load_many(lit, **kwargs):
  File "/home/mbocus/miniconda3/envs/iodata_env/lib/python3.11/site-packages/iodata/formats/fchk.py", line 329, in load_many
    raise LoadError("Cannot find IRC or Optimization trajectory in FCHK file.", lit)
iodata.utils.LoadError: Cannot find IRC or Optimization trajectory in FCHK file. (atom.fchk:1044)

I am using the latest iodata release (v1.0.0a8) and, attached, you may find the two fchk files I used for testing.

Thanks a lot in advance,
Massimo

water.fchk.txt
atom.fchk.txt

[tovrstra]

Thanks for bringing this up. I guess the check for the IRC is too strict and should be changed so that the IRC is loaded when available and skipped otherwise.

[Ao-chuba]

Hi, @tovrstra Gaussian omits the Optimization Number of geometries block for single step optimizations, even though the file is still marked as an optimization run (FOpt). This causes load_many to fail early.
If the trajectory block lookup fails (LoadError), we shouldn’t give up immediately. Instead, we check whether the file is marked as an optimization run. If it is, we can safely treat it as a trajectory of length 1 and construct a single frame using the standard Current Cartesian Coordinates, Total Energy, and Cartesian Gradient fields, which are already present.

This recovers valid data for single step optimizations without affecting the behavior for standard multi-step optimization files.
If this approach sounds good to you, could you please assign this issue to me? I’d be happy to open a PR with the fix.

[tovrstra]

Sorry for the slow response. I'm a bit snowed under with exams. GH issues are on the back burner a.t.m.

[Ao-chuba]

@tovrstra No worries sir, I’m glad to be contributing here. Thank you for the reply as well.
I’ve also created a PR that refactors iodata/formats/fchk.py to remove the redundancy in load_many, as mentioned in your request during the review of my PR.