duartegroup / autodE

automated reaction profile generation
https://duartegroup.github.io/autodE/
MIT License
173 stars 52 forks source link

Extract coordinates failed #358

Closed javialra97 closed 4 weeks ago

javialra97 commented 1 month ago

Describe the bug Sometimes the Gaussian doesn't print the molecule's geometry as 'Input orientation' and only prints the 'Standard orientation'. In these cases, autodE is not able to extract the geometry (L549, autode/wrappers/G09, func coordinates_from) and returns an AssertionError. I have seen these errors mainly during the initial frequency calculation, before the optimization of a TS guess.

To Reproduce Script to launch it:

import autode as ade
if __name__ == "__main__":
    ade.Config.n_cores=32
    ade.Config.max_core=1000
    ade.Config.hcode="G16"
    ade.Config.lcode ="xtb"
    rxn=ade.Reaction(r"C1#Cc2ccccc2CCc2ccccc21.COC(=O)c1nnc(C(C)(C)C)nn1>>COC(=O)[C@]12N=N[C@](C(C)(C)C)(N=N1)C1=C2c2ccccc2CCc2ccccc21")
    ade.Config.G16.keywords.set_functional('cam-b3lyp')
    ade.Config.G16.keywords.opt.basis_set = '6-311++G**'
    ade.Config.G16.keywords.opt_ts.basis_set = '6-311++G**'
    ade.Config.G16.keywords.hess.basis_set = '6-311++G**'
    ade.Config.G16.keywords.low_opt.basis_set = '6-31G*'
    ade.Config.G16.keywords.low_opt.max_opt_cycles = 15
    ade.Config.num_conformers=1000
    ade.Config.rmsd_threshold=0.1
    ade.Config.hmethod_conformers=True
    rxn.calculate_reaction_profile(free_energy=True)

Error:

File "/home/javialra/anaconda3/envs/autode_original/lib/python3.10/site-packages/autode/wrappers/G09.py", line 675, in hessian_from atoms=self.atoms_from(calc), File "/home/javialra/anaconda3/envs/autode_original/lib/python3.10/site-packages/autode/wrappers/methods.py", line 250, in atoms_from atoms.coordinates = self.coordinates_from(calc) File "/home/javialra/anaconda3/envs/autode_original/lib/python3.10/site-packages/autode/atoms.py", line 641, in coordinates assert value.shape == (len(self), 3) AssertionError

If it is necessary I can provide the full output.

Expected behavior Successful termination

Environment

Additional context To work around, I added the 'Standard orientation' in the if statement (L549, autode/wrappers/G09, func coordinates_from) and then a break at the end.

            if "Input orientation" in line or "Standard orientation" in line:
                coords.clear()
                xyz_lines = calc.output.file_lines[
                    i + 5 : i + 5 + calc.molecule.n_atoms
                ]

                for xyz_line in xyz_lines:
                    _, _, _, x, y, z = xyz_line.split()
                    coords.append([float(x), float(y), float(z)])
                break

PS: Amazing work with autodE!!

t-young31 commented 1 month ago

Hi @javialra97 – thanks a lot for reporting this 😄

Is the break necessary? As with it I think coordinate extraction from optimisation calculations won't grab the last (i.e. optimised) set of coordinates?

javialra97 commented 1 month ago

Hi @t-young31

True, it is not necessary. I put the break because I was only associating it with the frequency calculation, but yes, you won't grab the last set of coordinates. Also, I thought that you extracted the geometry from the final Gaussian block.