Quasi-RRHO Thermochemistry Analysis Module

[arepstein]

Summary

Added quasirrho.py to the analysis subpackage to calculate the Quasi-RRHO free energy from a Gaussian or QChem frequency calculation.

• Calculates Grimme's Quasi-RRHO free energy
• Option for also correcting for solvent concentration

Additional dependencies introduced (if any)

*None

TODO (if any)

• The rotational symmetry number is required as an input. To calculate this on the fly, will consider updates to PointGroupAnalyzer

Before a pull request can be merged, the following items must be checked:

• [ x] Code is in the standard Python style. Run pycodestyle and flake8 on your local machine.
• [x ] Docstrings have been added in the Google docstring format. Run pydocstyle on your code.
• [ ] Type annotations are highly encouraged. Run mypy to type check your code.
• [x ] Tests have been added for any new functionality or bug fixes.
• [ x] All existing tests pass.

Note that the CI system will run all the above checks. But it will be much more efficient if you already fix most errors prior to submitting the PR. It is highly recommended that you use the pre-commit hook provided in the pymatgen repository. Simply cp pre-commit .git/hooks and a check will be run prior to allowing commits.

[mkhorton]

@arepstein is this PR ready for review to be merged?

[arepstein]

I believe so. The only pending update in my mind is automated detection of the rotational symmetry number, but I think it's okay as an input parameter for now.

[htz1992213]

Hey @arepstein , I get to know this PR in today's subgroup meeting. I am just wondering how the functionalities here comparing to those in the Goodvibes package https://github.com/bobbypaton/GoodVibes?

[arepstein]

Hey @arepstein , I get to know this PR in today's subgroup meeting. I am just wondering how the functionalities here comparing to those in the Goodvibes package https://github.com/bobbypaton/GoodVibes?

Thanks for pointing this out, I didn't know about Goodvibes when I put in this PR. This PR should be identical to the Grimme Quasi-RRHO approximation for the entropy that's implemented in GoodVibes, but does not include any other methods that GoodVibes implements. One big difference I see in terms of implementation into pymatgen infrastructure is that this PR can calculate Quasi-RRHO entropeis for Q-Chem output files, Gaussian output files, or manual input parameters, which is useful for Atomate integration. It would be a good idea to check that this agrees with GoodVibes.

[coveralls]

Coverage decreased (-0.7%) to 83.428% when pulling d08a1ca5c4a2d5e3f29703767642b4bfa5fd7ef0 on arepstein:readytoPR into 3376f27a9c5a4cac1590169f2212361ab6c94388 on materialsproject:master.

[rkingsbury]

See also some small edits I made when I used the code, in a PR against your fork: https://github.com/arepstein/pymatgen/pull/1

[rkingsbury]

I will also note for other reviewers that although the GoodVibes package contains similar functionality, it only accepts Gaussian output files, and hence is not useful with our high-throughput Q-Chem infrastructure. So this PR is valuable because it gives us those capabilities.

[janosh]

@arepstein Is this still WIP?

[arepstein]

@janosh Yes, still WIP, I haven't incorporated Ryan's comments yet. Thanks @rkingsbury for the comments!

[rkingsbury]

@janosh Yes, still WIP, I haven't incorporated Ryan's comments yet. Thanks @rkingsbury for the comments!

No problem! Recall that I also opened a small PR against your branch with some other edits, notably exposing h_corrected as an attribute.

I had another thought about the call signature which I think would be better than what I initially proposed in that PR. Currently you have

def __init__(self, output: GaussianOutput | QCOutput | dict, sigma_r=1, temp=298.15, press=101317, conc=1, v0=100):

I think it would be cleaner to make the regular __init__ method generic by using explicit kwargs, e.g.

def __init__(self, mol: Molecule, frequencies: List, energy: float, sigma_r=1, temp=298.15, press=101317, conc=1, v0=100):

Where mol, frequencies, and energy are the items currently extracted from the GaussianOutput / QCOutput / dict. Then, you would add two @classmethod: .from_gaussian(self, output: GaussianOutput) and .from_qchem(self, output: QCOutput) to allow instantiating from the respective classes.

This structure allows for 1) use of the module with output from arbitrary codes (i.e., any molecule and frequency combo) and 2) leaves the door open for additional .from_xxxx methods to cleanly expand support to new codes in the future. You could even add some error handling that ensure the number of frequencies is consistent with the size of the molecule

[samblau]

Funny timing - I now would very much like to use this capability! One thing I want to note is that it would be very nice if, in addition to an actual Q-Chem output file, we could pass it a summary document that comes out of the new Molecules Emmet infrastructure. Happy to provide an example JSON.

[arepstein]

Hey @samblau , I'm getting to this now. Could you please provide an example JSON? Right now it can take a manual dictionary input, so it should be easy to incorporate.

[arepstein]

Hey @rkingsbury and @janosh, I've incorporated Ryan's edits and recommendations. While making these edits it became clear that a Class might not be the best choice for QuasiRRHO. As it stands now, it might be better as a simple function. Ideally, I think it could be nice to have a Thermochemistry class that is a Molecule plus some extra information like frequencies. We could then have different treatments of thermochemistry, like GoodVibes implements, but in a way works well with our ecosystem.

That being said, hopefully QuasiRRHO is still useful for now and can be a starting point for later changes if desired.

[codecov-commenter]

Codecov Report

:exclamation: No coverage uploaded for pull request base (master@57d8a2f). Click here to learn what that means. Patch has no changes to coverable lines.

Additional details and impacted files

@@            Coverage Diff            @@
##             master    #2028   +/-   ##
=========================================
  Coverage          ?   74.06%           
=========================================
  Files             ?      230           
  Lines             ?    69403           
  Branches          ?    16161           
=========================================
  Hits              ?    51403           
  Misses            ?    14957           
  Partials          ?     3043           

:umbrella: View full report in Codecov by Sentry.
:loudspeaker: Have feedback on the report? Share it here.

[rkingsbury]

I think this is ready to merge, right @arepstein ? @janosh can you remove the 'stale' label?

[arepstein]

Yes, I believe ready to merge!

[arepstein]

Thanks @janosh and @rkingsbury for all the help, edits, and revisions!