Fix kinetics fitting with Tmin & Tmax in Arakne

arkane/kinetics.py

@@ -70,35 +70,19 @@ def __init__(self, reaction,
                  Tcount=0,
                  sensitivity_conditions=None):
         self.usedTST = False
-        if Tmin is not None:
-            self.Tmin = quantity.Quantity(Tmin)
-        else:
-            self.Tmin = None
-
-        if Tmax is not None:
-            self.Tmax = quantity.Quantity(Tmax)
-        else:
-            self.Tmax = None
-
-        self.Tcount = Tcount
+        self.Tmin = Tmin if Tmin is not None else (298, 'K')
+        self.Tmax = Tmax if Tmax is not None else (2500, 'K')
+        self.Tcount = Tcount if Tcount > 3 else 50
 
         if Tlist is not None:
-            self.Tlist = quantity.Quantity(Tlist)
-            self.Tmin = quantity.Quantity(numpy.min(self.Tlist.value_si), "K")
-            self.Tmax = quantity.Quantity(numpy.max(self.Tlist.value_si), "K")
+            self.Tlist = Tlist
+            self.Tmin = (min(self.Tlist.value_si), 'K')
+            self.Tmax = (max(self.Tlist.value_si), 'K')
             self.Tcount = len(self.Tlist.value_si)
         else:
-            if Tmin and Tmax is not None:
-
-                if self.Tcount <= 3.:
-                    self.Tcount = 50
-
-                stepsize = (self.Tmax.value_si - self.Tmin.value_si) / self.Tcount
-
-                self.Tlist = quantity.Quantity(numpy.arange(self.Tmin.value_si,
-                                                            self.Tmax.value_si + stepsize, stepsize), 'K')
-            else:
-                self.Tlist = None
+            self.Tlist = (1 / numpy.linspace(1 / self.Tmax.value_si,
+                                            1 / self.Tmin.value_si,
+                                            self.Tcount), 'K')
 
         self.reaction = reaction
         self.kunits = None
