ReactionMechanismGenerator / AutoTST

AutoTST: A framework to perform automated transition state theory calculations
Other
32 stars 16 forks source link

Updated symmetry calculations #50

Closed nateharms closed 4 years ago

nateharms commented 4 years ago

Updated the way that AutoTST determines symmetry numbers. We were using the outdated symmetry package to calculate this but it was resulting in errors in determination of symmetry numbers (see #49). This PR updates the autotst.species.Conformer.calculate_symmetry_number method to use the symmetry number calculated by RMG's Species objects. Tests were updated accordingly as well.

codecov[bot] commented 4 years ago

Codecov Report

Merging #50 into master will increase coverage by 0.04%. The diff coverage is 100%.

Impacted file tree graph

@@            Coverage Diff             @@
##           master      #50      +/-   ##
==========================================
+ Coverage    63.2%   63.25%   +0.04%     
==========================================
  Files          27       27              
  Lines        4626     4629       +3     
==========================================
+ Hits         2924     2928       +4     
+ Misses       1702     1701       -1
Impacted Files Coverage Δ
autotst/species.py 71.35% <100%> (+0.07%) :arrow_up:
autotst/reactionTest.py 99.5% <100%> (+0.01%) :arrow_up:
autotst/speciesTest.py 98.75% <100%> (ø) :arrow_up:
autotst/calculator/statmechTest.py 95% <0%> (-1.25%) :arrow_down:
autotst/reaction.py 80.08% <0%> (+0.21%) :arrow_up:

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 e417241...8238aa5. Read the comment docs.

nateharms commented 4 years ago

Symmetry numbers calculated for TSs in previous checks pass... But I believe it breaks down a bit at that point. This symmetry method works via the connectivity of the RMGMolecule rather than on the positions. So, if we want to get it to work we have to do something like the following:

In [1]: from rmgpy.molecule import Molecule as RMGMolecule 
   ...: from rmgpy.species import Species as RMGSpecies 
   ...: from autotst.reaction import Reaction, TS                                                                            
WARNING:root:No kyotocabinet available

In [2]: rxn = Reaction("[CH3]+[OH]_C+[O]") #Expected symmetry number of 3  
   ...: rxn.get_labeled_reaction() 
   ...: ts = rxn.ts["forward"][0]                                                                                            

In [3]: numbers = ts.ase_molecule.numbers 
   ...: positions = ts.ase_molecule.positions                                                                                

In [4]: mol = RMGMolecule() 
   ...: mol.from_xyz(numbers, positions)                                                                                     

In [5]: species = RMGSpecies(molecule=[mol]) 
   ...: species.get_symmetry_number()                                                                                        
Out[5]: 3.0

I make these updates shortly.

nateharms commented 4 years ago

@skrsna @davidfarinajr “mol” is different than “self.rmg_molecule”. And the reason why I’m doing it this way is because the “from_xyz” method will add a bond between reacting atoms because they are positioned close to each other. e.g. 1 will be bound to 2 and 2 will be bound to 3. This connectivity is needed for the symmetry estimation because the estimations are based on the graph not the 3D geometry. If we were to use the “self.rmg_molecule” we would have the 3D geometry, but we would lack that bond created with “from_xyz” which is needed for our TS symmetry estimation. Here are some options:

Thoughts?

On Thu, Nov 21, 2019 at 16:56 davidfarinajr notifications@github.com wrote:

@davidfarinajr commented on this pull request.

In autotst/species.py https://github.com/ReactionMechanismGenerator/AutoTST/pull/50#discussion_r349338178 :

  • return symmetry_number

  • numbers = self.ase_molecule.numbers

  • positions = self.ase_molecule.positions

  • mol = RMGMolecule()

If ‘self.rmg_molecule’ has coords, I think we should use that. If not, maybe use the ‘update_coords_from’ method to create rmg mol from ase mol?

— You are receiving this because you were assigned. Reply to this email directly, view it on GitHub https://github.com/ReactionMechanismGenerator/AutoTST/pull/50?email_source=notifications&email_token=AE24N275ZUNNTBKRZPMS4JTQU372XA5CNFSM4JP2BGR2YY3PNVWWK3TUL52HS4DFWFIHK3DMKJSXC5LFON2FEZLWNFSXPKTDN5WW2ZLOORPWSZGOCMSVT3A#discussion_r349338178, or unsubscribe https://github.com/notifications/unsubscribe-auth/AE24N27OTXRQQRYPLSOGOIDQU372XANCNFSM4JP2BGRQ .

-- Please excuse any errors, the previous message was sent from a mobile device.

skrsna commented 4 years ago

sounds good to me to leave it as it is if making a RMGMolecule object doesn't take up lots of run time memory. I guess I got confused with RMGSpecies object cause @rwest commented that it's taking up lot of memory during resonance structure generation for catalysis. If @davidfarinajr agrees someone merge this.