Multi-Dimensional Rotors #1640

mjohnson541 · 2019-06-22T03:02:55Z

This PR enables ND empirical and semi-empirical rotor calculations along with 2D quantum mechanical rotor calculations.

ND empirical and semi-empirical rotors rely on scipy integration and interpolation algorithms. Frequencies are estimated by finding the hessian of a polynomial fit around the origin point of the rotor scans.

This PR enables 2D dimensional coupled rotor calculations using the 2D-NS approximation within Arkane using the external software Q2DTor written by David Ferro Costas and Antonio Fernandez Ramos. It adds a wrapper class in Arkane for computing the statistical mechanical properties of 2D coupled rotors using Q2DTor and adds rotor processing for 2D coupled rotors.

Installing my fork of Q2DTor:

git clone https://github.com/mjohnson541/Q2DTor.git
git checkout arkane

add to your .bashrc or .bash_profile (as appropriate):
export Q2DTor=/.../Q2DTor/src/Q2DTor.py (replacing ... with appropriate absolute path)

I've added an example for calculating the thermo of CH2CHOOH. However the full example takes ~5-8 hours of time to run. The wrapper (HinderedRotor2D) automatically checks for the presence of the .evals output file from Q2DTor that contains the eigenvalues of the Hamiltonian. If it finds this file it automatically skips the calculations and uses the eigenvalues to compute statistical mechanical properties. To shorten the example I've added all of the Q2DTor output files to the example allowing it to find the .evals file and finish rapidly.

alongd · 2019-07-11T18:39:19Z

@mjohnson541, a reminder to squash and rebase.

mjohnson541 · 2019-07-12T20:20:14Z

There's a reason I hadn't made a request for review on this one yet. There a tests issue that I need to resolve.

mjohnson541 · 2019-07-12T20:21:02Z

I'll try to have it ready soon.

mjohnson541 · 2019-07-18T18:06:18Z

Ok, this should be ready for review!

mjohnson541 · 2019-07-29T16:51:35Z

@alongd are you going to be able to review this?

alongd · 2019-07-29T18:03:02Z

Sure, it's just 1,400 files. I'm counting on the Hackathon...
I assume the comments from #1585 were addressed?
Meanwhile, do you mind rebasing?

mjohnson541 · 2019-07-29T19:03:16Z

I think I forgot those, I'll take a look tonight and get it setup by tomorrow.

mjohnson541 · 2019-07-29T19:03:42Z

Don't worry 99.9% of those files are Gaussian files

codecov · 2019-07-29T21:15:54Z

Codecov Report

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

@@            Coverage Diff             @@
##           master    #1640      +/-   ##
==========================================
- Coverage   41.87%   41.71%   -0.16%     
==========================================
  Files         176      176              
  Lines       29369    29359      -10     
  Branches     6059     6053       -6     
==========================================
- Hits        12297    12246      -51     
- Misses      16192    16244      +52     
+ Partials      880      869      -11

Impacted Files	Coverage Δ
rmgpy/tools/merge_models.py	`14.28% <0%> (-32.63%)`	⬇️
rmgpy/rmg/model.py	`38.46% <0%> (-3.25%)`	⬇️
rmgpy/data/kinetics/groups.py	`15.15% <0%> (-0.26%)`	⬇️
rmgpy/data/thermo.py	`61.86% <0%> (-0.07%)`	⬇️
rmgpy/reaction.py	`0% <0%> (ø)`	⬆️
rmgpy/quantity.py	`0% <0%> (ø)`	⬆️
rmgpy/molecule/molecule.py	`0% <0%> (ø)`	⬆️
rmgpy/data/kinetics/family.py	`52.9% <0%> (+0.23%)`	⬆️
rmgpy/rmg/react.py	`87.03% <0%> (+0.24%)`	⬆️
rmgpy/rmg/input.py	`35.77% <0%> (+0.38%)`	⬆️

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 69122e6...2ae9005. Read the comment docs.

codecov · 2019-07-29T21:15:54Z

Codecov Report

Merging #1640 into master will increase coverage by 0.22%.
The diff coverage is 59.51%.

@@            Coverage Diff             @@
##           master    #1640      +/-   ##
==========================================
+ Coverage   41.56%   41.78%   +0.22%     
==========================================
  Files         176      177       +1     
  Lines       29840    30207     +367     
  Branches     6192     6252      +60     
==========================================
+ Hits        12403    12622     +219     
- Misses      16529    16661     +132     
- Partials      908      924      +16

Impacted Files	Coverage Δ
rmgpy/quantity.py	`0% <0%> (ø)`	⬆️
rmgpy/statmech/ndTorsions.py	`59.67% <59.67%> (ø)`

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 114bd34...912f7b7. Read the comment docs.

