Determine the rotor symmetry number from the potential scan #1526

alongd · 2018-12-19T16:44:43Z

Motivation or Problem

We'd like to allow Arkane to automatically determine rotor symmetry numbers. It is essential for automated thermo and rate codes using Arkane, like ARC. This feature is purposely added to Arkane and not ARC to benefit all users (if it works correctly, it could assist users to avoid making errors).

Description of Changes

A function for rotor symmetry determination was added to statmech.py.
The primary criterion used for determining whether a scan has symmetry or not is a comparison between peak/valley height difference and the worst peak/valley resolution.
This feature is able to differentiate between the following challenging cases, the one on the left is not symmetric (although it might look like it is at first sight), and the second has a symmetry number of 3:

Testing

Tests for the rotors plotted above for ethylamine were added. A statmechTest.py file was created.

Reviewer Tips

Try running Arkane with two arbitrary rotor scans, one with symmetry and one without, and see (in the log file) whether the determined symmetry is correct.

codecov · 2018-12-19T18:01:44Z

Codecov Report

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

@@            Coverage Diff             @@
##           master    #1526      +/-   ##
==========================================
- Coverage   42.08%   42.07%   -0.02%     
==========================================
  Files         165      165              
  Lines       27824    27824              
  Branches     5668     5668              
==========================================
- Hits        11711    11706       -5     
- Misses      15328    15332       +4     
- Partials      785      786       +1

Impacted Files	Coverage Δ
rmgpy/data/kinetics/family.py	`58.29% <0%> (-0.3%)`	⬇️

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 fb10a7b...383395c. Read the comment docs.

amarkpayne

Excellent PR, @alongd! There a couple of questions I have about how your code will function (see my review comments), and a couple suggestions for cleaning up the git history with this PR. I am also in the process of testing out the code, but I want to get these comments to you first.

amarkpayne · 2019-01-17T18:45:18Z

arkane/statmech.py

