Skip to content

Commit

Permalink
Introduced getNumberDegreesOfFreedom into conformer class
Browse files Browse the repository at this point in the history
This partially addressed issue ReactionMechanismGenerator#213.

This function gives the total molecular degrees of freedom for a
species. I think it belongs in the conformer class, but please let me
know if you think otherwise. I've also included a couple unittests for
ethane and ethylene. I've recompiled RMG-Py and ensured that it's
working. That it's working is one thing, but apologies in advance for
the crappy programming style :)

fixup! Introduced `getNumberDegreesOfFreedom` into conformer class

Rather than look up self.conformer.modes[i] loads of times,  we look it up once
and store it in a local variable 'mode'. It is lookups that take time in Python,
so this would be faster (as well as hopefully easier to read and debug)

fixup! Introduced `getNumberDegreesOfFreedom` into conformer class

Better to explicitly check for what you are looking for, than assume that
the exception you catch in a try/catch block is the right one. Makes code
easier to read too.

My question is should we be checking for `Torsion` or is `HinderedRotor` the correct class?

fixup! Introduced `getNumberDegreesOfFreedom` into conformer class

Try/except could catch other exceptions, and suggests that the try ought to work.
Maybe this should actually be a type check on mode (is it an instance of `Vibration`, for example)?

fixup! Introduced `getNumberDegreesOfFreedom` into conformer class

The "else: pass" blocks are redundant.

You may consider chaining all the "if" blocks as "elif", and end with a final
"else:" that raises an "I don't know what this type of mode is!" type of error.

reformat code style, missing colon in `conformer.pyx`

This commit is in response to rwest's suggestions for your pull request
  • Loading branch information
enochd authored and connie committed Jul 17, 2014
1 parent d754af8 commit c96772e
Show file tree
Hide file tree
Showing 3 changed files with 51 additions and 1 deletion.
2 changes: 2 additions & 0 deletions rmgpy/statmech/conformer.pxd
Original file line number Diff line number Diff line change
Expand Up @@ -68,6 +68,8 @@ cdef class Conformer:
cpdef getSymmetricTopRotors(self)

cpdef list getActiveModes(self, bint activeJRotor=?, bint activeKRotor=?)

cpdef getNumberDegreesOfFreedom(self)

################################################################################

Expand Down
38 changes: 38 additions & 0 deletions rmgpy/statmech/conformer.pyx
Original file line number Diff line number Diff line change
Expand Up @@ -237,6 +237,44 @@ cdef class Conformer:

return center

cpdef getNumberDegreesOfFreedom(self):
"""
Return the number of degrees of freedom in a species object, which should be 3N,
and raises an exception if it is not.
"""
cdef int N, Natoms, degreesOfFreedom, numberDegreesOfFreedom
cdef Mode mode
N = 0
Natoms = len(self.conformer.mass.value)
# what the total number of degrees of freedom for the species should be
degreesOfFreedom = Natoms*3
for i in range(len(self.conformer.modes)):
mode = self.conformer.modes[i]
if isinstance(mode, HinderedRotor):
N += 1
if hasattr(mode,'frequencies'):
N += len(mode.frequencies.value) # found the harmonic frequencies
try:
mode.mass # found the translational degrees of freedom
N += 3
except:
pass
if type(mode) == NonlinearRotor:
N += 3
elif type(mode) == LinearRotor:
N += 2
elif type(mode) == KRotor:
N += 1
elif type(mode) == SphericalTopRotor:
N += 3
else:
raise TypeError("Mode type not supported")
numberDegreesOfFreedom = N
if N != degreesOfFreedom:
raise ValueError('The total degrees of molecular freedom for this species should be ' + str(degreesOfFreedom))

return numberDegreesOfFreedom

@cython.boundscheck(False)
@cython.wraparound(False)
cpdef numpy.ndarray getMomentOfInertiaTensor(self):
Expand Down
12 changes: 11 additions & 1 deletion rmgpy/statmech/conformerTest.py
Original file line number Diff line number Diff line change
Expand Up @@ -39,7 +39,7 @@
import scipy.interpolate

from rmgpy.statmech import Conformer, IdealGasTranslation, NonlinearRotor, HarmonicOscillator, \
LinearRotor, HinderedRotor
LinearRotor, HinderedRotor
import rmgpy.constants as constants

################################################################################
Expand Down Expand Up @@ -295,3 +295,13 @@ def test_getInternalReducedMomentOfInertia(self):
"""
I = self.conformer.getInternalReducedMomentOfInertia(pivots=[1,5], top1=[1,2,3,4])
self.assertAlmostEqual(I*constants.Na*1e23, 1.56768, 4)
def test_getNumberDegreesOfFreedom(self):
"""
Test the Conformer.getNumberDegreesOfFreedom() method.
"""
#this is for ethane:
numberDegreesOfFreedom = self.conformer.getNumberDegreesOfFreedom
self.assertTrue(numberDegreesOfFreedom, 24)
#this is for ethylene:
numberDegreesOfFreedom = self.ethylene.getNumberDegreesOfFreedom
self.assertTrue(numberDegreesOfFreedom, 18)

0 comments on commit c96772e

Please sign in to comment.