mjohnson541 · 2019-07-29T23:53:40Z

This should be ready for review now.

alongd

Thanks for this PR!
Please see my comments below

rmgpy/statmech/conformer.pyx

arkane/statmech.py

alongd · 2019-07-31T19:40:29Z

arkane/statmech.py

@@ -1135,3 +1136,58 @@ def determine_rotor_symmetry(energies, label, pivots):
        logging.info('Determined a symmetry number of {0} for rotor of species {1} between pivots {2}'
                     ' based on the {3}.'.format(symmetry, label, pivots, reason))
    return symmetry
+
+class HinderedRotor2D(Mode):


Seems more natural to have this in rmgpy/statmech/torsion.pyx
We can think whether we want it to inherit from Torsion or not (and if not, drop a line explaining why)

So I invested a lot of time trying to get this to work in torsion.pyx, but there's a circular import that occurs because HinderedRotor2D and HinderedRotorClassicalND unlike the other rotor functions both need quantum file reading related imports so I decided to locate them in statmech.py where this didn't create a circular import.

OK, let's try to make this move as we discussed offline

alongd · 2019-07-31T19:40:51Z

arkane/statmech.py

@@ -1135,3 +1136,58 @@ def determine_rotor_symmetry(energies, label, pivots):
        logging.info('Determined a symmetry number of {0} for rotor of species {1} between pivots {2}'
                     ' based on the {3}.'.format(symmetry, label, pivots, reason))
    return symmetry
+
+class HinderedRotor2D(Mode):


Do we have tests for this class?

alongd · 2019-07-31T19:43:37Z

arkane/statmech.py

