Refactor energy corrections in Arkane and implement new bond additivity correction type by cgrambow · Pull Request #1605 · ReactionMechanismGenerator/RMG-Py

cgrambow · 2019-05-22T22:54:38Z

Motivation or Problem

Applying energy corrections as part of a statmech job in Arkane takes place in a really messy function and does not allow for easy extension to new bond additivity correction types.

Description of Changes

This PR refactors energy corrections into separate modules and stores the data separately. Furthermore, it implements a new type of bond additivity correction scheme.

Testing

Existing unit tests pass.

Discussion

This PR provides the functionality to select either type of bond additivity correction, but it doesn't yet provide data for the new BAC type. My plan is to completely replace the current atom energies using just the calculated values at the respective levels of theory and to provide newly-derived bond additivity corrections for all levels of theory for both the existing type of BACs and the newly implemented type of BACs once we have selected a reference set of molecules and performed the quantum calculations at all the levels of theory.

This PR can be merged without the new data because it doesn't change current Arkane behavior. Alternatively, we can wait with merging until the new data is available. What do you think? In either case, the current changes are ready for review.

codecov · 2019-05-23T00:24:12Z

Codecov Report

Merging #1605 into master will decrease coverage by 0.13%.
The diff coverage is n/a.

@@            Coverage Diff             @@
##           master    #1605      +/-   ##
==========================================
- Coverage    41.8%   41.67%   -0.14%     
==========================================
  Files         177      176       -1     
  Lines       29455    29374      -81     
  Branches     6059     6059              
==========================================
- Hits        12314    12241      -73     
+ Misses      16271    16263       -8     
  Partials      870      870

Impacted Files	Coverage Δ
rmgpy/exceptions.py

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 f6faa9f...cdca610. Read the comment docs.

alongd

Thanks for this awesome PR. I really like the organization and methodology. I added relatively minor comments.

alongd · 2019-05-23T13:19:27Z

arkane/encorr/corr.py

+        try:
+            atom_energies = data.atom_energies[model_chemistry]
+        except KeyError:
+            raise Exception('Missing atom energies for model chemistry {}'.format(model_chemistry))


Could you make the exceptions (throughout the PR) more specific? For example, if we raise an AtomEnergyCorrectionException here, then a software like ARC calling Akane might know what the problem is and in this case be able to spawn AEC jobs for that model chemistry.

alongd · 2019-05-23T13:23:13Z

arkane/encorr/corr.py

+
+    for symbol, count in atoms.items():
+        if symbol in atom_energies:
+            corr -= count * atom_energies[symbol] * 4.35974394e-18 * constants.Na


Could you add a comment re the unit conversion factor (e.g., Hartree to J)?
Alternatively, and although it's not a "constant", would we like to add frequently used unit conversions to constants.py?

I added comments everywhere but just kept the numerical values in there for now. I'm inclined to not add unit conversions to constants.py, but am open to other suggestions.

alongd · 2019-05-23T13:25:10Z

arkane/encorr/corr.py

+            corr -= count * atom_energies[symbol] * 4.35974394e-18 * constants.Na
+        else:
+            raise Exception(
+                'Unknown element "{}". Turn off atom corrections if only running a kinetics jobs '


Instead of "unknown element", perhaps say "An energy correction for element X is unavailable in the data for model chemistry Y. Turn off..."

alongd · 2019-05-23T13:30:26Z

arkane/encorr/data.py

+
+# Atom energy corrections to reach gas-phase reference state
+# Experimental enthalpy of formation at 0 K, 1 bar for gas phase
+# See Gaussian thermo whitepaper at http://www.gaussian.com/g_whitepap/thermo.htm)


remove ) at the end of the line

Turns out the link has actually changed, so I replaced it with a more up-to-date link.

alongd · 2019-05-23T13:33:00Z

arkane/encorr/data.py

+           'Rb': 17.86, 'Ag': 66.61, 'Cd': 25.240, 'Sn': 70.50, 'I': 24.04, 'Xe': -1.481,
+           'Cs': 16.80, 'Hg': 13.19, 'Pb': 15.17}
+
+# Thermal contribution to enthalpy Hss(298 K) - Hss(0 K) reported by Gaussian thermo whitepaper