@@ -145,10 +129,7 @@ def execute(self, output_directory=None, plot=True):
         If `plot` is True, then plots of the raw and fitted values for the kinetics
         will be saved.
         """
-        if self.Tlist is not None:
-            self.generateKinetics(self.Tlist.value_si)
-        else:
-            self.generateKinetics()
+        self.generateKinetics()
         if output_directory is not None:
             try:
                 self.write_output(output_directory)
@@ -176,7 +157,7 @@ def execute(self, output_directory=None, plot=True):
         logging.debug('Finished kinetics job for reaction {0}.'.format(self.reaction))
         logging.debug(repr(self.reaction))
 
-    def generateKinetics(self, Tlist=None):
+    def generateKinetics(self):
         """
         Generate the kinetics data for the reaction and fit it to a modified Arrhenius model.
         """
@@ -203,19 +184,16 @@ def generateKinetics(self, Tlist=None):
             else:
                 raise ValueError('Unknown tunneling model {0!r} for reaction {1}.'.format(tunneling, self.reaction))
         logging.debug('Generating {0} kinetics model for {1}...'.format(kineticsClass, self.reaction))
-        if Tlist is None:
-            Tlist = 1000.0 / numpy.arange(0.4, 3.35, 0.05)
-        klist = numpy.zeros_like(Tlist)
-        for i in range(Tlist.shape[0]):
-            klist[i] = self.reaction.calculateTSTRateCoefficient(Tlist[i])
-
+        klist = numpy.zeros_like(self.Tlist.value_si)
+        for i, t in enumerate(self.Tlist.value_si):
+            klist[i] = self.reaction.calculateTSTRateCoefficient(t)
         order = len(self.reaction.reactants)
         klist *= 1e6 ** (order - 1)
         self.kunits = {1: 's^-1', 2: 'cm^3/(mol*s)', 3: 'cm^6/(mol^2*s)'}[order]
         self.Kequnits = {2: 'mol^2/cm^6', 1: 'mol/cm^3', 0: '       ', -1: 'cm^3/mol', -2: 'cm^6/mol^2'}[
             len(self.reaction.products) - len(self.reaction.reactants)]
         self.krunits = {1: 's^-1', 2: 'cm^3/(mol*s)', 3: 'cm^6/(mol^2*s)'}[len(self.reaction.products)]
-        self.reaction.kinetics = Arrhenius().fitToData(Tlist, klist, kunits=self.kunits)
+        self.reaction.kinetics = Arrhenius().fitToData(self.Tlist.value_si, klist, kunits=self.kunits)
         self.reaction.elementary_high_p = True
 
     def write_output(self, output_directory):

arkane/kineticsTest.py

@@ -0,0 +1,94 @@
+#!/usr/bin/env python
+# -*- coding: utf-8 -*-
+
+###############################################################################
+#                                                                             #
+# RMG - Reaction Mechanism Generator                                          #
+#                                                                             #
+# Copyright (c) 2002-2019 Prof. William H. Green (whgreen@mit.edu),           #
+# Prof. Richard H. West (r.west@neu.edu) and the RMG Team (rmg_dev@mit.edu)   #
+#                                                                             #
+# Permission is hereby granted, free of charge, to any person obtaining a     #
+# copy of this software and associated documentation files (the 'Software'),  #
+# to deal in the Software without restriction, including without limitation   #
+# the rights to use, copy, modify, merge, publish, distribute, sublicense,    #
+# and/or sell copies of the Software, and to permit persons to whom the       #
+# Software is furnished to do so, subject to the following conditions:        #
+#                                                                             #
+# The above copyright notice and this permission notice shall be included in  #
+# all copies or substantial portions of the Software.                         #
+#                                                                             #
+# THE SOFTWARE IS PROVIDED 'AS IS', WITHOUT WARRANTY OF ANY KIND, EXPRESS OR  #
+# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,    #
+# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE #
+# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER      #
+# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING     #
+# FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER         #
+# DEALINGS IN THE SOFTWARE.                                                   #
+#                                                                             #
+###############################################################################
+
+import unittest
+
+from rmgpy.species import TransitionState
+from rmgpy.reaction import Reaction
+
+from arkane.kinetics import KineticsJob
+
+################################################################################
+
+
+class KineticsTest(unittest.TestCase):
+    """
+    Contains unit tests for the Arkane Kinetics module
+    """
+
+    def test_give_tlist_for_kineticsjob(self):
+        """
+        Ensures that the proper temperature ranges are set when Tlist is specified
+        """
+        rxn = Reaction(transitionState=TransitionState())
+        Tlist = [50.7, 100, 300, 800, 1255]
+        kjob = KineticsJob(rxn, Tlist=(Tlist, 'K'))
+        self.assertEqual(min(Tlist), kjob.Tmin.value_si)
+        self.assertEqual(max(Tlist), kjob.Tmax.value_si)
+        self.assertEqual(len(Tlist), kjob.Tcount)
+
+    def test_give_Trange_for_kineticsjob(self):
+        """
+        Ensures that Tlist is set when a range of temperatures is specified
+        """
+        rxn = Reaction(transitionState=TransitionState())
+        kjob = KineticsJob(rxn, Tmin=(50, 'K'), Tmax=(4000, 'K'), Tcount=5)
+        self.assertEqual(5, len(kjob.Tlist.value_si))
+        self.assertEqual(50, min(kjob.Tlist.value_si))
+        self.assertAlmostEqual(4000, max(kjob.Tlist.value_si))
+        inverse_tlist = 1 / kjob.Tlist.value_si
+        self.assertAlmostEqual(inverse_tlist[1] - inverse_tlist[0],
+                               inverse_tlist[-1] - inverse_tlist[-2],
+                               msg='The points for fitting do not appear 1/T spaced. '
+                                   'Obtained values of {0} and {1}'.format(inverse_tlist[1] - inverse_tlist[0],
+                                                                           inverse_tlist[-1] - inverse_tlist[-2]))
+
+    def test_get_tlist_for_kineticsjob(self):
+        """
+        Ensures that Tlist is set when no range is specified
+        """
+        rxn = Reaction(transitionState=TransitionState())
+        kjob = KineticsJob(rxn)
+        self.assertAlmostEqual(298, kjob.Tmin.value_si)
+        self.assertAlmostEqual(2500, kjob.Tmax.value_si)
+        self.assertEqual(50, kjob.Tcount)
+        self.assertEqual(50, len(kjob.Tlist.value_si))
+        self.assertAlmostEqual(298, min(kjob.Tlist.value_si))
+        self.assertAlmostEqual(2500, max(kjob.Tlist.value_si))
+        inverse_tlist = 1 / kjob.Tlist.value_si
+        self.assertAlmostEqual(inverse_tlist[1] - inverse_tlist[0],
+                               inverse_tlist[-1] - inverse_tlist[-2],
+                               msg='The points for fitting do not appear 1/T spaced. '
+                                   'Obtained values of {0} and {1}'.format(inverse_tlist[1] - inverse_tlist[0],
+                                                                           inverse_tlist[-1] - inverse_tlist[-2]))
+
+
+if __name__ == '__main__':
+    unittest.main(testRunner=unittest.TextTestRunner(verbosity=2))

documentation/source/users/arkane/input.rst

@@ -818,24 +818,28 @@ High Pressure Limit Kinetics
 
 Use a ``kinetics()`` function to make Arkane execute the high-pressure limit kinetic
 parameters computation for a reaction. The ``'label'`` string must correspond to that of
-a defined ``reaction()`` function. If desired, define a temperature range and number
-of temperatures at which the high-pressure rate coefficient will be tabulated and saved to
-the outupt file. The 3-parameter modified Arrhenius coefficients will automatically be fit
-to the computed rate coefficients. The quantum tunneling factor will also be displayed.
+a defined ``reaction()`` function.
 
-Below is a typical ``kinetics()`` function::
+You have three options for specifying the temperature to which a modified Arrhenius
+expression will be fit.
+
+Give an explicit list of temperatures to fit::
 
     kinetics(
     label = 'H + C2H4 <=> C2H5',
-    Tmin = (400,'K'), Tmax = (1200,'K'), Tcount = 6,
+    Tlist = ([400,500,700,900,1100,1200],'K'),
     )
 
-If specific temperatures are desired, you may specify a list
-(``Tlist = ([400,500,700,900,1100,1200],'K')``) instead of Tmin, Tmax, and Tcount.
+Give minimmum and maximum temperatures to fit::
+
+    kinetics(
+    label = 'H + C2H4 <=> C2H5',
+    Tmin = (400,'K'), Tmax = (1200,'K'), Tcount = 6,
+    )
 
-This is also acceptable::
+Use the default range to fit (298 K to 2500 K at 50 points spaced equally in 1/T space)::
 
-    kinetics('H + C2H4 <=> C2H5')
+    kinetics(label = 'H + C2H4 <=> C2H5')
 
 If a sensitivity analysis is desired, simply add the conditions at which to calculate sensitivity coefficients
 in the following format, e.g.::
@@ -846,6 +850,9 @@ in the following format, e.g.::
         sensitivity_conditions = [(1000, 'K'), (2000, 'K')]
     )
 
+The ``output.py`` file will show the rates at various temperatures including the quantum tunneling factor.
+Kinetic rates will also be added to the ``chem.inp`` file.
+
 The output of a sensitivity analysis is saved into a ``sensitivity`` folder in the output directory. A text file, named
 with the reaction label, delineates the semi-normalized sensitivity coefficients ``dln(k)/dE0`` in units of ``mol/J``
 at all requested conditions. A horizontal bar figure is automatically generated per reaction with subplots for both the