ReactionMechanismGenerator / AutoTST

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

Symmetry number estimations not always accurate #49

Open nateharms opened 4 years ago

nateharms commented 4 years ago

Symmetry numbers calculated using RMG's symmetry package tend to result in incorrect estimations of symmetry numbers. Below is an example of ethane which should have a symmetry number of 6. This issue arose when creating PR #48. Some investigation into alternate packages may remedy this.

In [1]: from autotst.species import Conformer                                                                                

In [2]: conf = Conformer("CC")                                                                                               

In [3]: conf.calculate_symmetry_number()                                                                                     
symmetry.py:253 write_input_file INFO Symmetry input file written to ./CC.symm
symmetry.py:219 parse INFO Point group: C1
Out[3]: 1
whgreen commented 4 years ago

Nate, I suggest you discuss with Mark Goldman. Bill

Sent from my phone.

On Nov 19, 2019, at 2:48 PM, Nate Harms notifications@github.com wrote:



Symmetry numbers calculated using RMG's symmetry package tend to result in incorrect estimations of symmetry numbers. Below is an example of ethane which should have a symmetry number of 6https://webbook.nist.gov/cgi/cbook.cgi?ID=C74840&Mask=800. This issue arose when creating PR #48https://github.com/ReactionMechanismGenerator/AutoTST/pull/48. Some investigation into alternate packages may remedy this.

In [1]: from autotst.species import Conformer

In [2]: conf = Conformer("CC")

In [3]: conf.calculate_symmetry_number() symmetry.py:253 write_input_file INFO Symmetry input file written to ./CC.symm symmetry.py:219 parse INFO Point group: C1 Out[3]: 1

— You are receiving this because you are subscribed to this thread. Reply to this email directly, view it on GitHubhttps://github.com/ReactionMechanismGenerator/AutoTST/issues/49?email_source=notifications&email_token=AAJL7MCWK5PEUZFUXXGQHQLQUQ7IFA5CNFSM4JPIO4UKYY3PNVWWK3TUL52HS4DFUVEXG43VMWVGG33NNVSW45C7NFSM4H2OKNAA, or unsubscribehttps://github.com/notifications/unsubscribe-auth/AAJL7MBCVXNL3JYTUMNETZ3QUQ7IFANCNFSM4JPIO4UA.

nateharms commented 4 years ago

Thanks for the recommendation, I just reached out to him.

goldmanm commented 4 years ago

RMG's symmetry algorithm includes both internal (from sigma bond rotation) and external symmetry (from rotation the whole molecule). Since the C-C bond in ethane can rotate a methyl group, it has an internal rotational symmetry of 3. Combined with the external symmetry you cite, the symmetry number RMG's symmetry algorithm should return is 18.

Running the current RMG package leads to:

>>> from rmgpy.species import Species
>>> s = Species().from_smiles('CC')
>>> s.get_symmetry_number()
18.0

This seems to work in RMG. I am not sure why it isn't working in AutoTST.

nateharms commented 4 years ago

Thanks, I'll switch it to this. We are currently doing this at the moment as a method of autotst.species.Conformer:

    def calculate_symmetry_number(self):
        from rmgpy.qm.symmetry import PointGroupCalculator
        from rmgpy.qm.qmdata import QMData

        atom_numbers = self.ase_molecule.get_atomic_numbers()
        coordinates = self.ase_molecule.get_positions()

        qmdata = QMData(
            groundStateDegeneracy=1,  # Only needed to check if valid QMData
            numberOfAtoms=len(atom_numbers),
            atomicNumbers=atom_numbers,
            atomCoords=(coordinates, str('angstrom')),
            energy=(0.0, str('kcal/mol'))  # Only needed to avoid error
        )
        settings = type(str(''), (), dict(symmetryPath=str(
            'symmetry'), scratchDirectory="."))()  # Creates anonymous class
        pgc = PointGroupCalculator(settings, self.smiles, qmdata)
        pg = pgc.calculate()
        #os.remove("{}.symm".format(self.smiles))

        if pg is not None:
            symmetry_number = pg.symmetry_number
        else:
            symmetry_number = 1

        return symmetry_number

This aspect comes from what was done in the first version of AutoTST, but it's looking like it's time for a change.