Skip to content

Commit

Permalink
Merge pull request #1672 from ReactionMechanismGenerator/fix_kinetic_…
Browse files Browse the repository at this point in the history
…fitting_arkane

Fix kinetics fitting with Tmin & Tmax in Arakne
  • Loading branch information
alongd authored Aug 1, 2019
2 parents 69122e6 + 46ae3d7 commit 589ba9a
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
Original file line number Diff line number Diff line change
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
Original file line number Diff line number Diff line change
@@ -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
Original file line number Diff line number Diff line change
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 589ba9a

Please sign in to comment.