Adds output of Hessians to xtb when acting as an external program for Gaussian 16

[clavigne]

This is a feature improvement related to #329 .

Basically, I've added Hessian output capability for xtb acting as an external program invoked by Gaussian 16. I've made the changes to writeResultsGaussianExternal that are necesserary for the Hessian output.

This would have been otherwise extremely straightforward, except for the fact that the raw Hessian matrix (before mass-weighting, which is what Gaussian expects) is not saved in numhess for (I assume) memory usage reasons. I didn't want to just have the code keep the raw Hessian every time for the vast majority of users not running under Gaussian, so I short-circuited numhess instead to immediately return when the input file is a Gaussian external file. (This additionally has the advantage that it avoids an extraneous Hessian diagonalization by xtb, which gets expensive for large molecules.)

This solution is far from elegant as it required adding some more global state. Feel free to suggest an alternative.

[clavigne]

In addition, I'd be very happy if there was an "easy" way to control runtype based on the mode flag in

xtb/src/io/reader/gaussian.f90

Line 47 in 7e8e21c

read(unit, '(4i10)', iostat=err) n, mode, chrg, spin

which is set by Gaussian to 1 for gradient and 2 for Hessian evaluation. I've considered calling set_runtyp in the reading subroutine but that's maybe a bit too much global state manipulation. The alternative I can see is to create an entirely new "under gaussian control" runtype.

Very interested in your opinion on this, thanks.

[codecov]

Codecov Report

Merging #363 (629fc59) into master (2dc7fba) will decrease coverage by 0.03%.
The diff coverage is 3.22%.

@@            Coverage Diff             @@
##           master     #363      +/-   ##
==========================================
- Coverage   40.80%   40.77%   -0.04%     
==========================================
  Files         304      304              
  Lines       50774    50813      +39     
==========================================
  Hits        20718    20718              
- Misses      30056    30095      +39     

Impacted Files	Coverage Δ
src/hessian.F90	0.00% <0.00%> (ø)
src/io/reader/gaussian.f90	0.00% <0.00%> (ø)
src/io/writer/gaussian.f90	0.00% <0.00%> (ø)
src/set_module.f90	27.82% <0.00%> (-0.03%)	⬇️
src/setparam.f90	23.80% <ø> (ø)
src/prog/main.F90	40.95% <20.00%> (-0.13%)	⬇️
src/relaxation_engine.f90	0.00% <0.00%> (ø)

---

Continue to review full report at Codecov.

Legend - Click here to learn more
Δ = absolute <relative> (impact), ø = not affected, ? = missing data
Powered by Codecov. Last update 2dc7fba...629fc59. Read the comment docs.

[awvwgk]

The read mode could be cached in the molecular structure type and the main routine can than decide if it is able to invoke set_runtyp. There is no guarantee that this will work if the user already provided a flag, which locks it. A safer way would be to allow the xtb-gaussian wrapper to query the first line of the input and set the flags accordingly:

flag=$(awk 'NR == 1 {if ($2 == 2) {print "--hess"} else {print "--grad"} }' $2)

[clavigne]

Thanks for this and you are right. awk is the right tool to use here and will allow me to remove any final Perl dependencies, moving the whole thing to a self contained two line Bash script.

Will update (and sign off commits! Sorry again)

Should I add a unit test too? Without Gaussian of course, just testing file io structure.

[awvwgk]

We don't have unit tests for the frequency calculation yet, since it is a bit costly in the CI. I'll have a look myself while checking the PR.

[clavigne]

I made all the suggested changes. The script with awk works beautifully, and the performance is much better than by passing through perl or python to do the required Hessian reformatting.

Thanks for your help!

[awvwgk]

Overall it still looks a bit hacky, but as long as it only triggers with Gaussian input we can be sure it won't mess up regular xtb calculations. I added the dipole gradient to the printout as well, storing it in the intensities, which is somewhat inline with the hessian being stored in the normal modes.

[clavigne]

I agree it is pretty hacky. I don't need the feature right now, as I can just compile my PR branch myself, so if you think it's best to refactor instead, I'm open to it.