--- \n\n\n\n\n
![\"overreact\"](\"https://raw.githubusercontent.com/geem-lab/overreact-guide/master/logo.png\")
\n
---
\n\noverreact is a library and a command-line tool for building and\nanalyzing homogeneous microkinetic models from first-principles\ncalculations:
\n\n\n
In [1]: from overreact import api # the api\n\nIn [2]: api.get_k("S -> E\u2021 -> S", # your model\n ...: {"S": "data/ethane/B97-3c/staggered.out", # your data\n ...: "E\u2021": "data/ethane/B97-3c/eclipsed.out"})\nOut[2]: array([8.16880917e+10]) # your results\n
\n
The user specifies a set of\nelementary reactions that are believed to be relevant for the overall chemical\nphenomena. overreact offers a hopefully complete but simple environment for\nhypothesis testing in first-principles chemical kinetics.
\n\n\n \n \ud83e\udd14 What is microkinetic modeling?\n
\n \n Microkinetic modeling is a technique used to predict the outcome\n of complex chemical reactions.\n It can be used\n to investigate the catalytic transformations\n of molecules.\n overreact makes it easy to create\n and analyze microkinetic models built\n from computational chemistry data.\n
\n
\n\n
\n\n\n \n \ud83e\uddd0 What do you mean by first-principles calculations?\n
\n \n We use the term first-principles calculations to refer to\n calculations performed using quantum chemical modern methods such as\n Wavefunction\n and\n Density Functional\n theories.\n For instance, the three-line example code above calculates the rate of methyl rotation in ethane (at\n B97-3c).\n (Rather surprisingly, the error found is less than 2%\n when compared to available experimental results.)\n
\n
\n\n
\n\noverreact uses precise thermochemical partition funtions, tunneling\ncorrections and data is parsed directly from computational chemistry\noutput files thanks to cclib (see the\nlist of its supported programs).
\n\nInstallation
\n\noverreact is a Python package, so you can easily install it with\npip:
\n\n\n
$ pip install "overreact[cli,fast]"\n
\n
See the\ninstallation guide\nfor more details.
\n\n\n \ud83d\ude80 Where to go from here? Take a look at the\n short introduction.\n Or see\n below\n for more guidance.
\n
\n\nCiting overreact
\n\nIf you use overreact in your research, please cite:
\n\n\n Schneider, F. S. S.; Caramori, G. F.\n _Overreact, an in Silico Lab: Automative Quantum Chemical Microkinetic Simulations for Complex Chemical Reactions_.\n Journal of Computational Chemistry 2022, 44 (3), 209\u2013217.\n doi:10.1002/jcc.26861.
\n
\n\nHere's the reference in BibTeX format:
\n\n\n
@article{overreact_paper2022,\n title = {Overreact, an in silico lab: Automative quantum chemical microkinetic simulations for complex chemical reactions},\n author = {Schneider, Felipe S. S. and Caramori, Giovanni F.},\n year = {2022},\n month = {Apr},\n journal = {Journal of Computational Chemistry},\n publisher = {Wiley},\n volume = {44},\n number = {3},\n pages = {209\u2013217},\n doi = {10.1002/jcc.26861},\n issn = {1096-987x},\n url = {http://dx.doi.org/10.1002/jcc.26861},\n}\n@software{overreact_software2021,\n title = {geem-lab/overreact: v1.2.0 \\vert{} Zenodo},\n author = {Felipe S. S. Schneider and Let\\'{\\i}cia M. P. Madureira and Giovanni F. Caramori},\n year = {2023},\n month = {Jan},\n publisher = {Zenodo},\n doi = {10.5281/zenodo.7504800},\n url = {https://doi.org/10.5281/zenodo.7504800},\n version = {v1.2.0},\n howpublished = {\\url{https://github.com/geem-lab/overreact}},\n}\n
\n
License
\n\noverreact is open-source, released under the permissive MIT license. See\nthe LICENSE agreement.
\n\nFunding
\n\nThis project was developed at the GEEM lab\n(Federal University of Santa Catarina, Brazil), and was\npartially funded by the\nBrazilian National Council for Scientific and Technological Development (CNPq),\ngrant number 140485/2017-1.
\n"}, "overreact.Scheme": {"fullname": "overreact.Scheme", "modulename": "overreact", "qualname": "Scheme", "kind": "class", "doc": "A descriptor of a chemical reaction network.
\n\nMostly likely, this comes from a parsed input file.\nSee overreact.io.parse_model.
\n", "bases": "typing.NamedTuple"}, "overreact.Scheme.__init__": {"fullname": "overreact.Scheme.__init__", "modulename": "overreact", "qualname": "Scheme.__init__", "kind": "function", "doc": "Create new instance of Scheme(compounds, reactions, is_half_equilibrium, A, B)
\n", "signature": "(\tcompounds: Sequence[str],\treactions: Sequence[str],\tis_half_equilibrium: Sequence[bool],\tA: numpy.ndarray,\tB: numpy.ndarray)"}, "overreact.Scheme.compounds": {"fullname": "overreact.Scheme.compounds", "modulename": "overreact", "qualname": "Scheme.compounds", "kind": "variable", "doc": "A descriptor of compounds.
\n", "annotation": ": Sequence[str]"}, "overreact.Scheme.reactions": {"fullname": "overreact.Scheme.reactions", "modulename": "overreact", "qualname": "Scheme.reactions", "kind": "variable", "doc": "A descriptor of reactions.
\n", "annotation": ": Sequence[str]"}, "overreact.Scheme.is_half_equilibrium": {"fullname": "overreact.Scheme.is_half_equilibrium", "modulename": "overreact", "qualname": "Scheme.is_half_equilibrium", "kind": "variable", "doc": "An indicator of whether a reaction is half-equilibrium.
\n", "annotation": ": Sequence[bool]"}, "overreact.Scheme.A": {"fullname": "overreact.Scheme.A", "modulename": "overreact", "qualname": "Scheme.A", "kind": "variable", "doc": "A matrix of stoichiometric coefficients between reactants and products.
\n", "annotation": ": numpy.ndarray"}, "overreact.Scheme.B": {"fullname": "overreact.Scheme.B", "modulename": "overreact", "qualname": "Scheme.B", "kind": "variable", "doc": "A matrix of stoichiometric coefficients between reactants and transition states.
\n", "annotation": ": numpy.ndarray"}, "overreact.change_reference_state": {"fullname": "overreact.change_reference_state", "modulename": "overreact", "qualname": "change_reference_state", "kind": "function", "doc": "Calculate an additive entropy correction to a change in reference states.
\n\n$$\\Delta G_\\text{corr} =\n R T \\ln \\left( \\frac{\\chi_\\text{new}}{\\chi_\\text{old}} \\right)$$
\n\nThe value returned can be directly multiplied by temperature and summed to\nthe old reference free energies to obtain free energies with respect to a\nnew reference. See notes below.
\n\nFor instance, the concentration correction to Gibbs free energy for a\ngas-to-liquid standard state change is simply\n(\\( c^\\circ = \\frac{\\text{1 atm}}{R T} \\)),
\n\n$$\\Delta G_\\text{conc} =\n R T \\ln \\left( \\frac{\\text{1 M}}{c^\\circ} \\right)$$
\n\nParameters
\n\nnew_reference : array-like, optional\n New reference state. Default value corresponds to 1 mol/liter.\nold_reference : array-like, optional\n Old reference state. Default value corresponds to the concentration of\n an ideal gas at the given temperature and 1 atm.\nsign : float, optional\n Sign of the change in reference state. Default value is 1. This only\n multiplies the final result.\ntemperature : array-like, optional\n Absolute temperature in Kelvin.\npressure : array-like, optional\n Reference gas pressure.\nvolume : float, optional\n Molar volume.
\n\nReturns
\n\ncorrection : array-like\n Entropy correction in J/mol\u00b7K.
\n\nNotes
\n\nThis function can be used to add any entropy correction in the form above.\nThe only drawback is that, sometimes, those corrections are written with a\nminus sign in front of them (this implies switching the roles of\nold_reference and new_reference). The easiest way to accomplish this is\nby using sign=-1 or multiplying the result by -1.
\n\nExamples
\n\nBy default, the correction returns a change in concentration from the gas\nphase standard concentration to the solvated-state standard concentration:
\n\n\n
>>> -rx.change_reference_state() / constants.calorie\n-6.4\n>>> 298.15 * rx.change_reference_state() / constants.kcal\n1.89\n>>> 273.15 * rx.change_reference_state(temperature=273.15) / constants.kcal\n1.69\n
\n
But this function can also be used to adjust symmetry effects from C1\ncalculations (symmetry number equals to one). For D7h, for instance, the\nsymmetry number is 14:
\n\n\n
>>> -298.15 * rx.change_reference_state(14, 1) / constants.kcal\n-1.56\n
\n
\n
>>> rx.change_reference_state(sign=-1) == -rx.change_reference_state()\nTrue\n
\n
Estimate a energy bias for a given set of reference data points.
\n\nParameters
\n\nscheme : Scheme\n A descriptor of the reaction scheme.\n Mostly likely, this comes from a parsed model input file.\n See overreact.io.parse_model.\ncompounds : dict-like\n A descriptor of the compounds.\n Mostly likely, this comes from a parsed model input file.\n See overreact.io.parse_model.\ndata : dict-like of array-like\ny0: array-like\ntunneling : str or None, optional\n Choose between \"eckart\", \"wigner\" or None (or \"none\").\nqrrho : bool or tuple-like, optional\n Apply both the quasi-rigid rotor harmonic oscillator (QRRHO)\n approximations of M. Head-Gordon and others (enthalpy correction, see\n J. Phys. Chem. C 2015, 119, 4, 1840-1850) and S. Grimme (entropy correction, see\n Theory. Chem. Eur. J., 2012, 18: 9955-9964) on top of the classical RRHO.\ntemperature : array-like, optional\n Absolute temperature in Kelvin.\npressure : array-like, optional\n Reference gas pressure.\ndelta_freeenergies : array-like, optional\n Use this instead of obtaining delta free energies from the compounds.\nmolecularity : array-like, optional\n Reaction order, i.e., number of molecules that come together to react.\n If set, this is used to calculate delta_moles for\n equilibrium_constant, which effectively calculates a solution\n equilibrium constant between reactants and the transition state for\n gas phase data. You should set this to None if your free energies\n were already adjusted for solution Gibbs free energies.\nvolume : float, optional\n Molar volume.
\n\nReturns
\n\narray-like
\n\nExamples
\n\n\n
>>> model = rx.parse_model("data/tanaka1996/UMP2/cc-pVTZ/model.jk")\n
\n
The following are some estimates on actual atmospheric concentrations:
\n\n\n
>>> y0 = [4.8120675684099e-5,\n... 2.8206357713029e-5,\n... 0.0,\n... 0.0,\n... 2.7426565371219e-5]\n>>> data = {"t": [1.276472128376942246e-6,\n... 1.446535794555581743e-4,\n... 1.717069678525567564e-2],\n... "CH3\u00b7": [9.694916853338366211e-9,\n... 1.066033349343709026e-6,\n... 2.632179124780495175e-5]}\n>>> get_bias(model.scheme, model.compounds, data, y0) / constants.kcal\n-1.4\n
\n
Calculate deltas according to reactions.
\n\nDelta properties are differences in a property between the final and\ninitial state of a chemical transformation. They are calculated from\nmatrices representing the transformation and the absolute properties.\nTransformation matrices are expected to have column-wise transformation\ndefined.
\n\nVery useful for the calculation of reaction and activation free energies\nfrom absolute free energies of compounds. Matrices A and B of a\nScheme represent the transformations associated with reaction and\nactivation free energies, respectively.
\n\nParameters
\n\ntransform : array-like\nproperty : array-like
\n\nReturns
\n\ndelta_property : array-like
\n\nExamples
\n\n\n
>>> get_delta([-1, 1], [-10, 10])\n20\n
\n
You must ensure the transformation is properly defined, as no test is made\nto ensure, e.g., conservation of matter:
\n\n\n
>>> get_delta([-1, 0], [-10, 20])\n10\n
\n
Normally, transformations are given as columns in a matrix:
\n\n\n
>>> get_delta([[-1, -2],\n... [ 1, 3]], [-5, 12])\narray([17, 46])\n
\n
Generate a rate function that models a reaction scheme.
\n\nParameters
\n\nscheme : Scheme\n A descriptor of the reaction scheme.\n Mostly likely, this comes from a parsed model input file.\n See overreact.io.parse_model.\nk : array-like\n Reaction rate constant(s). Units match the concentration units given to\n the returned function dydt.\nef : float, optional\n Equilibrium factor. This is a parameter that can be used to scale the\n reaction rates associated to half-equilibrium reactions such that they\n are faster than the other reactions.
\n\nReturns
\n\ndydt : callable\n Reaction rate function. The actual reaction rate constants employed\n are stored in the attribute k of the returned function. If JAX is\n available, the attribute jac will hold the Jacobian function of\n dydt.
\n\nNotes
\n\nThe returned function is suited to be used by ODE solvers such as\nscipy.integrate.solve_ivp or the older scipy.integrate.ode (see\nexamples below). This is actually what the function get_y from the\ncurrent module does.
\n\nExamples
\n\n\n
>>> import numpy as np\n>>> import overreact as rx\n
\n
\n
>>> scheme = rx.parse_reactions("A <=> B")\n>>> dydt = get_dydt(scheme, np.array([1, 1]))\n>>> dydt(0.0, np.array([1., 1.]))\nArray([0., 0.], ...)\n
\n
If available, JAX is used for JIT compilation. This will make dydt\ncomplain if given lists instead of numpy arrays. So stick to the safer,\nfaster side as above.
\n\nThe actually used reaction rate constants can be inspected with the k\nattribute of dydt:
\n\n\n
>>> dydt.k\nArray([1., 1.], ...)\n
\n
If JAX is available, the Jacobian function will be available as\ndydt.jac:
\n\n\n
>>> dydt.jac(0.0, np.array([1., 1.]))\nArray([[-1., 1.],\n [ 1., -1.]], ...)\n
\n
Obtain enthalpies for compounds at a given temperature.
\n\nParameters
\n\ncompounds : dict-like\n A descriptor of the compounds.\n Mostly likely, this comes from a parsed input file.\n See overreact.io.parse_model.\nqrrho : bool, optional\n Apply the quasi-rigid rotor harmonic oscillator (QRRHO) approximation of\n M. Head-Gordon and others (see\n J. Phys. Chem. C 2015, 119, 4, 1840-1850)\n on top of the classical RRHO.\ntemperature : array-like, optional\n Absolute temperature in Kelvin.
\n\nReturns
\n\narray-like
\n\nExamples
\n\n\n
>>> import overreact as rx\n>>> from overreact import _constants as constants\n>>> model = rx.parse_model("data/ethane/B97-3c/model.k")\n>>> enthalpies = get_enthalpies(model.compounds)\n>>> (enthalpies - enthalpies.min()) / constants.kcal\narray([0. , 2.20053981])\n
\n
The enthalpies at absolute zero can easily be obtained (this is used,\ne.g., in the calculation of the Eckart tunneling coefficient, see\noverreact.tunnel.eckart). We can use this to calculate, for instance,\nthe thermal contributions to the enthalpy:
\n\n\n
>>> zero_enthalpies = get_enthalpies(model.compounds, temperature=0)\n>>> (enthalpies - zero_enthalpies) / constants.kcal\narray([2.78, 2.50])\n
\n
Obtain entropies for compounds at a given temperature and pressure.
\n\nParameters
\n\ncompounds : dict-like\n A descriptor of the compounds.\n Mostly likely, this comes from a parsed input file.\n See overreact.io.parse_model.\nenvironment : str or None, optional\n Choose between \"gas\" and a solvent. This is chosen for you by default,\n based on the names of each compound (e.g. A(g) or A is gas,\n A(w) or A(...) is solvated). In case this is given, all compounds\n will have the same behavior.\nmethod : str, optional\n This is a placeholder for future functionality.\n There are plans to implement more sophisticated methods for calculating\n entropies such as in\n Phys. Chem. Chem. Phys., 2019, 21, 18920-18929\n and\n J. Chem. Theory Comput. 2019, 15, 5, 3204-3214.\n Head over to the\n discussions if\n you're interested and would like to contribute.\n Leave this as \"standard\" for now.\nqrrho : bool, optional\n Apply the quasi-rigid rotor harmonic oscillator (QRRHO) approximation of\n S. Grimme (see\n Theory. Chem. Eur. J., 2012, 18: 9955-9964)\n on top of the classical RRHO.\ntemperature : array-like, optional\n Absolute temperature in Kelvin.\npressure : array-like, optional\n Reference gas pressure.
\n\nReturns
\n\narray-like
\n\nExamples
\n\n\n
>>> import overreact as rx\n>>> from overreact import _constants as constants\n>>> model = rx.parse_model("data/ethane/B97-3c/model.k")\n>>> entropies = get_entropies(model.compounds)\n>>> (entropies - entropies.min()) / constants.calorie\narray([1.4, 0. ])\n
\n
You can consider all compounds as solvated if you want:
\n\n\n
>>> sol_entropies = get_entropies(model.compounds, environment="solvent")\n>>> (sol_entropies - entropies) / constants.calorie\narray([-6.35360874, -6.35360874])\n
\n
Generate an alternative scheme with some concentrations fixed.
\n\nThis function returns data that allow the microkinetic simulation of a\nreaction network under constraints, namely when some compounds have fixed\nconcentrations. This works by 1. removing all references to the fixed\ncompounds and by 2. properly multiplying the reaction rate constants by\nthe respective concentrations.
\n\nParameters
\n\nscheme : Scheme\n A descriptor of the reaction scheme.\n Mostly likely, this comes from a parsed model input file.\n See overreact.io.parse_model.\nk : array-like\n Reaction rate constant(s). Units match the concentration units given to\n the returned function dydt.\nfixed_y0 : dict-like\n Fixed initial state. Units match the concentration units given to\n the returned function dydt.
\n\nReturns
\n\nscheme : Scheme\n Associated reaction scheme with all references to fixed compounds\n removed.\nk : array-like\n Associated (effective) reaction rate constants that model the fixed\n concentrations.
\n\nNotes
\n\nKeep in mind that when a compound get its concentration fixed, the\nreaction scheme no longer conserves matter. You can think of it as\nreacting close to an infinite source of the compound, but it accumulates\nin the milleu at the given concentration.
\n\nExamples
\n\n\n
>>> import numpy as np\n>>> import overreact as rx\n
\n
Equilibria under a specific pH can be easily modeled:
\n\n\n
>>> pH = 7\n>>> scheme = rx.parse_reactions("AH <=> A- + H+")\n>>> k = np.array([1, 1])\n>>> scheme, k = rx.get_fixed_scheme(scheme, k, {"H+": 10**-pH})\n>>> scheme\nScheme(compounds=('AH', 'A-'),\n reactions=('AH -> A-',\n 'A- -> AH'),\n is_half_equilibrium=(True, True),\n A=((-1.0, 1.0),\n (1.0, -1.0)),\n B=((-1.0, 0.0),\n (1.0, 0.0)))\n>>> k\narray([1.e+00, 1.e-07])\n
\n
It is also possible to model the fixed activity of a solvent, for\ninstance:
\n\n\n
>>> scheme = rx.parse_reactions("A + 2H2O -> B")\n>>> k = np.array([1.0])\n>>> scheme, k = rx.get_fixed_scheme(scheme, k, {"H2O": 55.6})\n>>> scheme\nScheme(compounds=('A', 'B'),\n reactions=('A -> B',),\n is_half_equilibrium=(False,),\n A=((-1.0,),\n (1.0,)),\n B=((-1.0,),\n (1.0,)))\n>>> k\narray([3091.36])\n
\n
Multiple reactions work fine, see both examples below:
\n\n\n
>>> pH = 12\n>>> scheme = rx.parse_reactions("B <- AH <=> A- + H+")\n>>> k = np.array([10.0, 1, 1])\n>>> scheme, k = rx.get_fixed_scheme(scheme, k, {"H+": 10**-pH})\n>>> scheme\nScheme(compounds=('AH', 'B', 'A-'),\n reactions=('AH -> B',\n 'AH -> A-',\n 'A- -> AH'),\n is_half_equilibrium=(False, True, True),\n A=((-1.0, -1.0, 1.0),\n (1.0, 0.0, 0.0),\n (0.0, 1.0, -1.0)),\n B=((-1.0, -1.0, 0.0),\n (1.0, 0.0, 0.0),\n (0.0, 1.0, 0.0)))\n>>> k\narray([1.e+01, 1.e+00, 1.e-12])\n
\n
\n
>>> pH = 2\n>>> scheme = rx.parse_reactions(["AH <=> A- + H+", "B- + H+ <=> BH"])\n>>> k = np.array([1, 1, 2, 2])\n>>> scheme, k = rx.get_fixed_scheme(scheme, k, {"H+": 10**-pH})\n>>> scheme\nScheme(compounds=('AH', 'A-', 'B-', 'BH'),\n reactions=('AH -> A-',\n 'A- -> AH',\n 'B- -> BH',\n 'BH -> B-'),\n is_half_equilibrium=(True, True, True, True),\n A=((-1.0, 1.0, 0.0, 0.0),\n (1.0, -1.0, 0.0, 0.0),\n (0.0, 0.0, -1.0, 1.0),\n (0.0, 0.0, 1.0, -1.0)),\n B=((-1.0, 0.0, 0.0, 0.0),\n (1.0, 0.0, 0.0, 0.0),\n (0.0, 0.0, -1.0, 0.0),\n (0.0, 0.0, 1.0, 0.0)))\n>>> k\narray([1. , 0.01, 0.02, 2. ])\n
\n
Multiple fixed compounds also work fine:
\n\n\n
>>> pH = 6\n>>> scheme = rx.parse_reactions("A + H2O -> B <=> B- + H+")\n>>> k = np.array([1.0, 100.0, 2.0])\n>>> scheme, k = rx.get_fixed_scheme(scheme, k, {"H+": 10**-pH, "H2O": 55.6})\n>>> scheme\nScheme(compounds=('A', 'B', 'B-'),\n reactions=('A -> B',\n 'B -> B-',\n 'B- -> B'),\n is_half_equilibrium=(False, True, True),\n A=((-1.0, 0.0, 0.0),\n (1.0, -1.0, 1.0),\n (0.0, 1.0, -1.0)),\n B=((-1.0, 0.0, 0.0),\n (1.0, -1.0, 0.0),\n (0.0, 1.0, 0.0)))\n>>> k\narray([5.56e+01, 1.00e+02, 2.00e-06])\n
\n
This function is a no-op if fixed_y0 is empty, which is very important\nfor overall code consistency:
\n\n\n
>>> scheme = rx.parse_reactions(["AH <=> A- + H+", "B- + H+ <=> BH"])\n>>> k = np.array([1, 1, 2, 2])\n>>> new_scheme, new_k = rx.get_fixed_scheme(scheme, k, {})\n>>> new_scheme == scheme\nTrue\n>>> np.allclose(new_k, k)\nTrue\n
\n
Obtain free energies for compounds at a given temperature and pressure.
\n\nParameters
\n\ncompounds : dict-like\n A descriptor of the compounds.\n Mostly likely, this comes from a parsed input file.\n See overreact.io.parse_model.\nbias : array-like, optional\n Energy to be added to free energies.\nenvironment : str or None, optional\n Choose between \"gas\" and a solvent. This is chosen for you by default,\n based on the names of each compound. If given, all compounds will\n have the same behavior.\nmethod : str, optional\n This is a placeholder for future functionality.\n There are plans to implement more sophisticated methods for calculating\n entropies such as in\n Phys. Chem. Chem. Phys., 2019, 21, 18920-18929\n and\n J. Chem. Theory Comput. 2019, 15, 5, 3204-3214.\n Head over to the\n discussions if\n you're interested and would like to contribute.\n Leave this as \"standard\" for now.\nqrrho : bool or tuple-like, optional\n Apply both the quasi-rigid rotor harmonic oscillator (QRRHO)\n approximations of M. Head-Gordon and others (enthalpy correction, see\n J. Phys. Chem. C 2015, 119, 4, 1840-1850)\n and S. Grimme (entropy correction, see\n Theory. Chem. Eur. J., 2012, 18: 9955-9964)\n on top of the classical RRHO.\ntemperature : array-like, optional\n Absolute temperature in Kelvin.\npressure : array-like, optional\n Reference gas pressure.
\n\nReturns
\n\narray-like
\n\nExamples
\n\n\n
>>> import overreact as rx\n>>> from overreact import _constants as constants\n>>> model = rx.parse_model("data/ethane/B97-3c/model.k")\n>>> freeenergies = get_freeenergies(model.compounds, qrrho=(False, True))\n>>> (freeenergies - freeenergies.min()) / constants.kcal\narray([0. , 2.62281461])\n>>> freeenergies = get_freeenergies(model.compounds)\n>>> (freeenergies - freeenergies.min()) / constants.kcal\narray([0. , 2.62862818])\n
\n
You can consider all compounds as solvated if you want:
\n\n\n
>>> sol_freeenergies = get_freeenergies(model.compounds, environment="solvent")\n>>> (sol_freeenergies - freeenergies) / constants.kcal\narray([1.89432845, 1.89432845])\n
\n
You can set a simple energy bias, either as a constant or compound-wise:
\n\n\n
>>> get_freeenergies(model.compounds, bias=1.0) - freeenergies\narray([1., 1.])\n>>> get_freeenergies(model.compounds, bias=-1.0) - freeenergies\narray([-1., -1.])\n>>> get_freeenergies(model.compounds, bias=[1.0, -1.0]) - freeenergies\narray([ 1., -1.])\n
\n
Obtain internal energies for compounds at a given temperature.
\n\nParameters
\n\ncompounds : dict-like\n A descriptor of the compounds.\n Mostly likely, this comes from a parsed input file.\n See overreact.io.parse_model.\nqrrho : bool, optional\n Apply the quasi-rigid rotor harmonic oscillator (QRRHO) approximation of\n M. Head-Gordon and others (see\n J. Phys. Chem. C 2015, 119, 4, 1840-1850)\n on top of the classical RRHO.\ntemperature : array-like, optional\n Absolute temperature in Kelvin.
\n\nReturns
\n\narray-like
\n\nExamples
\n\n\n
>>> import overreact as rx\n>>> from overreact import _constants as constants\n>>> model = rx.parse_model("data/ethane/B97-3c/model.k")\n>>> internal_energies = get_internal_energies(model.compounds)\n>>> (internal_energies - internal_energies.min()) / constants.kcal\narray([0. , 2.20053981])\n
\n
Obtain reaction rate constants for a given reaction scheme.
\n\nParameters
\n\nscheme : Scheme\n A descriptor of the reaction scheme.\n Mostly likely, this comes from a parsed input file.\n See overreact.io.parse_model.\ncompounds : dict-like, optional\n A descriptor of the compounds.\n Mostly likely, this comes from a parsed input file.\n See overreact.io.parse_model.\nbias : array-like, optional\n Energy to be added to free energies.\ntunneling : str or None, optional\n Choose between \"eckart\", \"wigner\" or None (or \"none\").\nqrrho : bool or tuple-like, optional\n Apply both the quasi-rigid rotor harmonic oscillator (QRRHO)\n approximations of M. Head-Gordonand others (enthalpy correction, see\n J. Phys. Chem. C 2015, 119, 4, 1840-1850)\n and S. Grimme (entropy correction, see\n Theory. Chem. Eur. J., 2012, 18: 9955-9964)\n on top of the classical RRHO.\nscale : str, optional\n Reaction rate units. Possible values are \"cm3 mol-1 s-1\",\n \"l mol-1 s-1\", \"m3 mol-1 s-1\", \"cm3 particle-1 s-1\", \"mmHg-1 s-1\" and\n \"atm-1 s-1\".\ntemperature : array-like, optional\n Absolute temperature in Kelvin.\npressure : array-like, optional\n Reference gas pressure.\ndelta_freeenergies : array-like, optional\n Use this instead of obtaining delta free energies from the compounds.\nmolecularity : array-like, optional\n Reaction order, i.e., number of molecules that come together to react.\n If set, this is used to calculate delta_moles for\n overreact.thermo.equilibrium_constant, which effectively calculates a solution\n equilibrium constant between reactants and the transition state for\n gas phase data. You should set this to None if your free energies\n were already adjusted for solution Gibbs free energies.\nvolume : float, optional\n Molar volume.
\n\nReturns
\n\narray-like
\n\nNotes
\n\nSome symbols are accepted as alternatives in \"scale\": \"M-1\", \"ml\" and\n\"torr-1\" are understood as \"l mol-1\", \"cm3\" and \"mmHg-1\", respectively.
\n\nExamples
\n\nBelow is an example of an estimate for the rate of methyl rotation in\nethane (a trivial attempt to reproduce\nScience, 2006, 313, 5795, 1951-1955).\nHow many turns it does per second?
\n\n\n
>>> model = rx.parse_model("data/ethane/B97-3c/model.k")\n>>> get_k(model.scheme, model.compounds)\narray([8.16e+10])\n>>> get_k(model.scheme, model.compounds, qrrho=(False, True))\narray([8.24968117e+10])\n>>> get_k(model.scheme, model.compounds, qrrho=False)\narray([8.26909266e+10])\n>>> get_k(model.scheme, model.compounds, tunneling="wigner")\narray([7.99e+10])\n>>> get_k(model.scheme, model.compounds, tunneling=None)\narray([7.35e+10])\n
\n
The calculated value is off by less than 2% from the experimental value\n(\\( \\frac{1}{12 \\times 10^{-12}} \\text{s}^{-1} = 8.33 \\times 10^{10} \\text{s}^{-1} \\)).\nWe use Eckart tunneling by default, but see the effect of changing it\nabove.
\n\nThe units of the returned reaction rate constants can be selected for\nnon-unimolecular processes. The following is an attempt to reproduce\nJ Atmos Chem, 1996 23, 37-49 for\nthe reaction of proton-withdrawal by a chloride radical from the methane\nmolecule\n\\( \\ce{CH4 + \\cdot Cl -> [H3C\\cdots H\\cdots Cl]^\\ddagger -> H3C\\cdot + HCl} \\):
\n\n\n
>>> model = rx.parse_model("data/tanaka1996/UMP2/cc-pVTZ/model.jk")\n>>> get_k(model.scheme, model.compounds, temperature=300,\n... scale="cm3 particle-1 s-1")\narray([9.60e-14])\n
\n
(By the way, according to the Jet Propulsion Laboratory,\nPublication No. 19-5,\nthe experimental reaction rate constant for this reaction is\n\\( 1.0 \\times 10^{-13} \\text{cm}^3 \\text{particle}^{-1} \\text{s}^{-1} \\).)
\n\nThe returned units are \"M-1 s-1\" by default:
\n\n\n
>>> get_k(model.scheme, model.compounds) \\\n... == get_k(model.scheme, model.compounds, scale="l mol-1 s-1")\narray([ True])\n
\n
You can also turn the tunneling correction off by using the string \"none\":
\n\n\n
>>> get_k(model.scheme, model.compounds, tunneling="none") \\\n... == get_k(model.scheme, model.compounds, tunneling=None)\narray([ True])\n
\n
You can set a simple energy bias, either as a constant or compound-wise:
\n\n\n
>>> get_k(model.scheme, model.compounds, bias=1.0 * constants.kcal,\n... temperature=300.0, scale="cm3 particle-1 s-1")\narray([5.14e-13])\n>>> get_k(model.scheme, model.compounds,\n... bias=np.array([0.0, 0.0, -1.4, 0.0, 0.0]) * constants.kcal,\n... temperature=300.0, scale="cm3 particle-1 s-1")\narray([1.1e-12])\n
\n
Obtain tunneling transmission coefficients at a given temperature.
\n\nOne tunneling transmission coefficient is calculated for each reaction. If\na reaction lacks a transition state (i.e., a half-equilibrium reaction),\nits transmission coefficient is set to unity.
\n\nParameters
\n\nscheme : Scheme\n A descriptor of the reaction scheme.\n Mostly likely, this comes from a parsed input file.\n See overreact.io.parse_model.\ncompounds : dict-like\n A descriptor of the compounds.\n Mostly likely, this comes from a parsed input file.\n See overreact.io.parse_model.\nmethod : str or None, optional\n Choose between \"eckart\", \"wigner\" or None (or \"none\").\nqrrho : bool, optional\n Apply both the quasi-rigid rotor harmonic oscillator (QRRHO)\n approximations of M. Head-Gordon and others (enthalpy correction, see\n J. Phys. Chem. C 2015, 119, 4, 1840-1850)\n and S. Grimme (entropy correction, see\n Theory. Chem. Eur. J., 2012, 18: 9955-9964)\n on top of the classical RRHO.\ntemperature : array-like, optional\n Absolute temperature in Kelvin.
\n\nReturns
\n\narray-like
\n\nRaises
\n\nValueError\n If method is not supported.
\n\nExamples
\n\nBelow is an example of an estimate of how much quantum tunneling\ncontributes to the rate of methyl rotation in ethane (see\nScience, 2006, 313, 5795, 1951-1955\nfor some interesting experimental data on this reaction).
\n\n\n
>>> model = rx.parse_model("data/ethane/B97-3c/model.k")\n>>> kappa = get_kappa(model.scheme, model.compounds)\n>>> kappa\narray([1.110])\n>>> get_kappa(model.scheme, model.compounds, method="none")\narray([1.0])\n>>> get_kappa(model.scheme, model.compounds, method="none") \\\n... == get_kappa(model.scheme, model.compounds, method=None)\narray([ True])\n
\n
You can calculate each piece of the reaction rate constant by hand,\nif you want. Just make sure that you don't calculate the tunneling\ncoefficient twice:
\n\n\n
>>> kappa * get_k(model.scheme, model.compounds, tunneling=None)\narray([8.e+10])\n
\n
Calculate entropy contributions from the overall reaction structure.
\n\nThis function currently implements the reaction translational entropy, a\nresult of the indistinguishability of reactants or products, i.e., a\ndifference in entropy of \\( R ln 2! \\) for the reactions
\n\n$$\\ce{A + B -> C}$$
\n\nand
\n\n$$\\ce{2A -> C}$$
\n\nParameters
\n\ntransform : array-like\ntemperature : array-like, optional\n Absolute temperature in Kelvin.\npressure : array-like, optional\n Reference gas pressure.
\n\nReturns
\n\ndelta_entropy : array-like
\n\nExamples
\n\n\n
>>> import overreact as rx\n
\n
\n
>>> scheme = rx.parse_reactions("A + B <=> C")\n>>> get_reaction_entropies(scheme.A)\narray([0.0, 0.0])\n>>> scheme = rx.parse_reactions("2A <=> C")\n>>> get_reaction_entropies(scheme.A)\narray([-5.763, 5.763])\n>>> get_reaction_entropies(scheme.A, temperature=400.0)\narray([-5.763, 5.763])\n
\n
\n
>>> scheme = rx.parse_reactions("A <=> B + C")\n>>> get_reaction_entropies(scheme.A)\narray([0.0, 0.0])\n>>> scheme = rx.parse_reactions("A <=> 2B")\n>>> get_reaction_entropies(scheme.A)\narray([ 5.763, -5.763])\n
\n
\n
>>> scheme = rx.parse_reactions("A + B -> E# -> C")\n>>> get_reaction_entropies(scheme.A)\narray([0.0])\n>>> get_reaction_entropies(scheme.B)\narray([0.0])\n>>> scheme = rx.parse_reactions("2A -> E# -> C")\n>>> get_reaction_entropies(scheme.A)\narray([-5.763])\n>>> get_reaction_entropies(scheme.B)\narray([-5.763])\n
\n
Return the indices of transition states for each reaction.
\n\nParameters
\n\nA, B : array-like\nis_half_equilibrium : sequence
\n\nReturns
\n\nsequence
\n\nExamples
\n\n\n
>>> scheme = parse_reactions("A -> B")\n>>> print(scheme)\nScheme(compounds=('A', 'B'),\n reactions=('A -> B',),\n is_half_equilibrium=(False,),\n A=((-1.,), (1.,)),\n B=((-1.,), (1.,)))\n>>> get_transition_states(scheme.A, scheme.B, scheme.is_half_equilibrium)\n(None,)\n
\n
\n
>>> scheme = parse_reactions("S -> E\u2021 -> S")\n>>> print(scheme)\nScheme(compounds=('S', 'E\u2021'),\n reactions=('S -> S',),\n is_half_equilibrium=(False,),\n A=((0.,), (0.,)),\n B=((-1.,), (1.,)))\n>>> get_transition_states(scheme.A, scheme.B, scheme.is_half_equilibrium)\n(1,)\n
\n
\n
>>> scheme = parse_reactions("E + S <=> ES -> ES\u2021 -> E + P")\n>>> print(scheme)\nScheme(compounds=('E', 'S', 'ES', 'ES\u2021', 'P'),\n reactions=('E + S -> ES', 'ES -> E + S', 'ES -> E + P'),\n is_half_equilibrium=(True, True, False),\n A=((-1., 1., 1.),\n (-1., 1., 0.),\n (1., -1., -1.),\n (0., 0., 0.),\n (0., 0., 1.)),\n B=((-1., 0., 0.),\n (-1., 0., 0.),\n (1., 0., -1.),\n (0., 0., 1.),\n (0., 0., 0.)))\n>>> get_transition_states(scheme.A, scheme.B, scheme.is_half_equilibrium)\n(None, None, 3)\n
\n
Simulate a reaction scheme from its rate function.
\n\nThis function provides two functions that calculate the concentrations and\nthe rates of formation at any point in time for any compound. It does that\nby solving an initial value problem (IVP) through scipy's solve_ivp\nunder the hood.
\n\nParameters
\n\ndydt : callable\n Right-hand side of the system.\ny0 : array-like\n Initial state.\nt_span : array-like, optional\n Interval of integration (t0, tf). The solver starts with t=t0 and\n integrates until it reaches t=tf. If not given, a conservative value\n is chosen based on the system at hand (the method of choice works for\n any zeroth-, first- or second-order reactions).\nmethod : str, optional\n Integration method to use. See scipy.integrate.solve_ivp for details.\n Kinetics problems are very often stiff and, as such, \"RK23\" and \"RK45\" may be\n unsuited. \"LSODA\", \"BDF\", and \"Radau\" are worth a try if things go bad.\nmax_step : float, optional\n Maximum step to be performed by the integrator.\n Defaults to half the total time span.\nfirst_step : float, optional\n First step size.\n Defaults to half the maximum step, or np.finfo(np.float64).eps,\n whichever is smallest.\nrtol, atol : array-like, optional\n See scipy.integrate.solve_ivp for details.\nmax_time : float, optional\n If t_span is not given, an interval will be estimated, but it can't\n be larger than this parameter.
\n\nReturns
\n\ny, r : callable\n Concentrations and reaction rates as functions of time. The y object\n is an OdeSolution and stores attributes t_min and t_max.
\n\nExamples
\n\n\n
>>> import numpy as np\n>>> import overreact as rx\n
\n
A toy simulation can be performed in just two lines:
\n\n\n
>>> scheme = rx.parse_reactions("A <=> B")\n>>> y, r = get_y(get_dydt(scheme, np.array([1, 1])), y0=[1, 0])\n
\n
The y object stores information about the simulation time, which can be\nused to produce a suitable vector of timepoints for, e.g., plotting:
\n\n\n
>>> y.t_min, y.t_max\n(0.0, 3.0)\n>>> t = np.linspace(y.t_min, y.t_max)\n>>> t\narray([0. , 0.06122449, ..., 2.93877551, 3. ])\n
\n
Both y and r can be used to check concentrations and rates in any\npoint in time. In particular, both are vectorized:
\n\n\n
>>> y(t)\narray([[1. , ...],\n [0. , ...]])\n>>> r(t)\narray([[-1. , ...],\n [ 1. , ...]])\n
\n
Check whether a name specifies a transition state.
\n\nParameters
\n\nname : str
\n\nReturns
\n\nbool
\n\nExamples
\n\n\n
>>> is_transition_state("A#")\nTrue\n>>> is_transition_state("A\u2021")\nTrue\n>>> is_transition_state("pyrrole#")\nTrue\n>>> is_transition_state("pyrrole\u2021")\nTrue\n>>> is_transition_state("A")\nFalse\n>>> is_transition_state("A~")\nFalse\n>>> is_transition_state("pyrrole")\nFalse\n>>> is_transition_state("pyrrole~")\nFalse\n
\n
This function also works for names that specify environment:
\n\n\n
>>> is_transition_state("A#(w)")\nTrue\n>>> is_transition_state("A\u2021(w)")\nTrue\n>>> is_transition_state("TS#(w)")\nTrue\n>>> is_transition_state("TS\u2021(w)")\nTrue\n>>> is_transition_state("A(w)")\nFalse\n>>> is_transition_state("A~(w)")\nFalse\n>>> is_transition_state("TS(w)")\nFalse\n>>> is_transition_state("TS~(w)")\nFalse\n
\n
Parse a set of compounds.
\n\nParameters
\n\ntext : str, sequence of str or dict-like\n Compound descriptions or sequence of lines of it.\npath : sequence of str, optional\n Paths for searching logfiles.\nselect : sequence of str, optional\n If defined, only those compounds will be returned.
\n\nReturns
\n\ncompounds : immutable dict-like
\n\nRaises
\n\nFileNotFoundError\n If a logfile is not found.
\n\nExamples
\n\n\n
>>> import overreact as rx\n
\n
\n
>>> compounds = rx.parse_compounds("S: data/ethane/B97-3c/staggered.out")\n>>> compounds\n{'S': {'logfile': 'data/ethane/B97-3c/staggered.out',\n 'energy': -209483812.77142256,\n 'mult': 1,\n 'atomnos': (6, 6, 1, 1, 1, 1, 1, 1),\n ...}}\n
\n
\n
>>> compounds = rx.parse_compounds({"S": "data/ethane/B97-3c/staggered.out"})\n>>> compounds\n{'S': {'logfile': 'data/ethane/B97-3c/staggered.out',\n 'energy': -209483812.77142256,\n 'mult': 1,\n 'atomnos': (6, 6, 1, 1, 1, 1, 1, 1),\n ...}}\n
\n
\n
>>> compounds = rx.parse_compounds(["S: data/ethane/B97-3c/staggered.out",\n... "E\u2021: data/ethane/B97-3c/eclipsed.out"])\n>>> compounds\n{'S': {'logfile': 'data/ethane/B97-3c/staggered.out',\n 'energy': -209483812.77142256,\n 'mult': 1,\n 'atomnos': (6, 6, 1, 1, 1, 1, 1, 1),\n ...},\n'E\u2021': {'logfile': 'data/ethane/B97-3c/eclipsed.out',\n 'energy': -209472585.3539883,\n 'mult': 1,\n 'atomnos': (6, 6, 1, 1, 1, 1, 1, 1),\n ...}}\n
\n
\n
>>> compounds = rx.parse_compounds('''S: data/ethane/B97-3c/staggered.out\n... E\u2021: data/ethane/B97-3c/eclipsed.out''')\n>>> compounds\n{'S': {'logfile': 'data/ethane/B97-3c/staggered.out',\n 'energy': -209483812.77142256,\n 'mult': 1,\n 'atomnos': (6, 6, 1, 1, 1, 1, 1, 1),\n ...},\n'E\u2021': {'logfile': 'data/ethane/B97-3c/eclipsed.out',\n 'energy': -209472585.3539883,\n 'mult': 1,\n 'atomnos': (6, 6, 1, 1, 1, 1, 1, 1),\n ...}}\n
\n
\n
>>> compounds = rx.parse_compounds({"S": "data/ethane/B97-3c/staggered.out",\n... "E\u2021": "data/ethane/B97-3c/eclipsed.out"})\n>>> compounds\n{'S': {'logfile': 'data/ethane/B97-3c/staggered.out',\n 'energy': -209483812.77142256,\n 'mult': 1,\n 'atomnos': (6, 6, 1, 1, 1, 1, 1, 1),\n ...},\n'E\u2021': {'logfile': 'data/ethane/B97-3c/eclipsed.out',\n 'energy': -209472585.3539883,\n 'mult': 1,\n 'atomnos': (6, 6, 1, 1, 1, 1, 1, 1),\n ...}}\n
\n
Parse either a source or model input file, whichever is available.
\n\nA source input file (also known as a .k file) contains all the information needed\nto create a model input file.\nA model input file (also known as a .jk file) is a JSON encoded file with all the\ninformation needed to study microkinetic simulations from first principles.
\n\nYou probably won't need to use model input files directly, they are\nautomatically created based on source input files.\nTake a look at our guide on how to write an source input file.
\n\nThis function attempts to parse a model input file if available. If not, a source\ninput file is parsed and a model input file is generated from it. Extensions are\nguessed if none given (i.e., if only the base name given).
\n\nParameters
\n\npath : str\n Path to the model or source input file.\n If the final extension is not .jk or .k, it is guessed.\nforce_compile : bool\n If True, a .k file will take precedence over any .jk file for reading. A\n .jk file is thus either generated or overwritten. This is sometimes\n needed to force an update with new data.
\n\nReturns
\n\nmodel : immutable dict-like
\n\nRaises
\n\nFileNotFoundError\n If the model or source input file is not found.
\n\nExamples
\n\nSome examples of how overreact \"sees\" your data below \ud83d\ude04:
\n\n\n
>>> model = parse_model("data/ethane/B97-3c/model.jk")\n>>> model.scheme\nScheme(compounds=('S', 'E\u2021'),\n reactions=('S -> S',),\n is_half_equilibrium=(False,),\n A=((0.0,), (0.0,)),\n B=((-1.0,), (1.0,)))\n>>> model.compounds["S"]\n{'logfile': 'data/ethane/B97-3c/staggered.out',\n 'energy': -209483812.77142256,\n 'mult': 1,\n 'atomnos': (6, 6, 1, 1, 1, 1, 1, 1),\n 'atommasses': (12.011, 12.011, 1.008, 1.008, 1.008, 1.008, 1.008, 1.008),\n 'atomcoords': ((-7.633588, 2.520693, -4.8e-05),\n ...,\n (-5.832852, 3.674431, 0.363239)),\n 'vibfreqs': (307.57, 825.42, ..., 3071.11, 3071.45),\n 'vibdisps': (((-1.7e-05, 3.4e-05, 5.4e-05),\n ...,\n (-0.011061, -0.030431, -0.027036)))}\n>>> model_from_source = parse_model("data/ethane/B97-3c/model.k",\n... force_compile=True)\n>>> model_from_source == model\nTrue\n>>> model_from_source = parse_model("data/ethane/B97-3c/model")\n>>> model_from_source == model\nTrue\n
\n
Parse a kinetic model as a chemical reaction scheme.
\n\nThis is an essential part of the parsing process.\nSee overreact.io.parse_model other details.
\n\nParameters
\n\ntext : str or sequence of str\n Model description or sequence of lines of it.
\n\nReturns
\n\nscheme : Scheme\n A descriptor of the reaction scheme.
\n\nNotes
\n\nThe model description should comply with the mini-language for systems of\nreactions. A semi-formal definition of the grammar in\nBackus-Naur form\nis given below:
\n\n equation ::= equation_side arrow equation_side\nequation_side ::= coefficient compound ['+' coefficient compound]*\n coefficient ::= [integers] (defaults to 1)\n compound ::= mix of printable characters\n arrow ::= '->' | '<=>' | '<-'\n
\n\nBlank lines and comments (starting with //) are ignored. Repeated\nreactions are ignored. Furthermore, reactions can be chained one after\nanother and, if a single compound (with either a \u2021 or a # at the end)\nappears alone on one side of a reaction, it's considered a transition\nstate. Transition states have zero lifetime during the simulation.
\n\nExamples
\n\nWhat follows is a rather long tour over the parsing process and its\noutput in general. You can skip it if you are not interested in the\ndetails.
\n\n\n
>>> scheme = parse_reactions("A -> B // a direct reaction")\n
\n
The reaction above is a direct one (observe that comments are ignored). The\nreturned object has the following attributes:
\n\n\n
>>> scheme.compounds\n('A', 'B')\n>>> scheme.reactions\n('A -> B',)\n>>> scheme.is_half_equilibrium\n(False,)\n>>> scheme.A\n((-1.,), (1.,))\n>>> scheme.B\n((-1.,), (1.,))\n
\n
The same reaction can be specified in reverse order:
\n\n\n
>>> parse_reactions("B <- A // reverse reaction of the above")\nScheme(compounds=('A', 'B'),\n reactions=('A -> B',),\n is_half_equilibrium=(False,),\n A=((-1.,), (1.,)),\n B=((-1.,), (1.,)))\n
\n
Equilibria produce twice as many direct reactions, while the $B$ matrix\ndefines an energy relationship for only one of each pair:
\n\n\n
>>> parse_reactions("A <=> B // an equilibrium")\nScheme(compounds=('A', 'B'),\n reactions=('A -> B', 'B -> A'),\n is_half_equilibrium=(True, True),\n A=((-1., 1.),\n (1., -1.)),\n B=((-1., 0.),\n (1., 0.)))\n
\n
Adding twice the same reaction results in a single reaction being added.\nThis of course also works with equilibria (extra whitespaces are ignored):
\n\n\n
>>> parse_reactions('''\n... A <=> B -> A\n... A -> B <=> A\n... A -> B <- A\n... B <- A -> B\n... ''')\nScheme(compounds=('A', 'B'),\n reactions=('A -> B', 'B -> A'),\n is_half_equilibrium=(True, True),\n A=((-1., 1.),\n (1., -1.)),\n B=((-1., 0.),\n (1., 0.)))\n
\n
Transition states are specified with a special symbol at the end (either\n\u2021 or #). They are shown among compounds, but the matrix $A$ ensures\nthey'll never have a non-zero rate of formation/consumption. On the other\nhand, they are needed in the $B$ matrix:
\n\n\n
>>> parse_reactions("A -> A\u2021 -> B")\nScheme(compounds=('A', 'A\u2021', 'B'),\n reactions=('A -> B',),\n is_half_equilibrium=(False,),\n A=((-1.,), (0.,), (1.,)),\n B=((-1.,), (1.,), (0.,)))\n
\n
This gives the same result as above:
\n\n\n
>>> parse_reactions("A -> A\u2021 -> B <- A\u2021 <- A")\nScheme(compounds=('A', 'A\u2021', 'B'),\n reactions=('A -> B',),\n is_half_equilibrium=(False,),\n A=((-1.,), (0.,), (1.,)),\n B=((-1.,), (1.,), (0.,)))\n
\n
It is possible to define a reaction whose product is the same as the\nreactant. This is found in isomerization processes (e.g., ammonia\ninversion or the methyl rotation in ethane):
\n\n\n
>>> parse_reactions("S -> E\u2021 -> S")\nScheme(compounds=('S', 'E\u2021'),\n reactions=('S -> S',),\n is_half_equilibrium=(False,),\n A=((0.,), (0.,)),\n B=((-1.,), (1.,)))\n
\n
As such, a column full of zeros in the $A$ matrix corresponds to a reaction\nwith zero net change. As can be seen, overreact allows for very general\nmodels. An interesting feature is that a single transition state can link\nmany different compounds (whether it is useful is a matter of debate):
\n\n\n
>>> parse_reactions('''\n... B -> B\u2021 -> C // chained reactions and transition states\n... B\u2021 -> D // this is a bifurcation\n... B -> B'\u2021 -> E // this is a classical competitive reaction\n... A -> B\u2021\n... ''')\nScheme(compounds=('B', 'B\u2021', 'C', 'D', "B'\u2021", 'E', 'A'),\n reactions=('B -> C', 'B -> D', 'B -> E', 'A -> C', 'A -> D'),\n is_half_equilibrium=(False, False, False, False, False),\n A=((-1., -1., -1., 0., 0.),\n (0., 0., 0., 0., 0.),\n (1., 0., 0., 1., 0.),\n (0., 1., 0., 0., 1.),\n (0., 0., 0., 0., 0.),\n (0., 0., 1., 0., 0.),\n (0., 0., 0., -1., -1.)),\n B=((-1., -1., -1., 0., 0.),\n (1., 1., 0., 1., 1.),\n (0., 0., 0., 0., 0.),\n (0., 0., 0., 0., 0.),\n (0., 0., 1., 0., 0.),\n (0., 0., 0., 0., 0.),\n (0., 0., 0., -1., -1.)))\n
\n
The following is a borderline case but both reactions should be considered\ndifferent since they define different processes:
\n\n\n
>>> parse_reactions('''\n... A -> A\u2021 -> B\n... A -> B\n... ''')\nScheme(compounds=('A', 'A\u2021', 'B'),\n reactions=('A -> B', 'A -> B'),\n is_half_equilibrium=(False, False),\n A=((-1., -1.),\n (0., 0.),\n (1., 1.)),\n B=((-1., -1.),\n (1., 0.),\n (0., 1.)))\n
\n
The following is correct behavior. In fact, the reactions are badly\ndefined: if more than one transition state are chained, the following\nhappens, which is correct since it's the most physically plausible model\nthat can be extracted. It can be seen as a feature that the product B is\nignored and not the reactant A, since the user would easily see the mistake\nin graphs of concentration over time (the alternative would be no\nreaction happening at all, which is rather cryptic to debug).
\n\n\n
>>> parse_reactions("A -> A\u2021 -> A'\u2021 -> B")\nScheme(compounds=('A', 'A\u2021', "A'\u2021", 'B'),\n reactions=("A -> A'\u2021",),\n is_half_equilibrium=(False,),\n A=((-1.,), (0.,), (1.,), (0.,)),\n B=((-1.,), (1.,), (0.,), (0.,)))\n
\n
In any case, it's not clear how a reaction barrier be defined in such a\ncase. If you have a use case, don't hesitate to\nopen an issue, we'll be\nhappy to hear from you.
\n", "signature": "(text: Union[str, Sequence[str]]) -> overreact.core.Scheme:", "funcdef": "def"}, "overreact.unparse_reactions": {"fullname": "overreact.unparse_reactions", "modulename": "overreact", "qualname": "unparse_reactions", "kind": "function", "doc": "Unparse a kinetic model.
\n\nParameters
\n\nscheme : Scheme\n A descriptor of the reaction scheme.\n Mostly likely, this comes from a parsed input file.\n See overreact.io.parse_model.
\n\nReturns
\n\ntext : str
\n\nNotes
\n\nThis function assumes complimentary half equilibria are located one after\nthe other in scheme.reactions, which is to be expected from\nparse_reactions.
\n\nExamples
\n\n\n
>>> unparse_reactions(Scheme(compounds=('A', 'B'), reactions=('A -> B',),\n... is_half_equilibrium=(False,),\n... A=((-1.,),\n... ( 1.,)),\n... B=((-1.,),\n... ( 1.,))))\n'A -> B'\n>>> unparse_reactions(Scheme(compounds=('A', 'B'),\n... reactions=('A -> B', 'B -> A'),\n... is_half_equilibrium=(True, True),\n... A=((-1., 1.),\n... ( 1., -1.)),\n... B=((-1., 0.),\n... ( 1., 0.))))\n'A <=> B'\n>>> unparse_reactions(Scheme(compounds=('A', 'A\u2021', 'B'),\n... reactions=('A -> B',),\n... is_half_equilibrium=(False,),\n... A=((-1.,),\n... ( 0.,),\n... ( 1.,)),\n... B=((-1.,),\n... ( 1.,),\n... ( 0.,))))\n'A -> A\u2021 -> B'\n>>> print(unparse_reactions(Scheme(compounds=('B', 'B\u2021', 'C', 'D', "B'\u2021",\n... 'E', 'A'),\n... reactions=('B -> C', 'B -> D', 'B -> E',\n... 'A -> C', 'A -> D'),\n... is_half_equilibrium=(False, False, False,\n... False, False),\n... A=((-1., -1., -1., 0., 0.),\n... ( 0., 0., 0., 0., 0.),\n... ( 1., 0., 0., 1., 0.),\n... ( 0., 1., 0., 0., 1.),\n... ( 0., 0., 0., 0., 0.),\n... ( 0., 0., 1., 0., 0.),\n... ( 0., 0., 0., -1., -1.)),\n... B=((-1., -1., -1., 0., 0.),\n... ( 1., 1., 0., 1., 1.),\n... ( 0., 0., 0., 0., 0.),\n... ( 0., 0., 0., 0., 0.),\n... ( 0., 0., 1., 0., 0.),\n... ( 0., 0., 0., 0., 0.),\n... ( 0., 0., 0., -1., -1.)))))\nB -> B\u2021 -> C\nB -> B\u2021 -> D\nB -> B'\u2021 -> E\nA -> B\u2021 -> C\nA -> B\u2021 -> D\n>>> print(unparse_reactions(Scheme(compounds=('A', 'A\u2021', 'B'),\n... reactions=('A -> B', 'A -> B'),\n... is_half_equilibrium=(False, False),\n... A=((-1., -1.),\n... ( 0., 0.),\n... ( 1., 1.)),\n... B=((-1., -1.),\n... ( 1., 0.),\n... ( 0., 1.)))))\nA -> A\u2021 -> B\nA -> B\n>>> unparse_reactions(Scheme(compounds=('A', 'A\u2021', "A'\u2021", 'B'),\n... reactions=("A -> A'\u2021",),\n... is_half_equilibrium=(False,),\n... A=((-1.,),\n... ( 0.,),\n... ( 1.,),\n... ( 0.,)),\n... B=((-1.,),\n... ( 1.,),\n... ( 0.,),\n... ( 0.,))))\n"A -> A\u2021 -> A'\u2021"\n>>> unparse_reactions(Scheme(compounds=('S', 'E\u2021'),\n... reactions=('S -> S',),\n... is_half_equilibrium=(False,),\n... A=((0.0,),\n... (0.0,)),\n... B=((-1.0,),\n... (1.0,))))\n'S -> E\u2021 -> S'\n
\n
---
\n\n\n\n\n
\n
---
\n\noverreact is a library and a command-line tool for building and\nanalyzing homogeneous microkinetic models from first-principles\ncalculations:
\n\n\n
In [1]: from overreact import api # the api\n\nIn [2]: api.get_k("S -> E\u2021 -> S", # your model\n ...: {"S": "data/ethane/B97-3c/staggered.out", # your data\n ...: "E\u2021": "data/ethane/B97-3c/eclipsed.out"})\nOut[2]: array([8.16880917e+10]) # your results\n
\n
The user specifies a set of\nelementary reactions that are believed to be relevant for the overall chemical\nphenomena. overreact offers a hopefully complete but simple environment for\nhypothesis testing in first-principles chemical kinetics.
\n\n\n \n \ud83e\udd14 What is microkinetic modeling?\n
\n \n Microkinetic modeling is a technique used to predict the outcome\n of complex chemical reactions.\n It can be used\n to investigate the catalytic transformations\n of molecules.\n overreact makes it easy to create\n and analyze microkinetic models built\n from computational chemistry data.\n
\n
\n\n
\n\n\n \n \ud83e\uddd0 What do you mean by first-principles calculations?\n
\n \n We use the term first-principles calculations to refer to\n calculations performed using quantum chemical modern methods such as\n Wavefunction\n and\n Density Functional\n theories.\n For instance, the three-line example code above calculates the rate of methyl rotation in ethane (at\n B97-3c).\n (Rather surprisingly, the error found is less than 2%\n when compared to available experimental results.)\n
\n
\n\n
\n\noverreact uses precise thermochemical partition funtions, tunneling\ncorrections and data is parsed directly from computational chemistry\noutput files thanks to cclib (see the\nlist of its supported programs).
\n\nInstallation
\n\noverreact is a Python package, so you can easily install it with\npip:
\n\n\n
$ pip install "overreact[cli,fast]"\n
\n
See the\ninstallation guide\nfor more details.
\n\n\n \ud83d\ude80 Where to go from here? Take a look at the\n short introduction.\n Or see\n below\n for more guidance.
\n
\n\nCiting overreact
\n\nIf you use overreact in your research, please cite:
\n\n\n Schneider, F. S. S.; Caramori, G. F.\n _Overreact, an in Silico Lab: Automative Quantum Chemical Microkinetic Simulations for Complex Chemical Reactions_.\n Journal of Computational Chemistry 2022, 44 (3), 209\u2013217.\n doi:10.1002/jcc.26861.
\n
\n\nHere's the reference in BibTeX format:
\n\n\n
@article{overreact_paper2022,\n title = {Overreact, an in silico lab: Automative quantum chemical microkinetic simulations for complex chemical reactions},\n author = {Schneider, Felipe S. S. and Caramori, Giovanni F.},\n year = {2022},\n month = {Apr},\n journal = {Journal of Computational Chemistry},\n publisher = {Wiley},\n volume = {44},\n number = {3},\n pages = {209\u2013217},\n doi = {10.1002/jcc.26861},\n issn = {1096-987x},\n url = {http://dx.doi.org/10.1002/jcc.26861},\n}\n@software{overreact_software2021,\n title = {geem-lab/overreact: v1.2.0 \\vert{} Zenodo},\n author = {Felipe S. S. Schneider and Let\\'{\\i}cia M. P. Madureira and Giovanni F. Caramori},\n year = {2023},\n month = {Jan},\n publisher = {Zenodo},\n doi = {10.5281/zenodo.7865357},\n url = {https://doi.org/10.5281/zenodo.7865357},\n version = {v1.2.0},\n howpublished = {\\url{https://github.com/geem-lab/overreact}},\n}\n
\n
License
\n\noverreact is open-source, released under the permissive MIT license. See\nthe LICENSE agreement.
\n\nFunding
\n\nThis project was developed at the GEEM lab\n(Federal University of Santa Catarina, Brazil), and was\npartially funded by the\nBrazilian National Council for Scientific and Technological Development (CNPq),\ngrant number 140485/2017-1.
\n"}, "overreact.Scheme": {"fullname": "overreact.Scheme", "modulename": "overreact", "qualname": "Scheme", "kind": "class", "doc": "A descriptor of a chemical reaction network.
\n\nMostly likely, this comes from a parsed input file.\nSee overreact.io.parse_model.
\n", "bases": "typing.NamedTuple"}, "overreact.Scheme.__init__": {"fullname": "overreact.Scheme.__init__", "modulename": "overreact", "qualname": "Scheme.__init__", "kind": "function", "doc": "Create new instance of Scheme(compounds, reactions, is_half_equilibrium, A, B)
\n", "signature": "(\tcompounds: Sequence[str],\treactions: Sequence[str],\tis_half_equilibrium: Sequence[bool],\tA: numpy.ndarray,\tB: numpy.ndarray)"}, "overreact.Scheme.compounds": {"fullname": "overreact.Scheme.compounds", "modulename": "overreact", "qualname": "Scheme.compounds", "kind": "variable", "doc": "A descriptor of compounds.
\n", "annotation": ": Sequence[str]"}, "overreact.Scheme.reactions": {"fullname": "overreact.Scheme.reactions", "modulename": "overreact", "qualname": "Scheme.reactions", "kind": "variable", "doc": "A descriptor of reactions.
\n", "annotation": ": Sequence[str]"}, "overreact.Scheme.is_half_equilibrium": {"fullname": "overreact.Scheme.is_half_equilibrium", "modulename": "overreact", "qualname": "Scheme.is_half_equilibrium", "kind": "variable", "doc": "An indicator of whether a reaction is half-equilibrium.
\n", "annotation": ": Sequence[bool]"}, "overreact.Scheme.A": {"fullname": "overreact.Scheme.A", "modulename": "overreact", "qualname": "Scheme.A", "kind": "variable", "doc": "A matrix of stoichiometric coefficients between reactants and products.
\n", "annotation": ": numpy.ndarray"}, "overreact.Scheme.B": {"fullname": "overreact.Scheme.B", "modulename": "overreact", "qualname": "Scheme.B", "kind": "variable", "doc": "A matrix of stoichiometric coefficients between reactants and transition states.
\n", "annotation": ": numpy.ndarray"}, "overreact.change_reference_state": {"fullname": "overreact.change_reference_state", "modulename": "overreact", "qualname": "change_reference_state", "kind": "function", "doc": "Calculate an additive entropy correction to a change in reference states.
\n\n$$\\Delta G_\\text{corr} =\n R T \\ln \\left( \\frac{\\chi_\\text{new}}{\\chi_\\text{old}} \\right)$$
\n\nThe value returned can be directly multiplied by temperature and summed to\nthe old reference free energies to obtain free energies with respect to a\nnew reference. See notes below.
\n\nFor instance, the concentration correction to Gibbs free energy for a\ngas-to-liquid standard state change is simply\n(\\( c^\\circ = \\frac{\\text{1 atm}}{R T} \\)),
\n\n$$\\Delta G_\\text{conc} =\n R T \\ln \\left( \\frac{\\text{1 M}}{c^\\circ} \\right)$$
\n\nParameters
\n\nnew_reference : array-like, optional\n New reference state. Default value corresponds to 1 mol/liter.\nold_reference : array-like, optional\n Old reference state. Default value corresponds to the concentration of\n an ideal gas at the given temperature and 1 atm.\nsign : float, optional\n Sign of the change in reference state. Default value is 1. This only\n multiplies the final result.\ntemperature : array-like, optional\n Absolute temperature in Kelvin.\npressure : array-like, optional\n Reference gas pressure.\nvolume : float, optional\n Molar volume.
\n\nReturns
\n\ncorrection : array-like\n Entropy correction in J/mol\u00b7K.
\n\nNotes
\n\nThis function can be used to add any entropy correction in the form above.\nThe only drawback is that, sometimes, those corrections are written with a\nminus sign in front of them (this implies switching the roles of\nold_reference and new_reference). The easiest way to accomplish this is\nby using sign=-1 or multiplying the result by -1.
\n\nExamples
\n\nBy default, the correction returns a change in concentration from the gas\nphase standard concentration to the solvated-state standard concentration:
\n\n\n
>>> -rx.change_reference_state() / constants.calorie\n-6.4\n>>> 298.15 * rx.change_reference_state() / constants.kcal\n1.89\n>>> 273.15 * rx.change_reference_state(temperature=273.15) / constants.kcal\n1.69\n
\n
But this function can also be used to adjust symmetry effects from C1\ncalculations (symmetry number equals to one). For D7h, for instance, the\nsymmetry number is 14:
\n\n\n
>>> -298.15 * rx.change_reference_state(14, 1) / constants.kcal\n-1.56\n
\n
\n
>>> rx.change_reference_state(sign=-1) == -rx.change_reference_state()\nTrue\n
\n
Estimate a energy bias for a given set of reference data points.
\n\nParameters
\n\nscheme : Scheme\n A descriptor of the reaction scheme.\n Mostly likely, this comes from a parsed model input file.\n See overreact.io.parse_model.\ncompounds : dict-like\n A descriptor of the compounds.\n Mostly likely, this comes from a parsed model input file.\n See overreact.io.parse_model.\ndata : dict-like of array-like\ny0: array-like\ntunneling : str or None, optional\n Choose between \"eckart\", \"wigner\" or None (or \"none\").\nqrrho : bool or tuple-like, optional\n Apply both the quasi-rigid rotor harmonic oscillator (QRRHO)\n approximations of M. Head-Gordon and others (enthalpy correction, see\n J. Phys. Chem. C 2015, 119, 4, 1840-1850) and S. Grimme (entropy correction, see\n Theory. Chem. Eur. J., 2012, 18: 9955-9964) on top of the classical RRHO.\ntemperature : array-like, optional\n Absolute temperature in Kelvin.\npressure : array-like, optional\n Reference gas pressure.\ndelta_freeenergies : array-like, optional\n Use this instead of obtaining delta free energies from the compounds.\nmolecularity : array-like, optional\n Reaction order, i.e., number of molecules that come together to react.\n If set, this is used to calculate delta_moles for\n equilibrium_constant, which effectively calculates a solution\n equilibrium constant between reactants and the transition state for\n gas phase data. You should set this to None if your free energies\n were already adjusted for solution Gibbs free energies.\nvolume : float, optional\n Molar volume.
\n\nReturns
\n\narray-like
\n\nExamples
\n\n\n
>>> model = rx.parse_model("data/tanaka1996/UMP2/cc-pVTZ/model.jk")\n
\n
The following are some estimates on actual atmospheric concentrations:
\n\n\n
>>> y0 = [4.8120675684099e-5,\n... 2.8206357713029e-5,\n... 0.0,\n... 0.0,\n... 2.7426565371219e-5]\n>>> data = {"t": [1.276472128376942246e-6,\n... 1.446535794555581743e-4,\n... 1.717069678525567564e-2],\n... "CH3\u00b7": [9.694916853338366211e-9,\n... 1.066033349343709026e-6,\n... 2.632179124780495175e-5]}\n>>> get_bias(model.scheme, model.compounds, data, y0) / constants.kcal\n-1.4\n
\n
Calculate deltas according to reactions.
\n\nDelta properties are differences in a property between the final and\ninitial state of a chemical transformation. They are calculated from\nmatrices representing the transformation and the absolute properties.\nTransformation matrices are expected to have column-wise transformation\ndefined.
\n\nVery useful for the calculation of reaction and activation free energies\nfrom absolute free energies of compounds. Matrices A and B of a\nScheme represent the transformations associated with reaction and\nactivation free energies, respectively.
\n\nParameters
\n\ntransform : array-like\nproperty : array-like
\n\nReturns
\n\ndelta_property : array-like
\n\nExamples
\n\n\n
>>> get_delta([-1, 1], [-10, 10])\n20\n
\n
You must ensure the transformation is properly defined, as no test is made\nto ensure, e.g., conservation of matter:
\n\n\n
>>> get_delta([-1, 0], [-10, 20])\n10\n
\n
Normally, transformations are given as columns in a matrix:
\n\n\n
>>> get_delta([[-1, -2],\n... [ 1, 3]], [-5, 12])\narray([17, 46])\n
\n
Generate a rate function that models a reaction scheme.
\n\nParameters
\n\nscheme : Scheme\n A descriptor of the reaction scheme.\n Mostly likely, this comes from a parsed model input file.\n See overreact.io.parse_model.\nk : array-like\n Reaction rate constant(s). Units match the concentration units given to\n the returned function dydt.\nef : float, optional\n Equilibrium factor. This is a parameter that can be used to scale the\n reaction rates associated to half-equilibrium reactions such that they\n are faster than the other reactions.
\n\nReturns
\n\ndydt : callable\n Reaction rate function. The actual reaction rate constants employed\n are stored in the attribute k of the returned function. If JAX is\n available, the attribute jac will hold the Jacobian function of\n dydt.
\n\nNotes
\n\nThe returned function is suited to be used by ODE solvers such as\nscipy.integrate.solve_ivp or the older scipy.integrate.ode (see\nexamples below). This is actually what the function get_y from the\ncurrent module does.
\n\nExamples
\n\n\n
>>> import numpy as np\n>>> import overreact as rx\n
\n
\n
>>> scheme = rx.parse_reactions("A <=> B")\n>>> dydt = get_dydt(scheme, np.array([1, 1]))\n>>> dydt(0.0, np.array([1., 1.]))\nArray([0., 0.], ...)\n
\n
If available, JAX is used for JIT compilation. This will make dydt\ncomplain if given lists instead of numpy arrays. So stick to the safer,\nfaster side as above.
\n\nThe actually used reaction rate constants can be inspected with the k\nattribute of dydt:
\n\n\n
>>> dydt.k\nArray([1., 1.], ...)\n
\n
If JAX is available, the Jacobian function will be available as\ndydt.jac:
\n\n\n
>>> dydt.jac(0.0, np.array([1., 1.]))\nArray([[-1., 1.],\n [ 1., -1.]], ...)\n
\n
Obtain enthalpies for compounds at a given temperature.
\n\nParameters
\n\ncompounds : dict-like\n A descriptor of the compounds.\n Mostly likely, this comes from a parsed input file.\n See overreact.io.parse_model.\nqrrho : bool, optional\n Apply the quasi-rigid rotor harmonic oscillator (QRRHO) approximation of\n M. Head-Gordon and others (see\n J. Phys. Chem. C 2015, 119, 4, 1840-1850)\n on top of the classical RRHO.\ntemperature : array-like, optional\n Absolute temperature in Kelvin.
\n\nReturns
\n\narray-like
\n\nExamples
\n\n\n
>>> import overreact as rx\n>>> from overreact import _constants as constants\n>>> model = rx.parse_model("data/ethane/B97-3c/model.k")\n>>> enthalpies = get_enthalpies(model.compounds)\n>>> (enthalpies - enthalpies.min()) / constants.kcal\narray([0. , 2.20053981])\n
\n
The enthalpies at absolute zero can easily be obtained (this is used,\ne.g., in the calculation of the Eckart tunneling coefficient, see\noverreact.tunnel.eckart). We can use this to calculate, for instance,\nthe thermal contributions to the enthalpy:
\n\n\n
>>> zero_enthalpies = get_enthalpies(model.compounds, temperature=0)\n>>> (enthalpies - zero_enthalpies) / constants.kcal\narray([2.78, 2.50])\n
\n
Obtain entropies for compounds at a given temperature and pressure.
\n\nParameters
\n\ncompounds : dict-like\n A descriptor of the compounds.\n Mostly likely, this comes from a parsed input file.\n See overreact.io.parse_model.\nenvironment : str or None, optional\n Choose between \"gas\" and a solvent. This is chosen for you by default,\n based on the names of each compound (e.g. A(g) or A is gas,\n A(w) or A(...) is solvated). In case this is given, all compounds\n will have the same behavior.\nmethod : str, optional\n This is a placeholder for future functionality.\n There are plans to implement more sophisticated methods for calculating\n entropies such as in\n Phys. Chem. Chem. Phys., 2019, 21, 18920-18929\n and\n J. Chem. Theory Comput. 2019, 15, 5, 3204-3214.\n Head over to the\n discussions if\n you're interested and would like to contribute.\n Leave this as \"standard\" for now.\nqrrho : bool, optional\n Apply the quasi-rigid rotor harmonic oscillator (QRRHO) approximation of\n S. Grimme (see\n Theory. Chem. Eur. J., 2012, 18: 9955-9964)\n on top of the classical RRHO.\ntemperature : array-like, optional\n Absolute temperature in Kelvin.\npressure : array-like, optional\n Reference gas pressure.
\n\nReturns
\n\narray-like
\n\nExamples
\n\n\n
>>> import overreact as rx\n>>> from overreact import _constants as constants\n>>> model = rx.parse_model("data/ethane/B97-3c/model.k")\n>>> entropies = get_entropies(model.compounds)\n>>> (entropies - entropies.min()) / constants.calorie\narray([1.4, 0. ])\n
\n
You can consider all compounds as solvated if you want:
\n\n\n
>>> sol_entropies = get_entropies(model.compounds, environment="solvent")\n>>> (sol_entropies - entropies) / constants.calorie\narray([-6.35360874, -6.35360874])\n
\n
Generate an alternative scheme with some concentrations fixed.
\n\nThis function returns data that allow the microkinetic simulation of a\nreaction network under constraints, namely when some compounds have fixed\nconcentrations. This works by 1. removing all references to the fixed\ncompounds and by 2. properly multiplying the reaction rate constants by\nthe respective concentrations.
\n\nParameters
\n\nscheme : Scheme\n A descriptor of the reaction scheme.\n Mostly likely, this comes from a parsed model input file.\n See overreact.io.parse_model.\nk : array-like\n Reaction rate constant(s). Units match the concentration units given to\n the returned function dydt.\nfixed_y0 : dict-like\n Fixed initial state. Units match the concentration units given to\n the returned function dydt.
\n\nReturns
\n\nscheme : Scheme\n Associated reaction scheme with all references to fixed compounds\n removed.\nk : array-like\n Associated (effective) reaction rate constants that model the fixed\n concentrations.
\n\nNotes
\n\nKeep in mind that when a compound get its concentration fixed, the\nreaction scheme no longer conserves matter. You can think of it as\nreacting close to an infinite source of the compound, but it accumulates\nin the milleu at the given concentration.
\n\nExamples
\n\n\n
>>> import numpy as np\n>>> import overreact as rx\n
\n
Equilibria under a specific pH can be easily modeled:
\n\n\n
>>> pH = 7\n>>> scheme = rx.parse_reactions("AH <=> A- + H+")\n>>> k = np.array([1, 1])\n>>> scheme, k = rx.get_fixed_scheme(scheme, k, {"H+": 10**-pH})\n>>> scheme\nScheme(compounds=('AH', 'A-'),\n reactions=('AH -> A-',\n 'A- -> AH'),\n is_half_equilibrium=(True, True),\n A=((-1.0, 1.0),\n (1.0, -1.0)),\n B=((-1.0, 0.0),\n (1.0, 0.0)))\n>>> k\narray([1.e+00, 1.e-07])\n
\n
It is also possible to model the fixed activity of a solvent, for\ninstance:
\n\n\n
>>> scheme = rx.parse_reactions("A + 2H2O -> B")\n>>> k = np.array([1.0])\n>>> scheme, k = rx.get_fixed_scheme(scheme, k, {"H2O": 55.6})\n>>> scheme\nScheme(compounds=('A', 'B'),\n reactions=('A -> B',),\n is_half_equilibrium=(False,),\n A=((-1.0,),\n (1.0,)),\n B=((-1.0,),\n (1.0,)))\n>>> k\narray([3091.36])\n
\n
Multiple reactions work fine, see both examples below:
\n\n\n
>>> pH = 12\n>>> scheme = rx.parse_reactions("B <- AH <=> A- + H+")\n>>> k = np.array([10.0, 1, 1])\n>>> scheme, k = rx.get_fixed_scheme(scheme, k, {"H+": 10**-pH})\n>>> scheme\nScheme(compounds=('AH', 'B', 'A-'),\n reactions=('AH -> B',\n 'AH -> A-',\n 'A- -> AH'),\n is_half_equilibrium=(False, True, True),\n A=((-1.0, -1.0, 1.0),\n (1.0, 0.0, 0.0),\n (0.0, 1.0, -1.0)),\n B=((-1.0, -1.0, 0.0),\n (1.0, 0.0, 0.0),\n (0.0, 1.0, 0.0)))\n>>> k\narray([1.e+01, 1.e+00, 1.e-12])\n
\n
\n
>>> pH = 2\n>>> scheme = rx.parse_reactions(["AH <=> A- + H+", "B- + H+ <=> BH"])\n>>> k = np.array([1, 1, 2, 2])\n>>> scheme, k = rx.get_fixed_scheme(scheme, k, {"H+": 10**-pH})\n>>> scheme\nScheme(compounds=('AH', 'A-', 'B-', 'BH'),\n reactions=('AH -> A-',\n 'A- -> AH',\n 'B- -> BH',\n 'BH -> B-'),\n is_half_equilibrium=(True, True, True, True),\n A=((-1.0, 1.0, 0.0, 0.0),\n (1.0, -1.0, 0.0, 0.0),\n (0.0, 0.0, -1.0, 1.0),\n (0.0, 0.0, 1.0, -1.0)),\n B=((-1.0, 0.0, 0.0, 0.0),\n (1.0, 0.0, 0.0, 0.0),\n (0.0, 0.0, -1.0, 0.0),\n (0.0, 0.0, 1.0, 0.0)))\n>>> k\narray([1. , 0.01, 0.02, 2. ])\n
\n
Multiple fixed compounds also work fine:
\n\n\n
>>> pH = 6\n>>> scheme = rx.parse_reactions("A + H2O -> B <=> B- + H+")\n>>> k = np.array([1.0, 100.0, 2.0])\n>>> scheme, k = rx.get_fixed_scheme(scheme, k, {"H+": 10**-pH, "H2O": 55.6})\n>>> scheme\nScheme(compounds=('A', 'B', 'B-'),\n reactions=('A -> B',\n 'B -> B-',\n 'B- -> B'),\n is_half_equilibrium=(False, True, True),\n A=((-1.0, 0.0, 0.0),\n (1.0, -1.0, 1.0),\n (0.0, 1.0, -1.0)),\n B=((-1.0, 0.0, 0.0),\n (1.0, -1.0, 0.0),\n (0.0, 1.0, 0.0)))\n>>> k\narray([5.56e+01, 1.00e+02, 2.00e-06])\n
\n
This function is a no-op if fixed_y0 is empty, which is very important\nfor overall code consistency:
\n\n\n
>>> scheme = rx.parse_reactions(["AH <=> A- + H+", "B- + H+ <=> BH"])\n>>> k = np.array([1, 1, 2, 2])\n>>> new_scheme, new_k = rx.get_fixed_scheme(scheme, k, {})\n>>> new_scheme == scheme\nTrue\n>>> np.allclose(new_k, k)\nTrue\n
\n
Obtain free energies for compounds at a given temperature and pressure.
\n\nParameters
\n\ncompounds : dict-like\n A descriptor of the compounds.\n Mostly likely, this comes from a parsed input file.\n See overreact.io.parse_model.\nbias : array-like, optional\n Energy to be added to free energies.\nenvironment : str or None, optional\n Choose between \"gas\" and a solvent. This is chosen for you by default,\n based on the names of each compound. If given, all compounds will\n have the same behavior.\nmethod : str, optional\n This is a placeholder for future functionality.\n There are plans to implement more sophisticated methods for calculating\n entropies such as in\n Phys. Chem. Chem. Phys., 2019, 21, 18920-18929\n and\n J. Chem. Theory Comput. 2019, 15, 5, 3204-3214.\n Head over to the\n discussions if\n you're interested and would like to contribute.\n Leave this as \"standard\" for now.\nqrrho : bool or tuple-like, optional\n Apply both the quasi-rigid rotor harmonic oscillator (QRRHO)\n approximations of M. Head-Gordon and others (enthalpy correction, see\n J. Phys. Chem. C 2015, 119, 4, 1840-1850)\n and S. Grimme (entropy correction, see\n Theory. Chem. Eur. J., 2012, 18: 9955-9964)\n on top of the classical RRHO.\ntemperature : array-like, optional\n Absolute temperature in Kelvin.\npressure : array-like, optional\n Reference gas pressure.
\n\nReturns
\n\narray-like
\n\nExamples
\n\n\n
>>> import overreact as rx\n>>> from overreact import _constants as constants\n>>> model = rx.parse_model("data/ethane/B97-3c/model.k")\n>>> freeenergies = get_freeenergies(model.compounds, qrrho=(False, True))\n>>> (freeenergies - freeenergies.min()) / constants.kcal\narray([0. , 2.62281461])\n>>> freeenergies = get_freeenergies(model.compounds)\n>>> (freeenergies - freeenergies.min()) / constants.kcal\narray([0. , 2.62862818])\n
\n
You can consider all compounds as solvated if you want:
\n\n\n
>>> sol_freeenergies = get_freeenergies(model.compounds, environment="solvent")\n>>> (sol_freeenergies - freeenergies) / constants.kcal\narray([1.89432845, 1.89432845])\n
\n
You can set a simple energy bias, either as a constant or compound-wise:
\n\n\n
>>> get_freeenergies(model.compounds, bias=1.0) - freeenergies\narray([1., 1.])\n>>> get_freeenergies(model.compounds, bias=-1.0) - freeenergies\narray([-1., -1.])\n>>> get_freeenergies(model.compounds, bias=[1.0, -1.0]) - freeenergies\narray([ 1., -1.])\n
\n
Obtain internal energies for compounds at a given temperature.
\n\nParameters
\n\ncompounds : dict-like\n A descriptor of the compounds.\n Mostly likely, this comes from a parsed input file.\n See overreact.io.parse_model.\nqrrho : bool, optional\n Apply the quasi-rigid rotor harmonic oscillator (QRRHO) approximation of\n M. Head-Gordon and others (see\n J. Phys. Chem. C 2015, 119, 4, 1840-1850)\n on top of the classical RRHO.\ntemperature : array-like, optional\n Absolute temperature in Kelvin.
\n\nReturns
\n\narray-like
\n\nExamples
\n\n\n
>>> import overreact as rx\n>>> from overreact import _constants as constants\n>>> model = rx.parse_model("data/ethane/B97-3c/model.k")\n>>> internal_energies = get_internal_energies(model.compounds)\n>>> (internal_energies - internal_energies.min()) / constants.kcal\narray([0. , 2.20053981])\n
\n
Obtain reaction rate constants for a given reaction scheme.
\n\nParameters
\n\nscheme : Scheme\n A descriptor of the reaction scheme.\n Mostly likely, this comes from a parsed input file.\n See overreact.io.parse_model.\ncompounds : dict-like, optional\n A descriptor of the compounds.\n Mostly likely, this comes from a parsed input file.\n See overreact.io.parse_model.\nbias : array-like, optional\n Energy to be added to free energies.\ntunneling : str or None, optional\n Choose between \"eckart\", \"wigner\" or None (or \"none\").\nqrrho : bool or tuple-like, optional\n Apply both the quasi-rigid rotor harmonic oscillator (QRRHO)\n approximations of M. Head-Gordonand others (enthalpy correction, see\n J. Phys. Chem. C 2015, 119, 4, 1840-1850)\n and S. Grimme (entropy correction, see\n Theory. Chem. Eur. J., 2012, 18: 9955-9964)\n on top of the classical RRHO.\nscale : str, optional\n Reaction rate units. Possible values are \"cm3 mol-1 s-1\",\n \"l mol-1 s-1\", \"m3 mol-1 s-1\", \"cm3 particle-1 s-1\", \"mmHg-1 s-1\" and\n \"atm-1 s-1\".\ntemperature : array-like, optional\n Absolute temperature in Kelvin.\npressure : array-like, optional\n Reference gas pressure.\ndelta_freeenergies : array-like, optional\n Use this instead of obtaining delta free energies from the compounds.\nmolecularity : array-like, optional\n Reaction order, i.e., number of molecules that come together to react.\n If set, this is used to calculate delta_moles for\n overreact.thermo.equilibrium_constant, which effectively calculates a solution\n equilibrium constant between reactants and the transition state for\n gas phase data. You should set this to None if your free energies\n were already adjusted for solution Gibbs free energies.\nvolume : float, optional\n Molar volume.
\n\nReturns
\n\narray-like
\n\nNotes
\n\nSome symbols are accepted as alternatives in \"scale\": \"M-1\", \"ml\" and\n\"torr-1\" are understood as \"l mol-1\", \"cm3\" and \"mmHg-1\", respectively.
\n\nExamples
\n\nBelow is an example of an estimate for the rate of methyl rotation in\nethane (a trivial attempt to reproduce\nScience, 2006, 313, 5795, 1951-1955).\nHow many turns it does per second?
\n\n\n
>>> model = rx.parse_model("data/ethane/B97-3c/model.k")\n>>> get_k(model.scheme, model.compounds)\narray([8.16e+10])\n>>> get_k(model.scheme, model.compounds, qrrho=(False, True))\narray([8.24968117e+10])\n>>> get_k(model.scheme, model.compounds, qrrho=False)\narray([8.26909266e+10])\n>>> get_k(model.scheme, model.compounds, tunneling="wigner")\narray([7.99e+10])\n>>> get_k(model.scheme, model.compounds, tunneling=None)\narray([7.35e+10])\n
\n
The calculated value is off by less than 2% from the experimental value\n(\\( \\frac{1}{12 \\times 10^{-12}} \\text{s}^{-1} = 8.33 \\times 10^{10} \\text{s}^{-1} \\)).\nWe use Eckart tunneling by default, but see the effect of changing it\nabove.
\n\nThe units of the returned reaction rate constants can be selected for\nnon-unimolecular processes. The following is an attempt to reproduce\nJ Atmos Chem, 1996 23, 37-49 for\nthe reaction of proton-withdrawal by a chloride radical from the methane\nmolecule\n\\( \\ce{CH4 + \\cdot Cl -> [H3C\\cdots H\\cdots Cl]^\\ddagger -> H3C\\cdot + HCl} \\):
\n\n\n
>>> model = rx.parse_model("data/tanaka1996/UMP2/cc-pVTZ/model.jk")\n>>> get_k(model.scheme, model.compounds, temperature=300,\n... scale="cm3 particle-1 s-1")\narray([9.60e-14])\n
\n
(By the way, according to the Jet Propulsion Laboratory,\nPublication No. 19-5,\nthe experimental reaction rate constant for this reaction is\n\\( 1.0 \\times 10^{-13} \\text{cm}^3 \\text{particle}^{-1} \\text{s}^{-1} \\).)
\n\nThe returned units are \"M-1 s-1\" by default:
\n\n\n
>>> get_k(model.scheme, model.compounds) \\\n... == get_k(model.scheme, model.compounds, scale="l mol-1 s-1")\narray([ True])\n
\n
You can also turn the tunneling correction off by using the string \"none\":
\n\n\n
>>> get_k(model.scheme, model.compounds, tunneling="none") \\\n... == get_k(model.scheme, model.compounds, tunneling=None)\narray([ True])\n
\n
You can set a simple energy bias, either as a constant or compound-wise:
\n\n\n
>>> get_k(model.scheme, model.compounds, bias=1.0 * constants.kcal,\n... temperature=300.0, scale="cm3 particle-1 s-1")\narray([5.14e-13])\n>>> get_k(model.scheme, model.compounds,\n... bias=np.array([0.0, 0.0, -1.4, 0.0, 0.0]) * constants.kcal,\n... temperature=300.0, scale="cm3 particle-1 s-1")\narray([1.1e-12])\n
\n
Obtain tunneling transmission coefficients at a given temperature.
\n\nOne tunneling transmission coefficient is calculated for each reaction. If\na reaction lacks a transition state (i.e., a half-equilibrium reaction),\nits transmission coefficient is set to unity.
\n\nParameters
\n\nscheme : Scheme\n A descriptor of the reaction scheme.\n Mostly likely, this comes from a parsed input file.\n See overreact.io.parse_model.\ncompounds : dict-like\n A descriptor of the compounds.\n Mostly likely, this comes from a parsed input file.\n See overreact.io.parse_model.\nmethod : str or None, optional\n Choose between \"eckart\", \"wigner\" or None (or \"none\").\nqrrho : bool, optional\n Apply both the quasi-rigid rotor harmonic oscillator (QRRHO)\n approximations of M. Head-Gordon and others (enthalpy correction, see\n J. Phys. Chem. C 2015, 119, 4, 1840-1850)\n and S. Grimme (entropy correction, see\n Theory. Chem. Eur. J., 2012, 18: 9955-9964)\n on top of the classical RRHO.\ntemperature : array-like, optional\n Absolute temperature in Kelvin.
\n\nReturns
\n\narray-like
\n\nRaises
\n\nValueError\n If method is not supported.
\n\nExamples
\n\nBelow is an example of an estimate of how much quantum tunneling\ncontributes to the rate of methyl rotation in ethane (see\nScience, 2006, 313, 5795, 1951-1955\nfor some interesting experimental data on this reaction).
\n\n\n
>>> model = rx.parse_model("data/ethane/B97-3c/model.k")\n>>> kappa = get_kappa(model.scheme, model.compounds)\n>>> kappa\narray([1.110])\n>>> get_kappa(model.scheme, model.compounds, method="none")\narray([1.0])\n>>> get_kappa(model.scheme, model.compounds, method="none") \\\n... == get_kappa(model.scheme, model.compounds, method=None)\narray([ True])\n
\n
You can calculate each piece of the reaction rate constant by hand,\nif you want. Just make sure that you don't calculate the tunneling\ncoefficient twice:
\n\n\n
>>> kappa * get_k(model.scheme, model.compounds, tunneling=None)\narray([8.e+10])\n
\n
Calculate entropy contributions from the overall reaction structure.
\n\nThis function currently implements the reaction translational entropy, a\nresult of the indistinguishability of reactants or products, i.e., a\ndifference in entropy of \\( R ln 2! \\) for the reactions
\n\n$$\\ce{A + B -> C}$$
\n\nand
\n\n$$\\ce{2A -> C}$$
\n\nParameters
\n\ntransform : array-like\ntemperature : array-like, optional\n Absolute temperature in Kelvin.\npressure : array-like, optional\n Reference gas pressure.
\n\nReturns
\n\ndelta_entropy : array-like
\n\nExamples
\n\n\n
>>> import overreact as rx\n
\n
\n
>>> scheme = rx.parse_reactions("A + B <=> C")\n>>> get_reaction_entropies(scheme.A)\narray([0.0, 0.0])\n>>> scheme = rx.parse_reactions("2A <=> C")\n>>> get_reaction_entropies(scheme.A)\narray([-5.763, 5.763])\n>>> get_reaction_entropies(scheme.A, temperature=400.0)\narray([-5.763, 5.763])\n
\n
\n
>>> scheme = rx.parse_reactions("A <=> B + C")\n>>> get_reaction_entropies(scheme.A)\narray([0.0, 0.0])\n>>> scheme = rx.parse_reactions("A <=> 2B")\n>>> get_reaction_entropies(scheme.A)\narray([ 5.763, -5.763])\n
\n
\n
>>> scheme = rx.parse_reactions("A + B -> E# -> C")\n>>> get_reaction_entropies(scheme.A)\narray([0.0])\n>>> get_reaction_entropies(scheme.B)\narray([0.0])\n>>> scheme = rx.parse_reactions("2A -> E# -> C")\n>>> get_reaction_entropies(scheme.A)\narray([-5.763])\n>>> get_reaction_entropies(scheme.B)\narray([-5.763])\n
\n
Return the indices of transition states for each reaction.
\n\nParameters
\n\nA, B : array-like\nis_half_equilibrium : sequence
\n\nReturns
\n\nsequence
\n\nExamples
\n\n\n
>>> scheme = parse_reactions("A -> B")\n>>> print(scheme)\nScheme(compounds=('A', 'B'),\n reactions=('A -> B',),\n is_half_equilibrium=(False,),\n A=((-1.,), (1.,)),\n B=((-1.,), (1.,)))\n>>> get_transition_states(scheme.A, scheme.B, scheme.is_half_equilibrium)\n(None,)\n
\n
\n
>>> scheme = parse_reactions("S -> E\u2021 -> S")\n>>> print(scheme)\nScheme(compounds=('S', 'E\u2021'),\n reactions=('S -> S',),\n is_half_equilibrium=(False,),\n A=((0.,), (0.,)),\n B=((-1.,), (1.,)))\n>>> get_transition_states(scheme.A, scheme.B, scheme.is_half_equilibrium)\n(1,)\n
\n
\n
>>> scheme = parse_reactions("E + S <=> ES -> ES\u2021 -> E + P")\n>>> print(scheme)\nScheme(compounds=('E', 'S', 'ES', 'ES\u2021', 'P'),\n reactions=('E + S -> ES', 'ES -> E + S', 'ES -> E + P'),\n is_half_equilibrium=(True, True, False),\n A=((-1., 1., 1.),\n (-1., 1., 0.),\n (1., -1., -1.),\n (0., 0., 0.),\n (0., 0., 1.)),\n B=((-1., 0., 0.),\n (-1., 0., 0.),\n (1., 0., -1.),\n (0., 0., 1.),\n (0., 0., 0.)))\n>>> get_transition_states(scheme.A, scheme.B, scheme.is_half_equilibrium)\n(None, None, 3)\n
\n
Simulate a reaction scheme from its rate function.
\n\nThis function provides two functions that calculate the concentrations and\nthe rates of formation at any point in time for any compound. It does that\nby solving an initial value problem (IVP) through scipy's solve_ivp\nunder the hood.
\n\nParameters
\n\ndydt : callable\n Right-hand side of the system.\ny0 : array-like\n Initial state.\nt_span : array-like, optional\n Interval of integration (t0, tf). The solver starts with t=t0 and\n integrates until it reaches t=tf. If not given, a conservative value\n is chosen based on the system at hand (the method of choice works for\n any zeroth-, first- or second-order reactions).\nmethod : str, optional\n Integration method to use. See scipy.integrate.solve_ivp for details.\n Kinetics problems are very often stiff and, as such, \"RK23\" and \"RK45\" may be\n unsuited. \"LSODA\", \"BDF\", and \"Radau\" are worth a try if things go bad.\nmax_step : float, optional\n Maximum step to be performed by the integrator.\n Defaults to half the total time span.\nfirst_step : float, optional\n First step size.\n Defaults to half the maximum step, or np.finfo(np.float64).eps,\n whichever is smallest.\nrtol, atol : array-like, optional\n See scipy.integrate.solve_ivp for details.\nmax_time : float, optional\n If t_span is not given, an interval will be estimated, but it can't\n be larger than this parameter.
\n\nReturns
\n\ny, r : callable\n Concentrations and reaction rates as functions of time. The y object\n is an OdeSolution and stores attributes t_min and t_max.
\n\nExamples
\n\n\n
>>> import numpy as np\n>>> import overreact as rx\n
\n
A toy simulation can be performed in just two lines:
\n\n\n
>>> scheme = rx.parse_reactions("A <=> B")\n>>> y, r = get_y(get_dydt(scheme, np.array([1, 1])), y0=[1, 0])\n
\n
The y object stores information about the simulation time, which can be\nused to produce a suitable vector of timepoints for, e.g., plotting:
\n\n\n
>>> y.t_min, y.t_max\n(0.0, 3.0)\n>>> t = np.linspace(y.t_min, y.t_max)\n>>> t\narray([0. , 0.06122449, ..., 2.93877551, 3. ])\n
\n
Both y and r can be used to check concentrations and rates in any\npoint in time. In particular, both are vectorized:
\n\n\n
>>> y(t)\narray([[1. , ...],\n [0. , ...]])\n>>> r(t)\narray([[-1. , ...],\n [ 1. , ...]])\n
\n
Check whether a name specifies a transition state.
\n\nParameters
\n\nname : str
\n\nReturns
\n\nbool
\n\nExamples
\n\n\n
>>> is_transition_state("A#")\nTrue\n>>> is_transition_state("A\u2021")\nTrue\n>>> is_transition_state("pyrrole#")\nTrue\n>>> is_transition_state("pyrrole\u2021")\nTrue\n>>> is_transition_state("A")\nFalse\n>>> is_transition_state("A~")\nFalse\n>>> is_transition_state("pyrrole")\nFalse\n>>> is_transition_state("pyrrole~")\nFalse\n
\n
This function also works for names that specify environment:
\n\n\n
>>> is_transition_state("A#(w)")\nTrue\n>>> is_transition_state("A\u2021(w)")\nTrue\n>>> is_transition_state("TS#(w)")\nTrue\n>>> is_transition_state("TS\u2021(w)")\nTrue\n>>> is_transition_state("A(w)")\nFalse\n>>> is_transition_state("A~(w)")\nFalse\n>>> is_transition_state("TS(w)")\nFalse\n>>> is_transition_state("TS~(w)")\nFalse\n
\n
Parse a set of compounds.
\n\nParameters
\n\ntext : str, sequence of str or dict-like\n Compound descriptions or sequence of lines of it.\npath : sequence of str, optional\n Paths for searching logfiles.\nselect : sequence of str, optional\n If defined, only those compounds will be returned.
\n\nReturns
\n\ncompounds : immutable dict-like
\n\nRaises
\n\nFileNotFoundError\n If a logfile is not found.
\n\nExamples
\n\n\n
>>> import overreact as rx\n
\n
\n
>>> compounds = rx.parse_compounds("S: data/ethane/B97-3c/staggered.out")\n>>> compounds\n{'S': {'logfile': 'data/ethane/B97-3c/staggered.out',\n 'energy': -209483812.77142256,\n 'mult': 1,\n 'atomnos': (6, 6, 1, 1, 1, 1, 1, 1),\n ...}}\n
\n
\n
>>> compounds = rx.parse_compounds({"S": "data/ethane/B97-3c/staggered.out"})\n>>> compounds\n{'S': {'logfile': 'data/ethane/B97-3c/staggered.out',\n 'energy': -209483812.77142256,\n 'mult': 1,\n 'atomnos': (6, 6, 1, 1, 1, 1, 1, 1),\n ...}}\n
\n
\n
>>> compounds = rx.parse_compounds(["S: data/ethane/B97-3c/staggered.out",\n... "E\u2021: data/ethane/B97-3c/eclipsed.out"])\n>>> compounds\n{'S': {'logfile': 'data/ethane/B97-3c/staggered.out',\n 'energy': -209483812.77142256,\n 'mult': 1,\n 'atomnos': (6, 6, 1, 1, 1, 1, 1, 1),\n ...},\n'E\u2021': {'logfile': 'data/ethane/B97-3c/eclipsed.out',\n 'energy': -209472585.3539883,\n 'mult': 1,\n 'atomnos': (6, 6, 1, 1, 1, 1, 1, 1),\n ...}}\n
\n
\n
>>> compounds = rx.parse_compounds('''S: data/ethane/B97-3c/staggered.out\n... E\u2021: data/ethane/B97-3c/eclipsed.out''')\n>>> compounds\n{'S': {'logfile': 'data/ethane/B97-3c/staggered.out',\n 'energy': -209483812.77142256,\n 'mult': 1,\n 'atomnos': (6, 6, 1, 1, 1, 1, 1, 1),\n ...},\n'E\u2021': {'logfile': 'data/ethane/B97-3c/eclipsed.out',\n 'energy': -209472585.3539883,\n 'mult': 1,\n 'atomnos': (6, 6, 1, 1, 1, 1, 1, 1),\n ...}}\n
\n
\n
>>> compounds = rx.parse_compounds({"S": "data/ethane/B97-3c/staggered.out",\n... "E\u2021": "data/ethane/B97-3c/eclipsed.out"})\n>>> compounds\n{'S': {'logfile': 'data/ethane/B97-3c/staggered.out',\n 'energy': -209483812.77142256,\n 'mult': 1,\n 'atomnos': (6, 6, 1, 1, 1, 1, 1, 1),\n ...},\n'E\u2021': {'logfile': 'data/ethane/B97-3c/eclipsed.out',\n 'energy': -209472585.3539883,\n 'mult': 1,\n 'atomnos': (6, 6, 1, 1, 1, 1, 1, 1),\n ...}}\n
\n
Parse either a source or model input file, whichever is available.
\n\nA source input file (also known as a .k file) contains all the information needed\nto create a model input file.\nA model input file (also known as a .jk file) is a JSON encoded file with all the\ninformation needed to study microkinetic simulations from first principles.
\n\nYou probably won't need to use model input files directly, they are\nautomatically created based on source input files.\nTake a look at our guide on how to write an source input file.
\n\nThis function attempts to parse a model input file if available. If not, a source\ninput file is parsed and a model input file is generated from it. Extensions are\nguessed if none given (i.e., if only the base name given).
\n\nParameters
\n\npath : str\n Path to the model or source input file.\n If the final extension is not .jk or .k, it is guessed.\nforce_compile : bool\n If True, a .k file will take precedence over any .jk file for reading. A\n .jk file is thus either generated or overwritten. This is sometimes\n needed to force an update with new data.
\n\nReturns
\n\nmodel : immutable dict-like
\n\nRaises
\n\nFileNotFoundError\n If the model or source input file is not found.
\n\nExamples
\n\nSome examples of how overreact \"sees\" your data below \ud83d\ude04:
\n\n\n
>>> model = parse_model("data/ethane/B97-3c/model.jk")\n>>> model.scheme\nScheme(compounds=('S', 'E\u2021'),\n reactions=('S -> S',),\n is_half_equilibrium=(False,),\n A=((0.0,), (0.0,)),\n B=((-1.0,), (1.0,)))\n>>> model.compounds["S"]\n{'logfile': 'data/ethane/B97-3c/staggered.out',\n 'energy': -209483812.77142256,\n 'mult': 1,\n 'atomnos': (6, 6, 1, 1, 1, 1, 1, 1),\n 'atommasses': (12.011, 12.011, 1.008, 1.008, 1.008, 1.008, 1.008, 1.008),\n 'atomcoords': ((-7.633588, 2.520693, -4.8e-05),\n ...,\n (-5.832852, 3.674431, 0.363239)),\n 'vibfreqs': (307.57, 825.42, ..., 3071.11, 3071.45),\n 'vibdisps': (((-1.7e-05, 3.4e-05, 5.4e-05),\n ...,\n (-0.011061, -0.030431, -0.027036)))}\n>>> model_from_source = parse_model("data/ethane/B97-3c/model.k",\n... force_compile=True)\n>>> model_from_source == model\nTrue\n>>> model_from_source = parse_model("data/ethane/B97-3c/model")\n>>> model_from_source == model\nTrue\n
\n
Parse a kinetic model as a chemical reaction scheme.
\n\nThis is an essential part of the parsing process.\nSee overreact.io.parse_model other details.
\n\nParameters
\n\ntext : str or sequence of str\n Model description or sequence of lines of it.
\n\nReturns
\n\nscheme : Scheme\n A descriptor of the reaction scheme.
\n\nNotes
\n\nThe model description should comply with the mini-language for systems of\nreactions. A semi-formal definition of the grammar in\nBackus-Naur form\nis given below:
\n\n equation ::= equation_side arrow equation_side\nequation_side ::= coefficient compound ['+' coefficient compound]*\n coefficient ::= [integers] (defaults to 1)\n compound ::= mix of printable characters\n arrow ::= '->' | '<=>' | '<-'\n
\n\nBlank lines and comments (starting with //) are ignored. Repeated\nreactions are ignored. Furthermore, reactions can be chained one after\nanother and, if a single compound (with either a \u2021 or a # at the end)\nappears alone on one side of a reaction, it's considered a transition\nstate. Transition states have zero lifetime during the simulation.
\n\nExamples
\n\nWhat follows is a rather long tour over the parsing process and its\noutput in general. You can skip it if you are not interested in the\ndetails.
\n\n\n
>>> scheme = parse_reactions("A -> B // a direct reaction")\n
\n
The reaction above is a direct one (observe that comments are ignored). The\nreturned object has the following attributes:
\n\n\n
>>> scheme.compounds\n('A', 'B')\n>>> scheme.reactions\n('A -> B',)\n>>> scheme.is_half_equilibrium\n(False,)\n>>> scheme.A\n((-1.,), (1.,))\n>>> scheme.B\n((-1.,), (1.,))\n
\n
The same reaction can be specified in reverse order:
\n\n\n
>>> parse_reactions("B <- A // reverse reaction of the above")\nScheme(compounds=('A', 'B'),\n reactions=('A -> B',),\n is_half_equilibrium=(False,),\n A=((-1.,), (1.,)),\n B=((-1.,), (1.,)))\n
\n
Equilibria produce twice as many direct reactions, while the $B$ matrix\ndefines an energy relationship for only one of each pair:
\n\n\n
>>> parse_reactions("A <=> B // an equilibrium")\nScheme(compounds=('A', 'B'),\n reactions=('A -> B', 'B -> A'),\n is_half_equilibrium=(True, True),\n A=((-1., 1.),\n (1., -1.)),\n B=((-1., 0.),\n (1., 0.)))\n
\n
Adding twice the same reaction results in a single reaction being added.\nThis of course also works with equilibria (extra whitespaces are ignored):
\n\n\n
>>> parse_reactions('''\n... A <=> B -> A\n... A -> B <=> A\n... A -> B <- A\n... B <- A -> B\n... ''')\nScheme(compounds=('A', 'B'),\n reactions=('A -> B', 'B -> A'),\n is_half_equilibrium=(True, True),\n A=((-1., 1.),\n (1., -1.)),\n B=((-1., 0.),\n (1., 0.)))\n
\n
Transition states are specified with a special symbol at the end (either\n\u2021 or #). They are shown among compounds, but the matrix $A$ ensures\nthey'll never have a non-zero rate of formation/consumption. On the other\nhand, they are needed in the $B$ matrix:
\n\n\n
>>> parse_reactions("A -> A\u2021 -> B")\nScheme(compounds=('A', 'A\u2021', 'B'),\n reactions=('A -> B',),\n is_half_equilibrium=(False,),\n A=((-1.,), (0.,), (1.,)),\n B=((-1.,), (1.,), (0.,)))\n
\n
This gives the same result as above:
\n\n\n
>>> parse_reactions("A -> A\u2021 -> B <- A\u2021 <- A")\nScheme(compounds=('A', 'A\u2021', 'B'),\n reactions=('A -> B',),\n is_half_equilibrium=(False,),\n A=((-1.,), (0.,), (1.,)),\n B=((-1.,), (1.,), (0.,)))\n
\n
It is possible to define a reaction whose product is the same as the\nreactant. This is found in isomerization processes (e.g., ammonia\ninversion or the methyl rotation in ethane):
\n\n\n
>>> parse_reactions("S -> E\u2021 -> S")\nScheme(compounds=('S', 'E\u2021'),\n reactions=('S -> S',),\n is_half_equilibrium=(False,),\n A=((0.,), (0.,)),\n B=((-1.,), (1.,)))\n
\n
As such, a column full of zeros in the $A$ matrix corresponds to a reaction\nwith zero net change. As can be seen, overreact allows for very general\nmodels. An interesting feature is that a single transition state can link\nmany different compounds (whether it is useful is a matter of debate):
\n\n\n
>>> parse_reactions('''\n... B -> B\u2021 -> C // chained reactions and transition states\n... B\u2021 -> D // this is a bifurcation\n... B -> B'\u2021 -> E // this is a classical competitive reaction\n... A -> B\u2021\n... ''')\nScheme(compounds=('B', 'B\u2021', 'C', 'D', "B'\u2021", 'E', 'A'),\n reactions=('B -> C', 'B -> D', 'B -> E', 'A -> C', 'A -> D'),\n is_half_equilibrium=(False, False, False, False, False),\n A=((-1., -1., -1., 0., 0.),\n (0., 0., 0., 0., 0.),\n (1., 0., 0., 1., 0.),\n (0., 1., 0., 0., 1.),\n (0., 0., 0., 0., 0.),\n (0., 0., 1., 0., 0.),\n (0., 0., 0., -1., -1.)),\n B=((-1., -1., -1., 0., 0.),\n (1., 1., 0., 1., 1.),\n (0., 0., 0., 0., 0.),\n (0., 0., 0., 0., 0.),\n (0., 0., 1., 0., 0.),\n (0., 0., 0., 0., 0.),\n (0., 0., 0., -1., -1.)))\n
\n
The following is a borderline case but both reactions should be considered\ndifferent since they define different processes:
\n\n\n
>>> parse_reactions('''\n... A -> A\u2021 -> B\n... A -> B\n... ''')\nScheme(compounds=('A', 'A\u2021', 'B'),\n reactions=('A -> B', 'A -> B'),\n is_half_equilibrium=(False, False),\n A=((-1., -1.),\n (0., 0.),\n (1., 1.)),\n B=((-1., -1.),\n (1., 0.),\n (0., 1.)))\n
\n
The following is correct behavior. In fact, the reactions are badly\ndefined: if more than one transition state are chained, the following\nhappens, which is correct since it's the most physically plausible model\nthat can be extracted. It can be seen as a feature that the product B is\nignored and not the reactant A, since the user would easily see the mistake\nin graphs of concentration over time (the alternative would be no\nreaction happening at all, which is rather cryptic to debug).
\n\n\n
>>> parse_reactions("A -> A\u2021 -> A'\u2021 -> B")\nScheme(compounds=('A', 'A\u2021', "A'\u2021", 'B'),\n reactions=("A -> A'\u2021",),\n is_half_equilibrium=(False,),\n A=((-1.,), (0.,), (1.,), (0.,)),\n B=((-1.,), (1.,), (0.,), (0.,)))\n
\n
In any case, it's not clear how a reaction barrier be defined in such a\ncase. If you have a use case, don't hesitate to\nopen an issue, we'll be\nhappy to hear from you.
\n", "signature": "(text: Union[str, Sequence[str]]) -> overreact.core.Scheme:", "funcdef": "def"}, "overreact.unparse_reactions": {"fullname": "overreact.unparse_reactions", "modulename": "overreact", "qualname": "unparse_reactions", "kind": "function", "doc": "Unparse a kinetic model.
\n\nParameters
\n\nscheme : Scheme\n A descriptor of the reaction scheme.\n Mostly likely, this comes from a parsed input file.\n See overreact.io.parse_model.
\n\nReturns
\n\ntext : str
\n\nNotes
\n\nThis function assumes complimentary half equilibria are located one after\nthe other in scheme.reactions, which is to be expected from\nparse_reactions.
\n\nExamples
\n\n\n
>>> unparse_reactions(Scheme(compounds=('A', 'B'), reactions=('A -> B',),\n... is_half_equilibrium=(False,),\n... A=((-1.,),\n... ( 1.,)),\n... B=((-1.,),\n... ( 1.,))))\n'A -> B'\n>>> unparse_reactions(Scheme(compounds=('A', 'B'),\n... reactions=('A -> B', 'B -> A'),\n... is_half_equilibrium=(True, True),\n... A=((-1., 1.),\n... ( 1., -1.)),\n... B=((-1., 0.),\n... ( 1., 0.))))\n'A <=> B'\n>>> unparse_reactions(Scheme(compounds=('A', 'A\u2021', 'B'),\n... reactions=('A -> B',),\n... is_half_equilibrium=(False,),\n... A=((-1.,),\n... ( 0.,),\n... ( 1.,)),\n... B=((-1.,),\n... ( 1.,),\n... ( 0.,))))\n'A -> A\u2021 -> B'\n>>> print(unparse_reactions(Scheme(compounds=('B', 'B\u2021', 'C', 'D', "B'\u2021",\n... 'E', 'A'),\n... reactions=('B -> C', 'B -> D', 'B -> E',\n... 'A -> C', 'A -> D'),\n... is_half_equilibrium=(False, False, False,\n... False, False),\n... A=((-1., -1., -1., 0., 0.),\n... ( 0., 0., 0., 0., 0.),\n... ( 1., 0., 0., 1., 0.),\n... ( 0., 1., 0., 0., 1.),\n... ( 0., 0., 0., 0., 0.),\n... ( 0., 0., 1., 0., 0.),\n... ( 0., 0., 0., -1., -1.)),\n... B=((-1., -1., -1., 0., 0.),\n... ( 1., 1., 0., 1., 1.),\n... ( 0., 0., 0., 0., 0.),\n... ( 0., 0., 0., 0., 0.),\n... ( 0., 0., 1., 0., 0.),\n... ( 0., 0., 0., 0., 0.),\n... ( 0., 0., 0., -1., -1.)))))\nB -> B\u2021 -> C\nB -> B\u2021 -> D\nB -> B'\u2021 -> E\nA -> B\u2021 -> C\nA -> B\u2021 -> D\n>>> print(unparse_reactions(Scheme(compounds=('A', 'A\u2021', 'B'),\n... reactions=('A -> B', 'A -> B'),\n... is_half_equilibrium=(False, False),\n... A=((-1., -1.),\n... ( 0., 0.),\n... ( 1., 1.)),\n... B=((-1., -1.),\n... ( 1., 0.),\n... ( 0., 1.)))))\nA -> A\u2021 -> B\nA -> B\n>>> unparse_reactions(Scheme(compounds=('A', 'A\u2021', "A'\u2021", 'B'),\n... reactions=("A -> A'\u2021",),\n... is_half_equilibrium=(False,),\n... A=((-1.,),\n... ( 0.,),\n... ( 1.,),\n... ( 0.,)),\n... B=((-1.,),\n... ( 1.,),\n... ( 0.,),\n... ( 0.,))))\n"A -> A\u2021 -> A'\u2021"\n>>> unparse_reactions(Scheme(compounds=('S', 'E\u2021'),\n... reactions=('S -> S',),\n... is_half_equilibrium=(False,),\n... A=((0.0,),\n... (0.0,)),\n... B=((-1.0,),\n... (1.0,))))\n'S -> E\u2021 -> S'\n
\n
+
+
![\"PyPI\"](\"https://img.shields.io/pypi/v/overreact\")
\n![\"Python](\"https://img.shields.io/pypi/pyversions/overreact\")
\n![\"CI\"](\"https://github.com/geem-lab/overreact/actions/workflows/python-package.yml/badge.svg\"){kind=icon}
\n\n![\"Coverage\"](\"https://codecov.io/gh/geem-lab/overreact/branch/main/graph/badge.svg?token=4WAVXCRXY8\"){kind=icon}
\n![\"License\"](\"https://img.shields.io/github/license/geem-lab/overreact\")
\n![\"User](\"https://img.shields.io/badge/user%20guide-available-blue\"){kind=icon}
\n![\"GitHub](\"https://img.shields.io/github/discussions/geem-lab/overreact\")
\n![\"GitHub](\"https://img.shields.io/github/issues-raw/geem-lab/overreact\")
\n![\"DOI\"](\"https://img.shields.io/badge/DOI-10.1002%2Fjcc.26861-blue\"){kind=icon}
\n![\"DOI\"](\"https://zenodo.org/badge/DOI/10.5281/zenodo.7504800.svg\"){kind=icon}
\n![\"Made](\"https://img.shields.io/badge/made%20in-Brazil-009c3b\"){kind=icon}
\n![\"DOI\"](\"https://zenodo.org/badge/DOI/10.5281/zenodo.7865357.svg\"){kind=icon}
\n