From 46ae3d705066cbd5096ffb6f2ac1337ce40d896e Mon Sep 17 00:00:00 2001 From: MarkGoldman Date: Sun, 28 Jul 2019 01:06:08 -0400 Subject: [PATCH] Set Tlist when creating KineticsJob Previously, there were two spots where a list of temperatures was created for fitting high pressure kinetics in Arkane. One fit using inverse temperature spacing and the other one fit with linear temperature spacing. See #1671 This commit moves the location for creating a listing of temperatures to the __init__ method, so there is only one spot where temperatures are initialized. The default min and max in `generateKinetics` was moved into the __init__ method. Three unittests were written for obtaining the proper temperatures in the __init__ method. `generateKinetics` previously accepted an arguement Tlist, which is already an attribute of KineticsJob, creating redundancy. The Tlist argument was removed and self.Tlist was used in this method instead, simplifying the code. Calls to quantity.Quantity in the __init__ method were unnecessary before setting them with quantity.Temperature, so the __init__ code was simplified by setting attributes with (value, 'K') instead. Documentation was updated to better describe the three methods of specifying temperature when calling `kinetics` in Arkane --- arkane/kinetics.py | 52 ++++-------- arkane/kineticsTest.py | 94 +++++++++++++++++++++ documentation/source/users/arkane/input.rst | 27 +++--- 3 files changed, 126 insertions(+), 47 deletions(-) create mode 100644 arkane/kineticsTest.py diff --git a/arkane/kinetics.py b/arkane/kinetics.py index 57adfab277..80ca56cd4a 100644 --- a/arkane/kinetics.py +++ b/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): diff --git a/arkane/kineticsTest.py b/arkane/kineticsTest.py new file mode 100644 index 0000000000..cf52745e21 --- /dev/null +++ b/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)) diff --git a/documentation/source/users/arkane/input.rst b/documentation/source/users/arkane/input.rst index ce6f0a996f..93977a8792 100644 --- a/documentation/source/users/arkane/input.rst +++ b/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