Skip to content

Commit

Permalink
Set Tlist when creating KineticsJob
Browse files Browse the repository at this point in the history
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
  • Loading branch information
goldmanm committed Jul 30, 2019
1 parent 69122e6 commit 46ae3d7
Show file tree
Hide file tree
Showing 3 changed files with 126 additions and 47 deletions.
52 changes: 15 additions & 37 deletions arkane/kinetics.py
Expand Up @@ -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
Expand Down Expand Up @@ -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)
Expand Down Expand Up @@ -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.
"""
Expand All @@ -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):
Expand Down
94 changes: 94 additions & 0 deletions 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))
27 changes: 17 additions & 10 deletions documentation/source/users/arkane/input.rst
Expand Up @@ -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.::
Expand All @@ -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
Expand Down

0 comments on commit 46ae3d7

Please sign in to comment.