t-young31
commented
3 years ago Hi!
So by default a Hessian evaluation is not performed at the end of an optimisation, as it's usually pretty expensive to calculate. Trying to get the Hessian from a calculation where one hasn't been calculated is going to be tricky!
This seems to work fine:
>>> import autode
>>> h2 = autode.Molecule(smiles='[H][H]')
>>> g16 = autode.methods.G16()
>>> calc = autode.Calculation(name='h2',molecule=h2,method=g16,keywords=g16.keywords.hess)
>>> calc.run()
>>> calc.get_hessian()
Hessian([[ 1.84704716 0. 0. -1.84704716 0. 0. ]
[ 0. -0.12221286 0. 0. 0.12221286 0. ]
[ 0. 0. -0.12221286 0. 0. 0.12221286]
[-1.84704716 0. 0. 1.84704716 0. 0. ]
[ 0. 0.12221286 0. 0. -0.12221286 0. ]
[ 0. 0. 0.12221286 0. 0. -0.12221286]] Ha Å^-2)Perhaps there should be a warning if a Hessian doesn't exist.
Describe the bug
I. When a Hessian matrix is to be calculated, it does not work. The error message is CouldNotGetProperty: Could not get the Hessian: could not convert string to float: '1\H'. I tried with other SMILEs string and they still do no work.
II. The bottom line block (i.e. 1\1\ ... @\n) of log file generated after running Gaussian16 has "3 * number_of_atoms". In this case 6. Due to this, ValueError: Array was not the correct shape to be broadcast into the lower triangle. Need N(N+1)/2 elements for an NxN array is made. Is the bottom block really hessian information? (I am not versed in QM calculation).
To Reproduce
import autode h2 = autode.Molecule(smiles='[H][H]') calc=autode.Calculation(name='h2',molecule=h2,method=autode.wrappers.G16.G16(),keywords=autode.wrappers.G16.G16().keywords.opt) calc.run() calc.get_energy() #Works! calc.get_gradients() #Works! calc.get_hessian() #Breaks with msg: CouldNotGetProperty: Could not get the Hessian: could not convert string to float: '1\H'
Expected behavior
Hessian matrix should be calculated.
Environment
Suggestion I. In autode/wrappers/G09.py file,
CHANGE from hess_str = "".join(hess_lines[::-1]).split(r'\')[-3] hess_values = [float(val) for val in hess_str.split(',')] TO hess_str = list(map(lambda inp: inp.split("\")[0], hess_str.split(',')[2:])) #CouldNotGetProperty error will be fixed as a result of removing \H, \C character strings, in addition to removing '1' and '1' characters. hess_values = [float(val) for val in hess_str.split(',')] #This will save values correctly into a float value list.
II. Is the bottom line block (i.e. 1\1\ ... @\n) really Hessian? In this case 6 elements (3 * num_atoms) exist. However, this information seems to be exactly the same as optimized Cartesian coordinate (not Hessian!), accessible by calc.get_final_atoms(). If so, this should be corrected and correct method to calculate (or read) Hessian matrix should be updated.
Additional context