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 :)
  • Loading branch information
enochd committed Jul 17, 2014
1 parent 48d24c8 commit 4599552
Show file tree
Hide file tree
Showing 3 changed files with 61 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
48 changes: 48 additions & 0 deletions rmgpy/statmech/conformer.pyx
Original file line number Diff line number Diff line change
Expand Up @@ -237,6 +237,54 @@ cdef class Conformer:

return center

cpdef getNumberDegreesOfFreedom(self):
"""
Return the number of degrees of freedom in a species object, which should be 3N
"""
cdef int N, Natoms, degreesOfFreedom, numberDegreesOfFreedom

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)):
try:
self.conformer.modes[i].semiclassical # found a HinderedRotor
N += 1
except:
pass
try:
self.conformer.modes[i].frequencies
N += len(self.conformer.modes[i].frequencies.value) # found the harmonic frequencies
except:
pass
try:
self.conformer.modes[i].mass # found the translational degrees of freedom
N += 3
except:
pass
if type(self.conformer.modes[i]) == NonlinearRotor:
N += 3
else:
pass
if type(self.conformer.modes[i]) == LinearRotor:
N += 2
else:
pass
if type(self.conformer.modes[i]) == KRotor:
N += 1
else:
pass
if type(self.conformer.modes[i]) == SphericalTopRotor:
N += 3
else:
pass
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 4599552

Please sign in to comment.