What's Hss? If it's standard state, then Hss should be T-independent. Perhaps it should be "standard pressure"? Or did I read it wrong?

These terms just describe the enthalpy changes for the elements going from 0 K to 298 K. I think the Hss notation was left over from a previous version of Arkane, so I just replaced it.

alongd · 2019-05-23T14:02:30Z

arkane/encorr/mbac.py

+    long bonds. Use RMG for hydrogen because Open Babel can't do it for
+    mysterious reasons.
+    """
+    if len(nums) == 2 and all(n == 1 for n in nums):


Perhaps if nums == [1, 1] instead?
Why hard code for H2, or am I not reading it correctly?

I wanted to make sure that the check works for all sequence types and not just lists.

I'm hardcoding for H2 because I've observed in the past that Open Babel would not recognize correctly that two hydrogens are bonded to each other.

Will if list(nums) == [1, 1] work?

OK. Could you add a small comment explaining why H2 was hardcoded?

It's mentioned in the docstring of the function. Is that not enough?

Perfect! (I missed it)

alongd · 2019-05-23T14:06:22Z

arkane/encorr/mbac.py

+        mol.fromXYZ(nums, coords)
+    else:
+        xyz = '{}\n\n'.format(len(nums))
+        xyz += '\n'.join('{0}  {1[0]: .10f}  {1[1]: .10f}  {1[2]: .10f}'.format(s, c) for s, c in zip(coords, nums))


Can you take another look here? Looks like coords and nums should be switched in the zip

Oops, great catch.

alongd · 2019-05-23T14:10:08Z

documentation/source/users/arkane/input.rst

 ``useHinderedRotors``       ``True`` (by default) if hindered rotors are used, ``False`` if not
 ``useAtomCorrections``      ``True`` (by default) if atom corrections are used, ``False`` if not
 ``useBondCorrections``      ``True`` if bond corrections are used, ``False`` (by default) if not
+``bondCorrectionType``      ``'p'`` for Petersson-type or ``'m'`` for Melius-type bond additivity corrections


Could you specify what the default is?

alongd · 2019-05-23T14:11:09Z

documentation/source/users/arkane/input.rst

 atom corrections can be turned off by setting ``useAtomCorrections`` to ``False``.

-The ``bond`` parameter is used to apply bond corrections (BC) for a given ``modelChemistry``.
+The ``bond`` parameter is used to apply bond corrections (BC) for a given ``modelChemistry`` if using Petersson-type


Can we change to apply bond corrections (BC) into to apply bond additivity corrections (BAC), and then just use BAC later on?

alongd · 2019-05-23T14:12:18Z

documentation/source/users/arkane/input.rst

-The ``bond`` parameter is used to apply bond corrections (BC) for a given ``modelChemistry``.
+The ``bond`` parameter is used to apply bond corrections (BC) for a given ``modelChemistry`` if using Petersson-type
+bond additivity corrections (BACs) (``bondCorrectionType = 'p'``). When using Melius-type bond additivity corrections
+(``bondCorrectionType = 'm'``), specifying bonds is not required because the molecular connectivity is automatically


Perhaps change specifying bonds is not required into specifying bond orders is not required

It means that the bonds dictionary is not required in the input file. To highlight that, I added `` around it so that it gets parsed as plain text in the documentation.

mliu49 · 2019-05-23T14:53:13Z

arkane/encorr/corr.py

+        else:
+            raise Exception(
+                'Element "{}" is not yet supported in Arkane.'
+                ' To include it, add its experimental heat of formation'.format(symbol)


Could consider indicating where one would add it.

@cgrambow, was this comment addressed?

Now it has been.

mliu49 · 2019-05-23T15:02:59Z

arkane/encorr/pbac.py

+        if symbol in params:
+            bac += count * params[symbol] * 4184.0
+        elif symbol[::-1] in params:
+            bac += count * params[symbol[::-1]] * 4184.0


This method will break once we start supporting elements with two letter symbols, right? Although it might be a moot point if we're deprecating this type of correction anyway.

Didn't even think of that. I replaced it with a regex which works for all possible scenarios.

cgrambow · 2019-05-23T17:50:25Z

I think I addressed all your comments and made some fixup commits. I just realized that I caused a circular import issue, so I'll fix that, too.

cgrambow · 2019-05-23T18:02:15Z

Added two more fixup commits. I introduced an exceptions.py file in the main Arkane directory and added the exceptions there to get rid of the circular imports.

mliu49 · 2019-05-23T18:05:06Z

arkane/encorr/corr.py

    it is consistent with the normal gas-phase reference states.
    Optionally, correct the energy using bond additivity corrections.
+
+    model_chemistry: The model chemistry, typically specified as method/basis.


You can format the arguments to get nice styling when we compile the API documentation. Minimum would be adding the Args keyword and indenting:

Args: model_chemistry: The model chemistry, typically specified as method/basis. ...

See https://sphinxcontrib-napoleon.readthedocs.io/en/latest/example_google.html for even more details.

Good to know!

alongd · 2019-07-11T18:33:36Z

@cgrambow, can you rebase?
What else is needed for this PR?
Let's get it in before the Py3 transision

cgrambow · 2019-07-11T21:49:25Z

Squashed and rebased. Nothing else is needed for this PR, except for BAC data, but we don't have that yet, so we can add it after Py3 transition.

alongd · 2019-07-12T00:26:03Z

Could you take a look at the Codacy comment?

cgrambow · 2019-07-12T14:32:23Z

The Codacy comment can be ignored. Apparently Codacy isn't smart enough to figure out that c in coords can be indexed.

Make new encorr package which handles all energy corrections. Store all correction parameters in a dedicated data file. Handle bond additivity corrections in its own module, so that new types of bond additivity corrections can easily be implemented in the future.

Instead of correcting the electronic energy using explicitly defined bond types, Anantharaman and Melius (JPCA, 2005) only use three fitted parameters per atom type to correct the energy. Implement the framework for this type of correction in the encorr package.

bondCorrectionType can be set to 'p' to use conventional bond additivity corrections or to 'm' to use the newly implemented Melius-type corrections.

…rrections

alongd

Thanks for this contribution!

cgrambow added Arkane Status: Ready for Review PR is complete and ready to be reviewed Status: Discussion Needed Discussion needed for further progress labels May 22, 2019

cgrambow requested review from alongd and mliu49 May 22, 2019 22:54

cgrambow self-assigned this May 22, 2019

alongd reviewed May 23, 2019

View reviewed changes

mliu49 reviewed May 23, 2019

View reviewed changes

cgrambow force-pushed the new_bacs branch from 274931d to fb947cf Compare July 11, 2019 19:57

alongd mentioned this pull request Jul 11, 2019

Arkane output fixes #1607

Merged

cgrambow force-pushed the new_bacs branch 2 times, most recently from 5468be1 to 4f792b4 Compare July 11, 2019 21:44

cgrambow force-pushed the new_bacs branch from 4f792b4 to d396d60 Compare July 12, 2019 14:29

Colin Grambow added 4 commits July 12, 2019 14:09

Add bondCorrectionType option to Arkane

2aadb6b

bondCorrectionType can be set to 'p' to use conventional bond additivity corrections or to 'm' to use the newly implemented Melius-type corrections.

Update Arkane documentation to reflect Melius-type bond additivity co…

cdca610

…rrections

cgrambow force-pushed the new_bacs branch from d396d60 to cdca610 Compare July 12, 2019 18:17

alongd approved these changes Jul 12, 2019

View reviewed changes

mliu49 merged commit 58c1a06 into master Jul 12, 2019

mliu49 deleted the new_bacs branch July 12, 2019 20:01

mliu49 mentioned this pull request Jul 23, 2019

RMG-Py v2.4.1 release #1664

Merged

Conversation

cgrambow commented May 22, 2019

Motivation or Problem

Description of Changes

Testing

Discussion

Uh oh!

codecov bot commented May 23, 2019 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Codecov Report

Uh oh!

alongd left a comment

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

cgrambow commented May 23, 2019

Uh oh!

cgrambow commented May 23, 2019

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

alongd commented Jul 11, 2019

Uh oh!

cgrambow commented Jul 11, 2019

Uh oh!

alongd commented Jul 12, 2019

Uh oh!

cgrambow commented Jul 12, 2019

Uh oh!

alongd left a comment

Choose a reason for hiding this comment

codecov bot commented May 23, 2019 •

edited

Loading