PageRenderTime 48ms CodeModel.GetById 19ms RepoModel.GetById 1ms app.codeStats 0ms

/rmgpy/statmech/conformerTest.py

https://github.com/jwallen/RMG-Py
Python | 281 lines | 166 code | 25 blank | 90 comment | 10 complexity | c9e0d2886baefa3f69ede3c36529d79c MD5 | raw file
    

  1. #!/usr/bin/env python
    2. 

  2. # encoding: utf-8
    3. 

  3. 

  4. ################################################################################
    5. 

  5. #
    6. 

  6. #   RMG - Reaction Mechanism Generator
    7. 

  7. #
    8. 

  8. #   Copyright (c) 2002-2009 Prof. William H. Green (whgreen@mit.edu) and the
    9. 

  9. #   RMG Team (rmg_dev@mit.edu)
    10. 

  10. #
    11. 

  11. #   Permission is hereby granted, free of charge, to any person obtaining a
    12. 

  12. #   copy of this software and associated documentation files (the "Software"),
    13. 

  13. #   to deal in the Software without restriction, including without limitation
    14. 

  14. #   the rights to use, copy, modify, merge, publish, distribute, sublicense,
    15. 

  15. #   and/or sell copies of the Software, and to permit persons to whom the
    16. 

  16. #   Software is furnished to do so, subject to the following conditions:
    17. 

  17. #
    18. 

  18. #   The above copyright notice and this permission notice shall be included in
    19. 

  19. #   all copies or substantial portions of the Software.
    20. 

  20. #
    21. 

  21. #   THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
    22. 

  22. #   IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
    23. 

  23. #   FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
    24. 

  24. #   THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
    25. 

  25. #   LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
    26. 

  26. #   FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
    27. 

  27. #   DEALINGS IN THE SOFTWARE.
    28. 

  28. #
    29. 

  29. ################################################################################
    30. 

  30. 

  31. """
    32. 

  32. This script contains unit tests of the :mod:`rmgpy.statmech.conformer` module.
    33. 

  33. """
    34. 

    35. 

  35. import unittest
    36. 

  36. import math
    37. 

  37. import numpy
    38. 

  38. import scipy.interpolate
    39. 

  39. 

  40. from rmgpy.statmech import Conformer, IdealGasTranslation, NonlinearRotor, HarmonicOscillator, \
    41. 

  41.                            LinearRotor, HinderedRotor
    42. 

  42. import rmgpy.constants as constants
    43. 

  43. 

  44. ################################################################################
    45. 

  45. 

  46. class TestConformer(unittest.TestCase):
    47. 

  47.     """
    48. 

  48.     Contains unit tests of the :class:`Conformer` class.
    49. 

  49.     """
    50. 

  50.     
    51. 

  51.     def setUp(self):
    52. 

  52.         """
    53. 

  53.         A function run before each unit test in this class.
    54. 

  54.         """
    55. 

  55.         self.ethylene = Conformer(
    56. 

  56.             E0 = (0.0,"kJ/mol"),
    57. 

  57.             modes = [
    58. 

  58.                 IdealGasTranslation(mass=(28.03,"amu")),
    59. 

  59.                 NonlinearRotor(inertia=([3.41526,16.6498,20.065],"amu*angstrom^2"), symmetry=4),
    60. 

  60.                 HarmonicOscillator(frequencies=([828.397,970.652,977.223,1052.93,1233.55,1367.56,1465.09,1672.25,3098.46,3111.7,3165.79,3193.54],"cm^-1")),
    61. 

  61.             ],
    62. 

  62.             spinMultiplicity = 1,
    63. 

  63.             opticalIsomers = 1,
    64. 

  64.         )
    65. 

  65.         self.oxygen = Conformer(
    66. 

  66.             E0 = (0.0,"kJ/mol"),
    67. 

  67.             modes = [
    68. 

  68.                 IdealGasTranslation(mass=(31.99,"amu")),
    69. 

  69.                 LinearRotor(inertia=(11.6056,"amu*angstrom^2"), symmetry=2),
    70. 

  70.                 HarmonicOscillator(frequencies=([1621.54],"cm^-1")),
    71. 

  71.             ],
    72. 

  72.             spinMultiplicity = 3,
    73. 

  73.             opticalIsomers = 1,
    74. 

  74.         )
    75. 

  75.         
    76. 

  76.         # The following data is for ethane at the CBS-QB3 level
    77. 

  77.         self.coordinates = numpy.array([
    78. 

  78.             [  0.0000,  0.0000,  0.0000],
    79. 

  79.             [ -0.0000, -0.0000,  1.0936],
    80. 

  80.             [  1.0430, -0.0000, -0.3288],
    81. 

  81.             [ -0.4484,  0.9417, -0.3288],
    82. 

  82.             [ -0.7609, -1.2051, -0.5580],
    83. 

  83.             [ -0.7609, -1.2051, -1.6516],
    84. 

  84.             [ -0.3125, -2.1468, -0.2292],
    85. 

  85.             [ -1.8039, -1.2051, -0.2293],
    86. 

  86.         ])
    87. 

  87.         self.number = numpy.array([6, 1, 1, 1, 6, 1, 1, 1])
    88. 

  88.         self.mass = numpy.array([12, 1.007825, 1.007825, 1.007825, 12, 1.007825, 1.007825, 1.007825])
    89. 

  89.         self.E0 = -93.5097
    90. 

  90.         self.conformer = Conformer(
    91. 

  91.             E0 = (self.E0,"kJ/mol"),
    92. 

  92.             modes = [
    93. 

  93.                 IdealGasTranslation(mass=(30.0469,"amu")),
    94. 

  94.                 NonlinearRotor(inertia=([6.27071,25.3832,25.3833],"amu*angstrom^2"), symmetry=6),
    95. 

  95.                 HarmonicOscillator(frequencies=([818.917,819.479,987.099,1206.76,1207.05,1396,1411.35,1489.73,1489.95,1492.49,1492.66,2995.36,2996.06,3040.77,3041,3065.86,3066.02],"cm^-1")),
    96. 

  96.                 HinderedRotor(inertia=(1.56768,"amu*angstrom^2"), symmetry=3, barrier=(2.69401,"kcal/mol"), quantum=False, semiclassical=False),
    97. 

  97.             ],
    98. 

  98.             spinMultiplicity = 1,
    99. 

  99.             opticalIsomers = 1,
    100. 

  100.             coordinates = (self.coordinates,"angstrom"),
    101. 

  101.             number = self.number,
    102. 

  102.             mass = (self.mass,"amu"),
    103. 

  103.         )
    104. 

  104.         
    105. 

  105.     def test_getPartitionFunction_ethylene(self):
    106. 

  106.         """
    107. 

  107.         Test the StatMech.getPartitionFunction() method for ethylene.
    108. 

  108.         """
    109. 

  109.         Tlist = numpy.array([300,500,1000,1500,2000])
    110. 

  110.         Qexplist = numpy.array([4.05311e+09, 4.19728e+10, 2.82309e+12, 7.51135e+13, 1.16538e+15])
    111. 

  111.         for T, Qexp in zip(Tlist, Qexplist):
    112. 

  112.             Qact = self.ethylene.getPartitionFunction(T)
    113. 

  113.             self.assertAlmostEqual(Qexp, Qact, delta=1e-4*Qexp)
    114. 

  114. 

  115.     def test_getHeatCapacity_ethylene(self):
    116. 

  116.         """
    117. 

  117.         Test the StatMech.getHeatCapacity() method for ethylene.
    118. 

  118.         """
    119. 

  119.         Tlist = numpy.array([300,500,1000,1500,2000])
    120. 

  120.         Cvexplist = numpy.array([5.11186, 7.40447, 11.1659, 13.1221, 14.1617]) * constants.R
    121. 

  121.         for T, Cvexp in zip(Tlist, Cvexplist):
    122. 

  122.             Cvact = self.ethylene.getHeatCapacity(T)
    123. 

  123.             self.assertAlmostEqual(Cvexp, Cvact, 3)
    124. 

  124.     
    125. 

  125.     def test_getEnthalpy_ethylene(self):
    126. 

  126.         """
    127. 

  127.         Test the StatMech.getEnthalpy() method for ethylene.
    128. 

  128.         """
    129. 

  129.         Tlist = numpy.array([300,500,1000,1500,2000])
    130. 

  130.         Hexplist = numpy.array([4.23129, 5.04826, 7.27337, 8.93167, 10.1223]) * constants.R * Tlist
    131. 

  131.         for T, Hexp in zip(Tlist, Hexplist):
    132. 

  132.             Hact = self.ethylene.getEnthalpy(T)
    133. 

  133.             self.assertAlmostEqual(Hexp, Hact, delta=1e-4*Hexp)
    134. 

  134.     
    135. 

  135.     def test_getEntropy_ethylene(self):
    136. 

  136.         """
    137. 

  137.         Test the StatMech.getEntropy() method for ethylene.
    138. 

  138.         """
    139. 

  139.         Tlist = numpy.array([300,500,1000,1500,2000])
    140. 

  140.         Sexplist = numpy.array([26.3540, 29.5085, 35.9422, 40.8817, 44.8142]) * constants.R
    141. 

  141.         for T, Sexp in zip(Tlist, Sexplist):
    142. 

  142.             Sact = self.ethylene.getEntropy(T)
    143. 

  143.             self.assertAlmostEqual(Sexp, Sact, 3)
    144. 

  144.     
    145. 

  145.     def test_getSumOfStates_ethylene(self):
    146. 

  146.         """
    147. 

  147.         Test the StatMech.getSumOfStates() method for ethylene.
    148. 

  148.         """
    149. 

  149.         Elist = numpy.arange(0, 5000*11.96, 2*11.96)
    150. 

  150.         sumStates = self.ethylene.getSumOfStates(Elist)
    151. 

  151.         densStates = self.ethylene.getDensityOfStates(Elist)
    152. 

  152.         for n in range(10, len(Elist)):
    153. 

  153.             self.assertTrue(0.8 < numpy.sum(densStates[0:n+1]) / sumStates[n] < 1.25, '{0} != {1}'.format(numpy.sum(densStates[0:n+1]), sumStates[n]))
    154. 

  154.             
    155. 

  155.     def test_getDensityOfStates_ethylene(self):
    156. 

  156.         """
    157. 

  157.         Test the StatMech.getDensityOfStates() method for ethylene.
    158. 

  158.         """
    159. 

  159.         Elist = numpy.arange(0, 5000*11.96, 2*11.96)
    160. 

  160.         densStates = self.ethylene.getDensityOfStates(Elist)
    161. 

  161.         T = 100
    162. 

  162.         Qact = numpy.sum(densStates * numpy.exp(-Elist / constants.R / T))
    163. 

  163.         Qexp = self.ethylene.getPartitionFunction(T)
    164. 

  164.         self.assertAlmostEqual(Qexp, Qact, delta=1e-1*Qexp)
    165. 

  165. 

  166.     def test_getPartitionFunction_oxygen(self):
    167. 

  167.         """
    168. 

  168.         Test the StatMech.getPartitionFunction() method for oxygen.
    169. 

  169.         """
    170. 

  170.         Tlist = numpy.array([300,500,1000,1500,2000])
    171. 

  171.         Qexplist = numpy.array([1.55584e+09, 9.38339e+09, 1.16459e+11, 5.51016e+11, 1.72794e+12])
    172. 

  172.         for T, Qexp in zip(Tlist, Qexplist):
    173. 

  173.             Qact = self.oxygen.getPartitionFunction(T)
    174. 

  174.             self.assertAlmostEqual(Qexp, Qact, delta=1e-4*Qexp)
    175. 

  175. 

  176.     def test_getHeatCapacity_oxygen(self):
    177. 

  177.         """
    178. 

  178.         Test the StatMech.getHeatCapacity() method for oxygen.
    179. 

  179.         """
    180. 

  180.         Tlist = numpy.array([300,500,1000,1500,2000])
    181. 

  181.         Cvexplist = numpy.array([3.52538, 3.70877, 4.14751, 4.32063, 4.39392]) * constants.R
    182. 

  182.         for T, Cvexp in zip(Tlist, Cvexplist):
    183. 

  183.             Cvact = self.oxygen.getHeatCapacity(T)
    184. 

  184.             self.assertAlmostEqual(Cvexp, Cvact, 3)
    185. 

  185.     
    186. 

  186.     def test_getEnthalpy_oxygen(self):
    187. 

  187.         """
    188. 

  188.         Test the StatMech.getEnthalpy() method for oxygen.
    189. 

  189.         """
    190. 

  190.         Tlist = numpy.array([300,500,1000,1500,2000])
    191. 

  191.         Hexplist = numpy.array([3.50326, 3.54432, 3.75062, 3.91623, 4.02765]) * constants.R * Tlist
    192. 

  192.         for T, Hexp in zip(Tlist, Hexplist):
    193. 

  193.             Hact = self.oxygen.getEnthalpy(T)
    194. 

  194.             self.assertAlmostEqual(Hexp, Hact, delta=1e-4*Hexp)
    195. 

  195.     
    196. 

  196.     def test_getEntropy_oxygen(self):
    197. 

  197.         """
    198. 

  198.         Test the StatMech.getEntropy() method for oxygen.
    199. 

  199.         """
    200. 

  200.         Tlist = numpy.array([300,500,1000,1500,2000])
    201. 

  201.         Sexplist = numpy.array([24.6685, 26.5065, 29.2314, 30.9513, 32.2056]) * constants.R
    202. 

  202.         for T, Sexp in zip(Tlist, Sexplist):
    203. 

  203.             Sact = self.oxygen.getEntropy(T)
    204. 

  204.             self.assertAlmostEqual(Sexp, Sact, 3)
    205. 

  205.     
    206. 

  206.     def test_getSumOfStates_oxygen(self):
    207. 

  207.         """
    208. 

  208.         Test the StatMech.getSumOfStates() method for oxygen.
    209. 

  209.         """
    210. 

  210.         Elist = numpy.arange(0, 5000*11.96, 2*11.96)
    211. 

  211.         sumStates = self.oxygen.getSumOfStates(Elist)
    212. 

  212.         densStates = self.oxygen.getDensityOfStates(Elist)
    213. 

  213.         for n in range(10, len(Elist)):
    214. 

  214.             self.assertTrue(0.8 < numpy.sum(densStates[0:n+1]) / sumStates[n] < 1.25, '{0} != {1}'.format(numpy.sum(densStates[0:n+1]), sumStates[n]))
    215. 

  215.             
    216. 

  216.     def test_getDensityOfStates_oxygen(self):
    217. 

  217.         """
    218. 

  218.         Test the StatMech.getDensityOfStates() method for oxygen.
    219. 

  219.         """
    220. 

  220.         Elist = numpy.arange(0, 5000*11.96, 2*11.96)
    221. 

  221.         densStates = self.oxygen.getDensityOfStates(Elist)
    222. 

  222.         T = 100
    223. 

  223.         Qact = numpy.sum(densStates * numpy.exp(-Elist / constants.R / T))
    224. 

  224.         Qexp = self.oxygen.getPartitionFunction(T)
    225. 

  225.         self.assertAlmostEqual(Qexp, Qact, delta=1e-1*Qexp)
    226. 

  226. 

  227.     def test_getTotalMass(self):
    228. 

  228.         """
    229. 

  229.         Test the Conformer.getTotalMass() method.
    230. 

  230.         """
    231. 

  231.         self.assertAlmostEqual(self.conformer.getTotalMass()*constants.Na*1000., numpy.sum(self.mass), 6)
    232. 

  232. 

  233.     def test_getCenterOfMass(self):
    234. 

  234.         """
    235. 

  235.         Test the Conformer.getCenterOfMass() method.
    236. 

  236.         """
    237. 

  237.         cm = self.conformer.getCenterOfMass()
    238. 

  238.         self.assertAlmostEqual(cm[0]*1e10, -0.38045, 4)
    239. 

  239.         self.assertAlmostEqual(cm[1]*1e10, -0.60255, 4)
    240. 

  240.         self.assertAlmostEqual(cm[2]*1e10, -0.27900, 4)
    241. 

  241. 

  242.     def test_getMomentOfInertiaTensor(self):
    243. 

  243.         """
    244. 

  244.         Test the Conformer.getMomentOfInertiaTensor() method.
    245. 

  245.         """
    246. 

  246.         I = self.conformer.getMomentOfInertiaTensor()
    247. 

  247.         self.assertAlmostEqual(I[0,0]*constants.Na*1e23, 20.65968, 4)
    248. 

  248.         self.assertAlmostEqual(I[0,1]*constants.Na*1e23, -7.48115, 4)
    249. 

  249.         self.assertAlmostEqual(I[0,2]*constants.Na*1e23, -3.46416, 4)
    250. 

  250.         self.assertAlmostEqual(I[1,0]*constants.Na*1e23, -7.48115, 4)
    251. 

  251.         self.assertAlmostEqual(I[1,1]*constants.Na*1e23, 13.53472, 4)
    252. 

  252.         self.assertAlmostEqual(I[1,2]*constants.Na*1e23, -5.48630, 4)
    253. 

  253.         self.assertAlmostEqual(I[2,0]*constants.Na*1e23, -3.46416, 4)
    254. 

  254.         self.assertAlmostEqual(I[2,1]*constants.Na*1e23, -5.48630, 4)
    255. 

  255.         self.assertAlmostEqual(I[2,2]*constants.Na*1e23, 22.84296, 4)
    256. 

  256. 

  257.     def test_getPrincipalMomentsOfInertia(self):
    258. 

  258.         """
    259. 

  259.         Test the Conformer.getPrincipalMomentsOfInertia() method.
    260. 

  260.         """
    261. 

  261.         I, V = self.conformer.getPrincipalMomentsOfInertia()
    262. 

  262.         self.assertAlmostEqual(I[0]*constants.Na*1e23,  6.27074, 4)
    263. 

  263.         self.assertAlmostEqual(I[1]*constants.Na*1e23, 25.38321, 3)
    264. 

  264.         self.assertAlmostEqual(I[2]*constants.Na*1e23, 25.38341, 3)
    265. 

  265.         print V
    266. 

  266.         self.assertAlmostEqual(V[0,0],  0.497140, 4)
    267. 

  267.         self.assertAlmostEqual(V[0,1], -0.610114, 4)
    268. 

  268.         self.assertAlmostEqual(V[0,2], -0.616938, 4)
    269. 

  269.         self.assertAlmostEqual(V[1,0],  0.787360, 4)
    270. 

  270.         self.assertAlmostEqual(V[1,1],  0.018454, 4)
    271. 

  271.         self.assertAlmostEqual(V[1,2],  0.616218, 4)
    272. 

  272.         self.assertAlmostEqual(V[2,0],  0.364578, 4)
    273. 

  273.         self.assertAlmostEqual(V[2,1],  0.792099, 4)
    274. 

  274.         self.assertAlmostEqual(V[2,2], -0.489554, 4)
    275. 

  275. 

  276.     def test_getInternalReducedMomentOfInertia(self):
    277. 

  277.         """
    278. 

  278.         Test the Conformer.getInternalReducedMomentOfInertia() method.
    279. 

  279.         """
    280. 

  280.         I = self.conformer.getInternalReducedMomentOfInertia(pivots=[1,5], top1=[1,2,3,4])
    281. 

  281.         self.assertAlmostEqual(I*constants.Na*1e23, 1.56768, 4)
    282. 

  282.