# -*- coding: utf-8 -*-
"""
Contains functions for the `Arrhenius equation
<https://en.wikipedia.org/wiki/Arrhenius_equation>`_
(:func:`arrhenius_equation`) and a convenience fitting routine
(:func:`fit_arrhenius_equation`).
"""
from __future__ import (absolute_import, division, print_function)
from .._util import get_backend
from ..util.regression import least_squares
from ..util.pyutil import defaultnamedtuple
from ..units import default_constants, default_units, format_string
try:
import numpy as np
except ImportError:
np = None
def _fit_linearized(backtransfm, lin_x, lin_y, lin_yerr):
if len(lin_x) != len(lin_y):
raise ValueError("k and T needs to be of equal length.")
if lin_yerr is not None:
if len(lin_yerr) != len(lin_y):
raise ValueError("kerr and T needs to be of equal length.")
lin_p, lin_vcv, lin_r2 = least_squares(lin_x, lin_y, lin_yerr)
return [cb(lin_p) for cb in backtransfm]
def _fit(T, k, kerr, func, lin_x, lin_y, backtransfm, linearized=False):
_lin_y = lin_y(T, k)
if kerr is None:
lin_yerr = 1
else:
lin_yerr = (abs(lin_y(k - kerr, T) - _lin_y) +
abs(lin_y(k + kerr, T) - _lin_y))/2
lopt = _fit_linearized(backtransfm, lin_x(T, k), _lin_y, lin_yerr)
if linearized:
return lopt
from scipy.optimize import curve_fit
popt, pcov = curve_fit(func, T, k, lopt, kerr)
return popt, pcov
def _get_R(constants=None, units=None):
if constants is None:
R = 8.314472
if units is not None:
J = units.Joule
K = units.Kelvin
mol = units.mol
R *= J/mol/K
else:
R = constants.molar_gas_constant
return R
[docs]def arrhenius_equation(A, Ea, T, constants=None, units=None, backend=None):
"""
Returns the rate coefficient according to the Arrhenius equation
Parameters
----------
A: float with unit
frequency factor
Ea: float with unit
activation energy
T: float with unit
temperature
constants: object (optional, default: None)
if None:
T assumed to be in Kelvin, Ea in J/(K mol)
else:
attributes accessed: molar_gas_constant
Tip: pass quantities.constants
units: object (optional, default: None)
attributes accessed: Joule, Kelvin and mol
backend: module (optional)
module with "exp", default: numpy, math
"""
be = get_backend(backend)
R = _get_R(constants, units)
try:
RT = (R*T).rescale(Ea.dimensionality)
except AttributeError:
RT = R*T
return A*be.exp(-Ea/RT)
[docs]def fit_arrhenius_equation(T, k, kerr=None, linearized=False, constants=None, units=None):
""" Curve fitting of the Arrhenius equation to data points
Parameters
----------
T : float
k : array_like
kerr : array_like (optional)
linearized : bool
"""
return _fit(T, k, kerr, arrhenius_equation, lambda T, k: 1/T, lambda T, k: np.log(k),
[lambda p: np.exp(p[0]), lambda p: -p[1]*_get_R(constants, units)], linearized=linearized)
def _fit_arrhenius_equation(T, k, kerr=None, linearized=False):
""" Curve fitting of the Arrhenius equation to data points
Parameters
----------
k : array_like
T : float
kerr : array_like (optional)
linearized : bool
"""
if len(k) != len(T):
raise ValueError("k and T needs to be of equal length.")
from math import exp
import numpy as np
p = np.polyfit(1/T, np.log(k), 1)
R = _get_R(constants=None, units=None)
Ea = -R*p[0]
A = exp(p[1])
if linearized:
return A, Ea
from scipy.optimize import curve_fit
if kerr is None:
weights = None
else:
weights = 1/kerr**2
popt, pcov = curve_fit(arrhenius_equation, T, k, [A, Ea], weights)
return popt, pcov
[docs]class ArrheniusParam(defaultnamedtuple('ArrheniusParam', 'A Ea ref', [None])):
""" Kinetic data in the form of an Arrhenius parameterisation
Parameters
----------
Ea: float
activation energy
A: float
preexponential prefactor (Arrhenius type eq.)
ref: object (default: None)
arbitrary reference (e.g. string representing citation key)
Examples
--------
>>> k = ArrheniusParam(1e13, 40e3)
>>> '%.5g' % k(298.15)
'9.8245e+05'
"""
@classmethod
[docs] def from_fit_of_data(cls, T, k, kerr=None, **kwargs):
args, vcv = fit_arrhenius_equation(T, k, kerr)
return cls(*args, **kwargs)
def __call__(self, T, constants=None, units=None, backend=None):
""" Evaluates the arrhenius equation for a specified state
Parameters
----------
T: float
constants: module (optional)
units: module (optional)
backend: module (default: math)
See also
--------
See :func:`chempy.arrhenius.arrhenius_equation`.
"""
return arrhenius_equation(self.A, self.Ea, T, constants=constants,
units=units, backend=backend)
def _as_RateExpr(self, rxn, arg_keys=None, constants=None, units=None, backend=None):
from .rates import ArrheniusMassAction as AMA
return AMA([self.A, self.Ea/_get_R(constants, units)], arg_keys, rxn=rxn, ref=self.ref)
[docs] def equation_as_string(self, precision, tex=False):
(str_A, str_A_unit), (str_Ea, str_Ea_unit) = self.format(precision, tex)
if tex:
return r"{}\exp \left(-\frac{{{}}}{{RT}} \right)".format(
str_A, str_Ea + ' ' + str_Ea_unit), str_A_unit
else:
return "{}*exp(-{}/(R*T))".format(str_A, str_Ea + ' ' + str_Ea_unit), str_A_unit
def __str__(self):
return ' '.join(self.equation_as_string('%.5g'))
[docs]class ArrheniusParamWithUnits(ArrheniusParam):
def __call__(self, state, constants=default_constants, units=default_units,
backend=None):
""" See :func:`chempy.arrhenius.arrhenius_equation`. """
return super(ArrheniusParamWithUnits, self).__call__(
state, constants, units, backend)
def _as_RateExpr(self, rxn, arg_keys=None, constants=default_constants, units=default_units):
return super(ArrheniusParamWithUnits, self)._as_RateExpr(rxn, arg_keys, constants, units)