Determine a species linearity from xyz #1601

alongd · 2019-05-14T02:50:04Z

Arkane uses the species linearity when projecting out rotor frequencies.
This PR automates the determination of the linearity attribute, so now it is optional in input files.
If the problem can be reduced in dimensionality, we first identify that. We then check that all atoms from the 3rd one on form a straight line (either in 2D space or in 3D space) within some precision with the first two atoms.
Examples and documentations were updated. Tests were added.

codecov · 2019-05-14T04:28:14Z

Codecov Report

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

@@            Coverage Diff             @@
##           master    #1601      +/-   ##
==========================================
+ Coverage   41.65%   41.67%   +0.01%     
==========================================
  Files         176      176              
  Lines       29102    29102              
  Branches     5988     5988              
==========================================
+ Hits        12122    12127       +5     
+ Misses      16142    16138       -4     
+ Partials      838      837       -1

Impacted Files	Coverage Δ
rmgpy/data/kinetics/family.py	`52.97% <0%> (+0.23%)`	⬆️

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 b944ddb...a66ffee. Read the comment docs.

cgrambow · 2019-05-14T17:25:06Z

arkane/statmech.py

+def nearly_equal(a, b, sig_fig=5):
+    """
+    A helper function to determine whether two floats are nearly equal.
+    Can be replaced by math.isclose() in Python3


You could also consider using numpy.isclose.

cgrambow · 2019-05-14T17:28:22Z

arkane/statmech.py

+
+    x, y, z = [], [], []
+
+    for coordinate in coordinates:


You can replace this using x, y, z = zip(*coordinates).

cgrambow · 2019-05-14T17:42:09Z

arkane/statmech.py

+        y.append(coordinate[1])
+        z.append(coordinate[2])
+
+    # first check whether the problem can be reduced into one or two dimensions


I think instead of doing all of this, it might be simpler to just check the angles between all atom triples and make sure they're all close to 180 degrees. To do so, you could calculate the distance vectors between all atoms and then check that the angle between any two of them is close to 180 degrees (or -180 degrees). A very efficient implementation of this should be possible with Numpy. Not sure what a suitable tolerance would be. I don't know if this is a better idea; let me know what you think.

For example, a tensor containing all distance vectors can be calculated with Numpy broadcasting in a single line:

d = -np.array([c[..., np.newaxis] - c[np.newaxis, ...] for c in coordinates.T])

which then gives the vector between atom 0 and atom 1 as d[:, 0, 1].

I think both methods are equivalent in essence, but your suggestion is significantly more friendly to implement / maintain / debug. I'd happily adopt it.

cgrambow · 2019-05-14T17:45:22Z

arkane/statmech.py

+    if linear is None:
+        linear = is_linear(coordinates, label)
+        not_txt = '' if linear else ' not'
+        logging.info('determined species {0}{1} to be linear'.format(label, not_txt))


Maybe only print this if you detect that the species is linear?

cgrambow · 2019-05-14T17:48:03Z

arkane/statmechTest.py

+            [-1.02600933, 2.12845307, 0.00000000],
+            [-0.77966935, 0.95278385, 0.00000000],
+            [-1.23666197, 3.17751246, 0.00000000],
+            [-0.56023545, -0.09447399, 0.00000000]])  # C2H2, just 0.5 degree off from linearity, so NOT linear


What do you think would be a good tolerance for determining linearity? Are there cases where we think that 0.5 degrees should be linear?

alongd · 2019-05-16T03:23:22Z

The tolerance is indeed a mystery.
I looked up CCCBDB for some linear molecules at a lousy (B3LYP/6-31G*) and a fancy method (CCSD(T)=FULL/aug-cc-pVQZ), but for CO2 and C2H2 both methods had exactly 180.0 (or 0.0) degrees between all triples. Thinking that they probably used Z-Mats, I ran CO2, C2H2, HCN, NCN(S), NCN(T), C#CC#C, and C#CC#CC#CC#C at wb97xd/6-311+G(d,p) using ARC (no ZMats), but got the same thing: all angles were dead on 180.0 (or 0.0) degrees (I used Gaussian).

So I plan on using 0.2 degrees for the tolerance, just based on my hunch, but I'm open for sugestions.

I added fixup commits and rebased, let me know when to squash and whether you have any other comments

Thanks for helping me improve the algorithm!

cgrambow · 2019-05-16T14:16:35Z

arkane/statmech.py

+        if i > 1:
+            u1 = d[:, 0, 1] / numpy.linalg.norm(d[:, 0, 1])  # unit vector between atoms 0 and 1
+            u2 = d[:, 1, i] / numpy.linalg.norm(d[:, 1, i])  # unit vector between atoms 1 and i
+            a = math.degrees(numpy.arccos(numpy.clip(numpy.dot(u1, u2), -1.0, 1.0)))  # angle between atoms 0, 1, i


Are you using numpy.clip to avoid invalid arguments for numpy.arccos if there are floating point errors in the calculation of the dot product?

Yes, I thought to use cautious although my examples ran fine without it

cgrambow · 2019-05-16T14:19:51Z

arkane/statmech.py

-        return False
+    # A tensor containing all distance vectors in the molecule
+    d = -numpy.array([c[:, numpy.newaxis] - c[numpy.newaxis, :] for c in coordinates.T])
+    for i, c in enumerate(coordinates):


You don't actually need to iterate over coordinates, correct? I'd just do
for i in range(2, len(coordinates))
and remove the if-statement on the next line.

absolutely!

cgrambow · 2019-05-16T14:32:25Z

That actually really interesting and good to know. 0.2 is probably ok for now, maybe even smaller based on your observations? We can adjust as we learn more. I requested one small change. After that you can squash.

alongd · 2019-05-17T13:17:24Z

Thanks for all the comments, the PR is in a much better shape now. Let me know if there are any other issues, otherwise this should be good to go

alongd · 2019-05-17T13:28:43Z

@cgrambow, I changed the tolerance to 0.1 degrees and added a descriptive comment (squashed)

cgrambow

Looks good. Just one more if-statement that you can remove.

cgrambow · 2019-05-17T14:33:21Z

arkane/statmech.py

+    # A tensor containing all distance vectors in the molecule
+    d = -numpy.array([c[:, numpy.newaxis] - c[numpy.newaxis, :] for c in coordinates.T])
+    for i in range(2, len(coordinates)):
+        if i > 1:


You don't need this if-statement anymore since i now starts at 2.

Thanks, removed

Left it in C2H6 for the example

alongd added the Arkane label May 14, 2019

alongd requested a review from cgrambow May 14, 2019 02:50

alongd self-assigned this May 14, 2019

cgrambow requested changes May 14, 2019

View reviewed changes

alongd force-pushed the linearity branch from 6ff92c3 to 5ff1e9e Compare May 16, 2019 03:00

cgrambow requested changes May 16, 2019

View reviewed changes

alongd force-pushed the linearity branch from 5ff1e9e to bfc2740 Compare May 17, 2019 13:16

alongd force-pushed the linearity branch from bfc2740 to e25198f Compare May 17, 2019 13:28

cgrambow requested changes May 17, 2019

View reviewed changes

alongd added 5 commits May 17, 2019 10:40

Determine species linearity from xyz if not given

c09abb4

Tests: is_linear() is Arkane statmech

1a096cf

Removed the 'linear' keyword from most of Arkane's examples

422ac40

Left it in C2H6 for the example

Updated documentation to reflect that 'linear' is optional

3095b87

Import numpy as np in Arkane statmech

a66ffee

alongd force-pushed the linearity branch from e25198f to a66ffee Compare May 17, 2019 14:41

cgrambow approved these changes May 17, 2019

View reviewed changes

cgrambow merged commit 19837cd into master May 17, 2019

cgrambow deleted the linearity branch May 17, 2019 15:29

mliu49 mentioned this pull request May 23, 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 a species linearity from xyz #1601

Determine a species linearity from xyz #1601

alongd commented May 14, 2019 •

edited

codecov bot commented May 14, 2019 •

edited

cgrambow May 14, 2019

cgrambow May 14, 2019

cgrambow May 14, 2019

alongd May 16, 2019

cgrambow May 14, 2019

cgrambow May 14, 2019

alongd commented May 16, 2019

cgrambow May 16, 2019

alongd May 17, 2019

cgrambow May 16, 2019

alongd May 17, 2019

cgrambow commented May 16, 2019

alongd commented May 17, 2019

alongd commented May 17, 2019

cgrambow left a comment

cgrambow May 17, 2019

alongd May 17, 2019

Determine a species linearity from xyz #1601

Determine a species linearity from xyz #1601

Conversation

alongd commented May 14, 2019 • edited

codecov bot commented May 14, 2019 • edited

Codecov Report

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 May 16, 2019

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

cgrambow commented May 16, 2019

alongd commented May 17, 2019

alongd commented May 17, 2019

cgrambow left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

alongd commented May 14, 2019 •

edited

codecov bot commented May 14, 2019 •

edited