+
+   def getTorsions(self):
+       """
+       determine torsions, not entirely necessary for 2D-NS, but important for E2DT


please add (in parentheses?) the definitions of NS ans E2DT

arkane/statmech.py

alongd · 2019-08-01T01:11:37Z

arkane/statmech.py

@@ -1175,6 +1175,32 @@ def determine_rotor_symmetry(energies, label, pivots):
                     ' based on the {3}.'.format(symmetry, label, pivots, reason))
    return symmetry

+def is_linear(coordinates):


We already had the is_linear function in Arkane. I think it got deleted somewhere in this PR, and from briefly looking at it the added function seems extremely similar if not identical. Please change the PR so the function doesn't get deleted and re-added

We did, it got deleted in a prior PR and I had to recreate it after I rebased.

(It's deleted on master)

Can you double-check? I see it on master here:
https://github.com/ReactionMechanismGenerator/RMG-Py/blob/master/arkane/statmech.py#L737

I was incorrect, but I've squashed them so it shouldn't show up as deleted and created in a different commit

rmgpy/statmech/torsion.pyx

alongd · 2019-08-01T01:14:21Z

arkane/mainTest.py

@@ -55,6 +56,11 @@ def setUpClass(cls):
        cls.base_path = os.path.join(os.path.dirname(os.path.dirname(rmgpy.__file__)), 'examples', 'arkane')
        cls.failed = []
        cls.example_types = ['species', 'reactions', 'explorer', 'networks']
+        ch2ooh_path = os.path.join(cls.base_path,'species','CH2CHOOH','CH2CHOOHscans')


space between arguments? (below as well)

Can you take a second look?
It's minor, but I don't think its resolved.
It should be os.path.join(cls.base_path, 'species', 'CH2CHOOH', 'CH2CHOOHscans'),
but I still see os.path.join(cls.base_path,'species','CH2CHOOH','CH2CHOOHscans').

same thing below, and possibly other places in this PR

arkane/mainTest.py

alongd

I went over just a couple of changes, but didn't see the changes in the code. Were the changes pushed? Am I not looking at it correctly? Let's resolve it offline and I'll continue reviewing. Thanks!

alongd · 2019-08-05T14:47:18Z

arkane/mainTest.py

@@ -55,6 +56,11 @@ def setUpClass(cls):
        cls.base_path = os.path.join(os.path.dirname(os.path.dirname(rmgpy.__file__)), 'examples', 'arkane')
        cls.failed = []
        cls.example_types = ['species', 'reactions', 'explorer', 'networks']
+        ch2ooh_path = os.path.join(cls.base_path,'species','CH2CHOOH','CH2CHOOHscans')


Can you take a second look?
It's minor, but I don't think its resolved.
It should be os.path.join(cls.base_path, 'species', 'CH2CHOOH', 'CH2CHOOHscans'),
but I still see os.path.join(cls.base_path,'species','CH2CHOOH','CH2CHOOHscans').

same thing below, and possibly other places in this PR

alongd · 2019-08-05T14:50:21Z

arkane/statmech.py

+            elif pivots[1] in top:
+                pivot1 = pivots[1]
+                pivot2 = pivots[0]
+            else: raise Exception('Could not determine pivot atom.')


I still see else: raise Exception('Could not determine pivot atom.'), all in one line. Were the changes pushed?

mjohnson541 · 2019-08-05T14:53:52Z

Sorry I'll let you know once I've pushed. I haven't finished the documentation and then I need to pack everything into commits. I just wanted to check off things so I knew I fixed them.

mjohnson541 · 2019-08-05T16:32:42Z

Okay, I've pushed the changes!

mjohnson541 · 2019-08-05T22:37:03Z

It isn't packaged, so it can't be install with pip or conda. I can add it to the travis build though.

alongd

I added a few more comment

alongd · 2019-08-06T00:00:37Z

arkane/statmech.py

@@ -1135,3 +1136,58 @@ def determine_rotor_symmetry(energies, label, pivots):
        logging.info('Determined a symmetry number of {0} for rotor of species {1} between pivots {2}'
                     ' based on the {3}.'.format(symmetry, label, pivots, reason))
    return symmetry
+
+class HinderedRotor2D(Mode):


OK, let's try to make this move as we discussed offline

alongd · 2019-08-06T00:01:40Z

arkane/statmech.py

+       """
+       use Q2DTor to generate a .ics file
+       """
+       out = subprocess.check_call(['python2',os.environ['Q2DTor'],self.name,'--init'],


OK, let me know where to look once it's up (or is it already?)

alongd · 2019-08-06T00:03:23Z

arkane/statmech.py

@@ -1618,6 +1622,17 @@ def readScan(self):
                    Es.append(Es[j])
                    xyzs.append(xyzs[j])

+            q = len(phis)


Can you double check the fix for the duplicate loops?
I still see two instances of

for i in xrange(N): # add the negative values to improve fit near 0.0 for j in xrange(q):

at the same level

alongd · 2019-08-06T00:04:29Z

examples/arkane/species/Toulene_Hindered_Rotor_SemiClassicalND/output.py

@@ -0,0 +1,74 @@
+# Coordinates for Toluene in Input Orientation (angstroms):


let me know the status of this comment

examples/arkane/species/Toulene_Hindered_Rotor_SemiClassicalND/toluene_HinderedRotor.py

alongd · 2019-08-06T00:09:28Z

rmgpy/quantity.py

@@ -405,7 +405,7 @@ def __str__(self):
            uncertainty = '[{0}]'.format(','.join(uncertainty))

        result = '{0}'.format(value)
-        if any(self.uncertainty != 0.0):
+        if (self.uncertainty>0).any():


Can you teach me what this did and what it does now?
I'm used to the format if any([a, b, c]):, but here it looks like we're checking whether [a, b, c] > 0 instead of entry >0 for entry in [a, b, c]

Basically this gave me an error related to comparing an array with a value whenever it was actually called. I think it just never got called in normal operation until I was trying to log stuff. Essentially the line is suppose to display uncertainties if they aren't all zero, but it resulted in an error until I fixed it.

This seems odd to me. You're still comparing an array to a value, so it doesn't make sense why this would work any differently...

You're right, but I can tell you for certain that the first one crashed everytime and the second one displays properly

alongd · 2019-08-06T00:10:16Z

rmgpy/statmech/conformer.pyx

@@ -262,6 +262,8 @@ cdef class Conformer(RMGObject):
                N += 1
            elif type(mode) == SphericalTopRotor:
                N += 3
+            elif hasattr(mode,'dof'):


please add a space between mode and 'dof'

mjohnson541 · 2019-08-06T00:30:40Z

@alongd Here's the problem with moving I need to calculate Internal reduced moment of inertia at each point in HR multi D (this improves the accuracy of the calculation relative to Arkane's original HR software) However Arkane is designed to have this computed through the Conformer object (which makes sense), and the Conformer object imports torsion.pyx.

Essentially in order to make this possible I either need to write a file of independent code that duplicates a lot of Conformer functionality and then redesign how HR multi D gets conformer information and sends it to the independent code or Arkane needs redesigned.

… rotors from projectRotors

…rtition functions fit the partition function to a spline and generate first and second derivatives

for semiclassical correct using rotor frequencies

get the hessian by fitting a ND quadratic polynomial around the minimum

add omnia to conda channels list add numdifftools to dependencies

unzip CH2CHOOH zip file if necessary before running tests don't delete key directories for arkane test jobs

…object

alongd

Thanks for this awesome contribution!

alongd self-requested a review July 11, 2019 18:38

mjohnson541 force-pushed the classicalrotors branch 2 times, most recently from 8f20ac1 to ec28610 Compare July 18, 2019 15:58

mjohnson541 self-assigned this Jul 19, 2019

mjohnson541 force-pushed the classicalrotors branch from c35dae5 to 2ae9005 Compare July 29, 2019 21:15

mjohnson541 force-pushed the classicalrotors branch 2 times, most recently from dcf0db8 to d79ae2f Compare July 29, 2019 23:12

alongd reviewed Aug 1, 2019

View reviewed changes

alongd added Arkane Before Py3 Should be merged before Python 3 transition labels Aug 1, 2019

alongd reviewed Aug 5, 2019

View reviewed changes

mjohnson541 force-pushed the classicalrotors branch from d79ae2f to 4a07e3e Compare August 5, 2019 16:32

mjohnson541 force-pushed the classicalrotors branch from 4a07e3e to 13af14e Compare August 5, 2019 20:43

alongd reviewed Aug 6, 2019

View reviewed changes

mjohnson541 force-pushed the classicalrotors branch from fbabdb0 to de84045 Compare August 6, 2019 00:50

mjohnson541 added 24 commits August 7, 2019 17:13

add example data for testing

7cccff4

make it possible to projectRotors from classical ND rotors

e22e13b

make it possible to get the frequencies corresponding to the internal…

ae2a646

… rotors from projectRotors

add HinderedRotorClassicalND class

e667dec

read in scan information

fce1291

fit scan data and reduced moment of inertia information, calculate pa…

0350864

…rtition functions fit the partition function to a spline and generate first and second derivatives

calculate partition function by integrating the classical integral

a4c4647

for semiclassical correct using rotor frequencies

add functions for computing thermodynamic properties of the ND rotor

f38bd3c

add function for running the rotor calculations

e7305dc

add HinderedRotor2D example in examples

9d1c826

add example for HinderedRotorClassicalND

e8a58a9

add frequency estimation for ND rotors

d9b0e24

get the hessian by fitting a ND quadratic polynomial around the minimum

make 1d-hr quantum by default

ed6d1d8

update dependencies

389561b

add omnia to conda channels list add numdifftools to dependencies

tweak pdep tests now that quantum HR is default

66a04b7

update mainTest.py

7174733

unzip CH2CHOOH zip file if necessary before running tests don't delete key directories for arkane test jobs

allow degrees of freedom to be supplied by a dof attribute on a Mode …

8ccb586

…object

add tests for HinderedRotor2D

35f7498

add test for HinderedRotorClassicalND

f767d38

fix quantity.py uncertainty display for ArrayQuantity

bad745e

add Q2DTor and Q2DTor installation documentation

c6f129d

add documentation for multidimensional hindered rotors

1b15227

install Q2DTor when running travis

e036ae3

add make command for q2dtor

912f7b7

mjohnson541 force-pushed the classicalrotors branch from 5aba167 to 912f7b7 Compare August 7, 2019 21:13

alongd approved these changes Aug 8, 2019

View reviewed changes

alongd merged commit fabc165 into master Aug 8, 2019

alongd deleted the classicalrotors branch August 8, 2019 00:32

mliu49 mentioned this pull request Nov 26, 2019

RMG v3.0.0 Release Planning #1830

Closed

mliu49 mentioned this pull request Dec 16, 2019

RMG-Py v3.0.0 Release #1852

Merged

		@@ -0,0 +1,74 @@
		# Coordinates for Toluene in Input Orientation (angstroms):

Multi-Dimensional Rotors #1640

Multi-Dimensional Rotors #1640

Conversation

mjohnson541 commented Jun 22, 2019

alongd commented Jul 11, 2019

mjohnson541 commented Jul 12, 2019

mjohnson541 commented Jul 12, 2019

mjohnson541 commented Jul 18, 2019

mjohnson541 commented Jul 29, 2019

alongd commented Jul 29, 2019

mjohnson541 commented Jul 29, 2019

mjohnson541 commented Jul 29, 2019

codecov bot commented Jul 29, 2019

Codecov Report

codecov bot commented Jul 29, 2019 • edited Loading

Codecov Report

mjohnson541 commented Jul 29, 2019

alongd left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

alongd Aug 5, 2019 • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

alongd left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

mjohnson541 commented Aug 5, 2019

mjohnson541 commented Aug 5, 2019

mjohnson541 commented Aug 5, 2019

alongd left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

mjohnson541 Aug 6, 2019 • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

mjohnson541 commented Aug 6, 2019

alongd left a comment

Choose a reason for hiding this comment

codecov bot commented Jul 29, 2019 •

edited

Loading

alongd Aug 5, 2019 •

edited

Loading

mjohnson541 Aug 6, 2019 •

edited

Loading