-                        rms_fourier = numpy.sqrt(numpy.sum((Vlist_fourier - Vlist) * (Vlist_fourier - Vlist))/ (len(Vlist) - 1)) / 4184.
-
+                    elif fit == 'best':
+                        rms_cosine = numpy.sqrt(numpy.sum((Vlist_cosine - v_list) * (Vlist_cosine - v_list)) /


You don't rename Vlist to v_list until the next commit though, right? If so, then if anyone were to reset their branch to this commit (without including your later commits) then this will throw an error. This isn't a big deal, but for consistent git history it might be worth moving the name change to this commit.

amarkpayne · 2019-01-17T20:45:39Z

arkane/statmech.py

+        # identify peaks and valleys, and determine worst resolutions in the scan
+        if i != 0 and i != len(energies) - 1:
+            # this is an intermediate point in the scan
+            if e > energies[i - 1] and e > energies[i + 1]:


A possible but optional thing you can do here would be to define variables for energies[i-1] and energies[i+1] since you reference these values a lot. The marginal code efficiency should be weighed against the readability of the code though, so this change is up to you.

amarkpayne · 2019-01-17T20:50:13Z

arkane/statmech.py

+                valleys.append(e)
+    # The number of peaks and valley must always be the same (what goes up must come down), if it isn't then there's
+    # something seriously wrong with the scan
+    if len(peaks) != len(valleys):


What happens though if there is a peak or valley right on one of the end points? You only check the intermediate points a few lines above this. Will your code catch this scenario?

Above we define peaks, valleys = list(), [energies[0]], so we take the first point into consideration, which is also represented by the last point (assuming a 360 deg scan). In a normal case both these points represent a valley, if it's a peak then this scan is problematic. This is checked elsewhere (only a warning is raised).

Gotcha, sorry that I missed this

amarkpayne · 2019-01-17T20:58:35Z

arkane/statmech.py

+    min_valley = min(valleys)
+    max_valley = max(valleys)
+    # Criterion 1: worst resolution
+    if max_peak - min_peak > worst_peak_resolution:


If my understanding of this is correct, you are assuming that for the symmetry number to be anything other than 1 that all of the peaks need to be the same size (within a tolerance). Is this true though? What if your scan had 4 peaks, with two larger peaks of height A and two smaller peaks of height B. If the peaks went ABAB, wouldn't this scan have a symmetry number of 2? Would your code catch this?

As we discussed offline, while this scenario would indeed have a symmetry number of two, we could not come up with a physical system that would have this property. Therefore, this technique should work in most if not all cases. No changes necessary.

amarkpayne · 2019-01-17T21:01:44Z

arkane/commonTest.py

@@ -270,24 +270,6 @@ def testTransitionStateStatmech(self):
        job.load()


-class TestStatmech(unittest.TestCase):


Can you move this change to the next commit where you add statmechTest.py? The rest of this commit can probably be combined with one of your earlier syntax commits as well.

amarkpayne · 2019-01-17T21:03:29Z

documentation/source/users/arkane/input.rst

-The ``symmetry`` parameter will usually equal either 1, 2 or 3. Below are examples of internal rotor scans with these
-commonly encountered symmetry numbers. First, ``symmetry = 3``:
+The ``symmetry`` parameter will usually equal either 1, 2 or 3. It could be determined automatically by Arkane
+(by simple not specifying it altogether), however it is always better to explicitly specify it if it is known. If it is


change by simple not specifying to by simply not specifying

amarkpayne · 2019-01-17T21:06:21Z

documentation/source/users/arkane/input_pdep.rst

@@ -351,8 +351,10 @@ scan in the following format::

 The ``Energy`` can be in units of ``kJ/mol``, ``J/mol``, ``cal/mol``, ``kcal/mol``, ``cm^-1`` or ``hartree``.

-The ``symmetry`` parameter will usually equal either 1, 2 or 3. Below are examples of internal rotor scans with these
-commonly encountered symmetry numbers. First, ``symmetry = 3``:
+The ``symmetry`` parameter will usually equal either 1, 2 or 3. It could be determined automatically by Arkane


Is there a way in .rst files to add a section of text from another file? That why, future updates to this in arkane/input.rst will automatically update the same section in input_pdep.rst. I am not sure if this is possible, but I will do some digging to find out.

Interesting, I don't know

It appears to be possible with the .. include :: statement, but it involves moving the text to another file. I am not sure if it is worth doing then, as it will clutter the documentation source files a little bit.

Let's leave it as is for now, I think we should change the pdep doc of Arkane to only include pdep stuff and avoid this repetition.

Also renamed Vlist as v_list

alongd · 2019-01-19T03:05:19Z

@amarkpayne, I addressed the comments, let me know what your thoughts are on the ABAB case.

amarkpayne

Everything looks good! I was also able to run the Arkane examples, and your technique got the symmetry right every time. Take a quick look at my final comments, then squash the fix-up commits and this is ready to be merged.

amarkpayne · 2019-01-20T22:09:32Z

arkane/statmech.py

+                valleys.append(e)
+    # The number of peaks and valley must always be the same (what goes up must come down), if it isn't then there's
+    # something seriously wrong with the scan
+    if len(peaks) != len(valleys):


Gotcha, sorry that I missed this

amarkpayne · 2019-01-20T22:11:35Z

arkane/statmech.py

+    min_valley = min(valleys)
+    max_valley = max(valleys)
+    # Criterion 1: worst resolution
+    if max_peak - min_peak > worst_peak_resolution:


As we discussed offline, while this scenario would indeed have a symmetry number of two, we could not come up with a physical system that would have this property. Therefore, this technique should work in most if not all cases. No changes necessary.

amarkpayne · 2019-01-20T22:21:27Z

arkane/statmechTest.py

+
+import unittest
+import os
+import shutil


My IDE suggests that this import statement is not needed, and I only see it in this file once. Can it be deleted?

amarkpayne · 2019-01-20T22:35:01Z

examples/arkane/species/C2H5/ethyl.py

-fit is fit of the scan data. It defaults to 'best', but can also be assigned as 'cosine' or 'fourier'"""
-rotors = [
-    HinderedRotor(scanLog=Log('ethyl_scan_72.log'), pivots=[1,2], top=[1,3,4], symmetry=6, fit='best')
-]


Is it worth doing the same thing to examples/arkane/species/Toluene_Hindered_Rotor?

I wanted to leave one example in with the explicit symmetry number (although it's well documented) for users who would like to specify it. But I'm open to other thoughts.

With test_rotor_symmetry_determination and a relocated test_gaussian_log_file_error Removed the TestStatmech class from commonTest.py

alongd · 2019-01-20T23:06:07Z

Done!
Thanks for the review and the helpful comments. Let me know if you have anything else in mind

alongd added the Arkane label Dec 19, 2018

alongd self-assigned this Dec 19, 2018

alongd requested a review from amarkpayne December 19, 2018 16:44

alongd force-pushed the rotsym branch from 0f15705 to 02f79cb Compare December 19, 2018 16:45

alongd force-pushed the rotsym branch 4 times, most recently from 55fad10 to 05f8196 Compare December 21, 2018 20:52

alongd force-pushed the rotsym branch 2 times, most recently from 530dbe6 to bea31c9 Compare January 3, 2019 15:28

amarkpayne requested changes Jan 17, 2019

View reviewed changes

alongd force-pushed the rotsym branch from bea31c9 to 7c488c5 Compare January 19, 2019 03:02

Minor: Whitespaces and linebreaks rearrangements in Arkane statmech.py

d68e314

Also renamed Vlist as v_list

alongd force-pushed the rotsym branch from 7c488c5 to 9491b19 Compare January 19, 2019 03:04

amarkpayne reviewed Jan 20, 2019

View reviewed changes

alongd added 5 commits January 20, 2019 18:01

Determine the rotor symmetry number from the potential scan

168a102

Minor: Corrected commented units in gaussian load scan

6beb888

Added statmechTest.py

a8abefb

With test_rotor_symmetry_determination and a relocated test_gaussian_log_file_error Removed the TestStatmech class from commonTest.py

Updated documentation with auto-rotor-symmetry determination

d81a679

Modified Arkane examples to reflect auto rotor symmetry determination

383395c

alongd force-pushed the rotsym branch from 4b7b777 to 383395c Compare January 20, 2019 23:04

amarkpayne approved these changes Jan 21, 2019

View reviewed changes

amarkpayne merged commit f14c7e1 into master Jan 21, 2019

amarkpayne deleted the rotsym branch January 21, 2019 00:35

mliu49 mentioned this pull request May 15, 2019

RMG v2.4.0 (?) Release Planning #1591

Closed

5 tasks

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Determine the rotor symmetry number from the potential scan #1526

Determine the rotor symmetry number from the potential scan #1526

alongd commented Dec 19, 2018 •

edited

Loading

codecov bot commented Dec 19, 2018 •

edited

Loading

amarkpayne left a comment

amarkpayne Jan 17, 2019

amarkpayne Jan 17, 2019

amarkpayne Jan 17, 2019

alongd Jan 18, 2019

amarkpayne Jan 20, 2019

amarkpayne Jan 17, 2019

amarkpayne Jan 20, 2019

amarkpayne Jan 17, 2019

amarkpayne Jan 17, 2019

amarkpayne Jan 17, 2019

alongd Jan 19, 2019

amarkpayne Jan 19, 2019

alongd Jan 19, 2019 •

edited

Loading

alongd commented Jan 19, 2019

amarkpayne left a comment •

edited

Loading

amarkpayne Jan 20, 2019

amarkpayne Jan 20, 2019

amarkpayne Jan 20, 2019

amarkpayne Jan 20, 2019

alongd Jan 20, 2019

alongd commented Jan 20, 2019

		@@ -270,24 +270,6 @@ def testTransitionStateStatmech(self):
		job.load()


		class TestStatmech(unittest.TestCase):

Determine the rotor symmetry number from the potential scan #1526

Determine the rotor symmetry number from the potential scan #1526

Conversation

alongd commented Dec 19, 2018 • edited Loading

Motivation or Problem

Description of Changes

Testing

Reviewer Tips

codecov bot commented Dec 19, 2018 • edited Loading

Codecov Report

amarkpayne 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

Choose a reason for hiding this comment

Choose a reason for hiding this comment

alongd Jan 19, 2019 • edited Loading

Choose a reason for hiding this comment

alongd commented Jan 19, 2019

amarkpayne left a comment • 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

Choose a reason for hiding this comment

alongd commented Jan 20, 2019

alongd commented Dec 19, 2018 •

edited

Loading

codecov bot commented Dec 19, 2018 •

edited

Loading

alongd Jan 19, 2019 •

edited

Loading

amarkpayne left a comment •

edited

Loading