overreact

19class Scheme(NamedTuple):
20    """
21    A descriptor of a chemical reaction network.
22
23    Mostly likely, this comes from a parsed input file.
24    See `overreact.io.parse_model`.
25    """
26
27    compounds: list[str]
28    """A descriptor of compounds."""
29    reactions: list[str]
30    """A descriptor of reactions."""
31    is_half_equilibrium: list[bool]
32    """An indicator of whether a reaction is half-equilibrium."""
33    A: np.ndarray
34    """A matrix of stoichiometric coefficients between reactants and products."""
35    B: np.ndarray
36    """A matrix of stoichiometric coefficients between reactants and transition states."""

759def change_reference_state(
760    new_reference: float = 1.0 / constants.liter,
761    old_reference: float | None = None,
762    sign: int = 1,
763    temperature: float | np.ndarray = 298.15,
764    pressure: float = constants.atm,
765    volume: float | None = None,
766):
767    r"""Calculate an additive entropy correction to a change in reference states.
768
769    .. math::
770        \Delta G_\text{corr} =
771            R T \ln \left( \frac{\chi_\text{new}}{\chi_\text{old}} \right)
772
773    The value returned can be directly multiplied by temperature and summed to
774    the old reference free energies to obtain free energies with respect to a
775    new reference. See notes below.
776
777    For instance, the concentration correction to Gibbs free energy for a
778    gas-to-liquid standard state change is simply
779    (:math:`c^\circ = \frac{\text{1 atm}}{R T}`),
780
781    .. math::
782        \Delta G_\text{conc} =
783            R T \ln \left( \frac{\text{1 M}}{c^\circ} \right)
784
785    Parameters
786    ----------
787    new_reference : array-like, optional
788        New reference state. Default value corresponds to 1 mol/liter.
789    old_reference : array-like, optional
790        Old reference state. Default value corresponds to the concentration of
791        an ideal gas at the given temperature and 1 atm.
792    sign : float, optional
793        Sign of the change in reference state. Default value is 1. This only
794        multiplies the final result.
795    temperature : array-like, optional
796        Absolute temperature in Kelvin.
797    pressure : array-like, optional
798        Reference gas pressure.
799    volume : float, optional
800        Molar volume.
801
802    Returns
803    -------
804    correction : array-like
805        Entropy correction in J/mol·K.
806
807    Notes
808    -----
809    This function can be used to add any entropy correction in the form above.
810    The only drawback is that, sometimes, those corrections are written with a
811    minus sign in front of them (this implies switching the roles of
812    `old_reference` and `new_reference`). The easiest way to accomplish this is
813    by using ``sign=-1`` or multiplying the result by ``-1``.
814
815    Examples
816    --------
817    By default, the correction returns a change in concentration from the gas
818    phase standard concentration to the solvated-state standard concentration:
819
820    >>> -rx.change_reference_state() / constants.calorie
821    -6.4
822    >>> 298.15 * rx.change_reference_state() / constants.kcal
823    1.89
824    >>> 273.15 * rx.change_reference_state(temperature=273.15) / constants.kcal
825    1.69
826
827    But this function can also be used to adjust symmetry effects from C1
828    calculations (symmetry number equals to one). For D7h, for instance, the
829    symmetry number is 14:
830
831    >>> -298.15 * rx.change_reference_state(14, 1) / constants.kcal
832    -1.56
833
834    >>> rx.change_reference_state(sign=-1) == -rx.change_reference_state()
835    True
836
837    """
838    temperature = np.asarray(temperature)
839
840    if old_reference is None:
841        if volume is None:
842            volume = molar_volume(temperature=temperature, pressure=pressure)
843        old_reference = 1.0 / volume
844
845    return sign * constants.R * np.log(new_reference / old_reference)

578def get_bias(
579    scheme,
580    compounds,
581    data,
582    y0,
583    tunneling="eckart",
584    qrrho=True,
585    temperature=298.15,
586    pressure=constants.atm,
587    method="RK23",
588    rtol=1e-3,
589    atol=1e-6,
590):
591    r"""Estimate a energy bias for a given set of reference data points.
592
593    Parameters
594    ----------
595    scheme : Scheme
596        A descriptor of the reaction scheme.
597        Mostly likely, this comes from a parsed model input file.
598        See `overreact.io.parse_model`.
599    compounds : dict-like
600        A descriptor of the compounds.
601        Mostly likely, this comes from a parsed model input file.
602        See `overreact.io.parse_model`.
603    data : dict-like of array-like
604    y0: array-like
605    tunneling : str or None, optional
606        Choose between "eckart", "wigner" or None (or "none").
607    qrrho : bool or tuple-like, optional
608        Apply both the quasi-rigid rotor harmonic oscillator (QRRHO)
609        approximations of M. Head-Gordon and others (enthalpy correction, see
610        [*J. Phys. Chem. C* **2015**, 119, 4, 1840-1850](http://dx.doi.org/10.1021/jp509921r)) and S. Grimme (entropy correction, see
611        [*Theory. Chem. Eur. J.*, **2012**, 18: 9955-9964](https://doi.org/10.1002/chem.201200497)) on top of the classical RRHO.
612    temperature : array-like, optional
613        Absolute temperature in Kelvin.
614    pressure : array-like, optional
615        Reference gas pressure.
616    delta_freeenergies : array-like, optional
617        Use this instead of obtaining delta free energies from the compounds.
618    molecularity : array-like, optional
619        Reaction order, i.e., number of molecules that come together to react.
620        If set, this is used to calculate `delta_moles` for
621        `equilibrium_constant`, which effectively calculates a solution
622        equilibrium constant between reactants and the transition state for
623        gas phase data. You should set this to `None` if your free energies
624        were already adjusted for solution Gibbs free energies.
625    volume : float, optional
626        Molar volume.
627
628    Returns
629    -------
630    array-like
631
632    Examples
633    --------
634    >>> model = rx.parse_model("data/tanaka1996/UMP2/cc-pVTZ/model.jk")
635
636    The following are some estimates on actual atmospheric concentrations:
637
638    >>> y0 = [4.8120675684099e-5,
639    ...       2.8206357713029e-5,
640    ...       0.0,
641    ...       0.0,
642    ...       2.7426565371219e-5]
643    >>> data = {"t": [1.276472128376942246e-6,
644    ...               1.446535794555581743e-4,
645    ...               1.717069678525567564e-2],
646    ...         "CH3·": [9.694916853338366211e-9,
647    ...                  1.066033349343709026e-6,
648    ...                  2.632179124780495175e-5]}
649    >>> get_bias(model.scheme, model.compounds, data, y0) / constants.kcal
650    -1.4
651    """
652    max_time = np.max(data["t"])
653
654    def f(bias):
655        k = rx.get_k(
656            scheme,
657            compounds,
658            bias=bias,
659            tunneling=tunneling,
660            qrrho=qrrho,
661            temperature=temperature,
662            pressure=pressure,
663        )
664
665        # TODO(schneiderfelipe): support schemes with fixed concentrations
666        dydt = rx.get_dydt(scheme, k)
667        y, _ = rx.get_y(
668            dydt,
669            y0=y0,
670            method=method,
671            rtol=rtol,
672            atol=atol,
673            max_time=max_time,
674        )
675
676        yhat = y(data["t"])
677        return np.sum(
678            [
679                (yhat[i] - data[name]) ** 2
680                for (i, name) in enumerate(compounds)
681                if name in data
682            ],
683        )
684
685    res = minimize_scalar(f)
686    return res.x

572def get_delta(transform, property):
573    """Calculate deltas according to reactions.
574
575    Delta properties are differences in a property between the final and
576    initial state of a chemical transformation. They are calculated from
577    matrices representing the transformation and the absolute properties.
578    Transformation matrices are expected to have column-wise transformation
579    defined.
580
581    Very useful for the calculation of reaction and activation free energies
582    from absolute free energies of compounds. Matrices ``A`` and ``B`` of a
583    `Scheme` represent the transformations associated with reaction and
584    activation free energies, respectively.
585
586    Parameters
587    ----------
588    transform : array-like
589    property : array-like
590
591    Returns
592    -------
593    delta_property : array-like
594
595    Examples
596    --------
597    >>> get_delta([-1, 1], [-10, 10])
598    20
599
600    You must ensure the transformation is properly defined, as no test is made
601    to ensure, e.g., conservation of matter:
602
603    >>> get_delta([-1, 0], [-10, 20])
604    10
605
606    Normally, transformations are given as columns in a matrix:
607
608    >>> get_delta([[-1, -2],
609    ...            [ 1,  3]], [-5, 12])
610    array([17, 46])
611    """
612    return np.asarray(transform).T @ np.asarray(property)

200def get_dydt(scheme, k, ef=EF):
201    """Generate a rate function that models a reaction scheme.
202
203    Parameters
204    ----------
205    scheme : Scheme
206        A descriptor of the reaction scheme.
207        Mostly likely, this comes from a parsed model input file.
208        See `overreact.io.parse_model`.
209    k : array-like
210        Reaction rate constant(s). Units match the concentration units given to
211        the returned function ``dydt``.
212    ef : float, optional
213        Equilibrium factor. This is a parameter that can be used to scale the
214        reaction rates associated to half-equilibrium reactions such that they
215        are faster than the other reactions.
216
217    Returns
218    -------
219    dydt : callable
220        Reaction rate function. The actual reaction rate constants employed
221        are stored in the attribute `k` of the returned function. If JAX is
222        available, the attribute `jac` will hold the Jacobian function of
223        `dydt`.
224
225    Notes
226    -----
227    The returned function is suited to be used by ODE solvers such as
228    `scipy.integrate.solve_ivp` or the older `scipy.integrate.ode` (see
229    examples below). This is actually what the function `get_y` from the
230    current module does.
231
232    Examples
233    --------
234    >>> import numpy as np
235    >>> import overreact as rx
236
237    >>> scheme = rx.parse_reactions("A <=> B")
238    >>> dydt = get_dydt(scheme, np.array([1, 1]))
239    >>> dydt(0.0, np.array([1., 1.]))
240    Array([0., 0.], ...)
241
242    If available, JAX is used for JIT compilation. This will make `dydt`
243    complain if given lists instead of numpy arrays. So stick to the safer,
244    faster side as above.
245
246    The actually used reaction rate constants can be inspected with the `k`
247    attribute of `dydt`:
248
249    >>> dydt.k
250    Array([1., 1.], ...)
251
252    If JAX is available, the Jacobian function will be available as
253    `dydt.jac`:
254
255    >>> dydt.jac(0.0, np.array([1., 1.]))
256    Array([[-1.,  1.],
257           [ 1., -1.]], ...)
258
259    """
260    scheme = rx.core._check_scheme(scheme)
261    A = jnp.asarray(scheme.A)
262    M = jnp.where(A > 0, 0, -A).T
263    k_adj = _adjust_k(scheme, k, ef=ef)
264
265    def _dydt(_t, y):
266        r = k_adj * jnp.prod(jnp.power(y, M), axis=1)
267        return jnp.dot(A, r)
268
269    if _found_jax:
270        # Using JAX for JIT compilation is much faster.
271        _dydt = jit(_dydt)
272
273        # NOTE(schneiderfelipe): the following function is defined
274        # such that _jac(t, y)[i, j] == d f_i / d y_j,
275        # with shape of (n_compounds, n_compounds).
276        def _jac(t, y):
277            logger.warning(f"\x1b[A@t = \x1b[94m{t:10.3f} \x1b[ms\x1b[K")
278            return jacfwd(lambda _y: _dydt(t, _y))(y)
279
280        _dydt.jac = _jac
281
282    _dydt.k = k_adj
283    return _dydt

 97def get_enthalpies(
 98    compounds: dict,
 99    qrrho: bool = True,
100    temperature: float = 298.15,
101):
102    """Obtain enthalpies for compounds at a given temperature.
103
104    Parameters
105    ----------
106    compounds : dict-like
107        A descriptor of the compounds.
108        Mostly likely, this comes from a parsed input file.
109        See `overreact.io.parse_model`.
110    qrrho : bool, optional
111        Apply the quasi-rigid rotor harmonic oscillator (QRRHO) approximation of
112        M. Head-Gordon and others (see
113        [*J. Phys. Chem. C* **2015**, 119, 4, 1840-1850](http://dx.doi.org/10.1021/jp509921r))
114        on top of the classical RRHO.
115    temperature : array-like, optional
116        Absolute temperature in Kelvin.
117
118    Returns
119    -------
120    array-like
121
122    Examples
123    --------
124    >>> import overreact as rx
125    >>> from overreact import _constants as constants
126    >>> model = rx.parse_model("data/ethane/B97-3c/model.k")
127    >>> enthalpies = get_enthalpies(model.compounds)
128    >>> (enthalpies - enthalpies.min()) / constants.kcal
129    array([0. , 2.20053981])
130
131    The enthalpies at absolute zero can easily be obtained (this is used,
132    e.g., in the calculation of the Eckart tunneling coefficient, see
133    `overreact.tunnel.eckart`). We can use this to calculate, for instance,
134    the thermal contributions to the enthalpy:
135
136    >>> zero_enthalpies = get_enthalpies(model.compounds, temperature=0)
137    >>> (enthalpies - zero_enthalpies) / constants.kcal
138    array([2.78, 2.50])
139    """
140    compounds = rx.io._check_compounds(compounds)
141    enthalpies = []
142    for name in compounds:
143        logger.info(f"calculate enthalpy: {name}")
144
145        # TODO(schneiderfelipe): inertia might benefit from caching
146        moments, _, _ = coords.inertia(
147            compounds[name].atommasses,
148            compounds[name].atomcoords,
149        )
150
151        enthalpy = rx.thermo.calc_enthalpy(
152            energy=compounds[name].energy,
153            degeneracy=compounds[name].mult,
154            moments=moments,
155            vibfreqs=compounds[name].vibfreqs,
156            qrrho=qrrho,
157            temperature=temperature,
158        )
159        enthalpies.append(enthalpy)
160    return np.array(enthalpies)

163def get_entropies(
164    compounds: dict,
165    environment: str | None = None,
166    method: str = "standard",
167    qrrho: bool = True,
168    temperature: float = 298.15,
169    pressure: float = constants.atm,
170):
171    """Obtain entropies for compounds at a given temperature and pressure.
172
173    Parameters
174    ----------
175    compounds : dict-like
176        A descriptor of the compounds.
177        Mostly likely, this comes from a parsed input file.
178        See `overreact.io.parse_model`.
179    environment : str or None, optional
180        Choose between "gas" and a solvent. This is chosen for you by default,
181        based on the names of each compound (e.g. `A(g)` or `A` is gas,
182        `A(w)` or `A(...)` is solvated). In case this is given, all compounds
183        will have the same behavior.
184    method : str, optional
185        This is a placeholder for future functionality.
186        There are plans to implement more sophisticated methods for calculating
187        entropies such as in
188        [*Phys. Chem. Chem. Phys.*, **2019**, 21, 18920-18929](https://doi.org/10.1039/C9CP03226F)
189        and
190        [*J. Chem. Theory Comput.* **2019**, 15, 5, 3204-3214](https://doi.org/10.1021/acs.jctc.9b00214).
191        Head over to the
192        [discussions](https://github.com/geem-lab/overreact/discussions) if
193        you're interested and would like to contribute.
194        Leave this as "standard" for now.
195    qrrho : bool, optional
196        Apply the quasi-rigid rotor harmonic oscillator (QRRHO) approximation of
197        S. Grimme (see
198        [*Theory. Chem. Eur. J.*, **2012**, 18: 9955-9964](https://doi.org/10.1002/chem.201200497))
199        on top of the classical RRHO.
200    temperature : array-like, optional
201        Absolute temperature in Kelvin.
202    pressure : array-like, optional
203        Reference gas pressure.
204
205    Returns
206    -------
207    array-like
208
209    Examples
210    --------
211    >>> import overreact as rx
212    >>> from overreact import _constants as constants
213    >>> model = rx.parse_model("data/ethane/B97-3c/model.k")
214    >>> entropies = get_entropies(model.compounds)
215    >>> (entropies - entropies.min()) / constants.calorie
216    array([1.4, 0. ])
217
218    You can consider all compounds as solvated if you want:
219
220    >>> sol_entropies = get_entropies(model.compounds, environment="solvent")
221    >>> (sol_entropies - entropies) / constants.calorie
222    array([-6.35360874, -6.35360874])
223    """
224    compounds = rx.io._check_compounds(compounds)
225    entropies = []
226    for name in compounds:
227        logger.info(f"calculate entropy: {name}")
228
229        if "point_group" in compounds[name]:
230            point_group = compounds[name].point_group
231        else:
232            point_group = coords.find_point_group(
233                compounds[name].atommasses,
234                compounds[name].atomcoords,
235            )
236        symmetry_number = coords.symmetry_number(point_group)
237
238        # TODO(schneiderfelipe): inertia might benefit from caching
239        moments, _, _ = coords.inertia(
240            compounds[name].atommasses,
241            compounds[name].atomcoords,
242        )
243
244        if environment is None:
245            environment = rx.core._get_environment(name)
246        entropy = rx.thermo.calc_entropy(
247            atommasses=compounds[name].atommasses,
248            atomcoords=compounds[name].atomcoords,
249            energy=compounds[name].energy,
250            degeneracy=compounds[name].mult,
251            moments=moments,
252            symmetry_number=symmetry_number,
253            vibfreqs=compounds[name].vibfreqs,
254            environment=environment,
255            method=method,
256            qrrho=qrrho,
257            temperature=temperature,
258            pressure=pressure,
259        )
260
261        if compounds[name].symmetry is not None:
262            # The negative sign here seems correct. See equations (9) and (10)
263            # of doi:10.1002/qua.25686.
264            entropy -= rx.change_reference_state(
265                compounds[name].symmetry,
266                1,
267                temperature=temperature,
268                pressure=pressure,
269            )
270
271        entropies.append(entropy)
272    return np.array(entropies)

365def get_fixed_scheme(scheme, k, fixed_y0):
366    """Generate an alternative scheme with some concentrations fixed.
367
368    This function returns data that allow the microkinetic simulation of a
369    reaction network under constraints, namely when some compounds have fixed
370    concentrations. This works by 1. removing all references to the fixed
371    compounds and by 2. properly multiplying the reaction rate constants by
372    the respective concentrations.
373
374    Parameters
375    ----------
376    scheme : Scheme
377        A descriptor of the reaction scheme.
378        Mostly likely, this comes from a parsed model input file.
379        See `overreact.io.parse_model`.
380    k : array-like
381        Reaction rate constant(s). Units match the concentration units given to
382        the returned function ``dydt``.
383    fixed_y0 : dict-like
384        Fixed initial state. Units match the concentration units given to
385        the returned function ``dydt``.
386
387    Returns
388    -------
389    scheme : Scheme
390        Associated reaction scheme with all references to fixed compounds
391        removed.
392    k : array-like
393        Associated (effective) reaction rate constants that model the fixed
394        concentrations.
395
396    Notes
397    -----
398    Keep in mind that when a compound get its concentration fixed, the
399    reaction scheme no longer conserves matter. You can think of it as
400    reacting close to an infinite source of the compound, but it accumulates
401    in the milleu at the given concentration.
402
403    Examples
404    --------
405    >>> import numpy as np
406    >>> import overreact as rx
407
408    Equilibria under a specific pH can be easily modeled:
409
410    >>> pH = 7
411    >>> scheme = rx.parse_reactions("AH <=> A- + H+")
412    >>> k = np.array([1, 1])
413    >>> scheme, k = rx.get_fixed_scheme(scheme, k, {"H+": 10**-pH})
414    >>> scheme
415    Scheme(compounds=('AH', 'A-'),
416           reactions=('AH -> A-',
417                      'A- -> AH'),
418           is_half_equilibrium=(True, True),
419           A=((-1.0, 1.0),
420              (1.0, -1.0)),
421           B=((-1.0, 0.0),
422              (1.0, 0.0)))
423    >>> k
424    array([1.e+00, 1.e-07])
425
426    It is also possible to model the fixed activity of a solvent, for
427    instance:
428
429    >>> scheme = rx.parse_reactions("A + 2H2O -> B")
430    >>> k = np.array([1.0])
431    >>> scheme, k = rx.get_fixed_scheme(scheme, k, {"H2O": 55.6})
432    >>> scheme
433    Scheme(compounds=('A', 'B'),
434           reactions=('A -> B',),
435           is_half_equilibrium=(False,),
436           A=((-1.0,),
437              (1.0,)),
438           B=((-1.0,),
439              (1.0,)))
440    >>> k
441    array([3091.36])
442
443    Multiple reactions work fine, see both examples below:
444
445    >>> pH = 12
446    >>> scheme = rx.parse_reactions("B <- AH <=> A- + H+")
447    >>> k = np.array([10.0, 1, 1])
448    >>> scheme, k = rx.get_fixed_scheme(scheme, k, {"H+": 10**-pH})
449    >>> scheme
450    Scheme(compounds=('AH', 'B', 'A-'),
451           reactions=('AH -> B',
452                      'AH -> A-',
453                      'A- -> AH'),
454           is_half_equilibrium=(False, True, True),
455           A=((-1.0, -1.0, 1.0),
456              (1.0, 0.0, 0.0),
457              (0.0, 1.0, -1.0)),
458           B=((-1.0, -1.0, 0.0),
459              (1.0, 0.0, 0.0),
460              (0.0, 1.0, 0.0)))
461    >>> k
462    array([1.e+01, 1.e+00, 1.e-12])
463
464    >>> pH = 2
465    >>> scheme = rx.parse_reactions(["AH <=> A- + H+", "B- + H+ <=> BH"])
466    >>> k = np.array([1, 1, 2, 2])
467    >>> scheme, k = rx.get_fixed_scheme(scheme, k, {"H+": 10**-pH})
468    >>> scheme
469    Scheme(compounds=('AH', 'A-', 'B-', 'BH'),
470           reactions=('AH -> A-',
471                      'A- -> AH',
472                      'B- -> BH',
473                      'BH -> B-'),
474           is_half_equilibrium=(True, True, True, True),
475           A=((-1.0, 1.0, 0.0, 0.0),
476              (1.0, -1.0, 0.0, 0.0),
477              (0.0, 0.0, -1.0, 1.0),
478              (0.0, 0.0, 1.0, -1.0)),
479           B=((-1.0, 0.0, 0.0, 0.0),
480              (1.0, 0.0, 0.0, 0.0),
481              (0.0, 0.0, -1.0, 0.0),
482              (0.0, 0.0, 1.0, 0.0)))
483    >>> k
484    array([1. , 0.01, 0.02, 2. ])
485
486    Multiple fixed compounds also work fine:
487
488    >>> pH = 6
489    >>> scheme = rx.parse_reactions("A + H2O -> B <=> B- + H+")
490    >>> k = np.array([1.0, 100.0, 2.0])
491    >>> scheme, k = rx.get_fixed_scheme(scheme, k, {"H+": 10**-pH, "H2O": 55.6})
492    >>> scheme
493    Scheme(compounds=('A', 'B', 'B-'),
494           reactions=('A -> B',
495                      'B -> B-',
496                      'B- -> B'),
497           is_half_equilibrium=(False, True, True),
498           A=((-1.0, 0.0, 0.0),
499              (1.0, -1.0, 1.0),
500              (0.0, 1.0, -1.0)),
501           B=((-1.0, 0.0, 0.0),
502              (1.0, -1.0, 0.0),
503              (0.0, 1.0, 0.0)))
504    >>> k
505    array([5.56e+01, 1.00e+02, 2.00e-06])
506
507    This function is a no-op if `fixed_y0` is empty, which is very important
508    for overall code consistency:
509
510    >>> scheme = rx.parse_reactions(["AH <=> A- + H+", "B- + H+ <=> BH"])
511    >>> k = np.array([1, 1, 2, 2])
512    >>> new_scheme, new_k = rx.get_fixed_scheme(scheme, k, {})
513    >>> new_scheme == scheme
514    True
515    >>> np.allclose(new_k, k)
516    True
517
518    """
519    new_k = np.asarray(k, dtype=np.float64).copy()
520    new_reactions = []
521    for i, (reaction, is_half_equilibrium) in enumerate(
522        zip(scheme.reactions, scheme.is_half_equilibrium),
523    ):
524        for reactants, products, _ in rx.core._parse_reactions(
525            reaction,
526        ):
527            new_reactants = tuple(
528                (coeff, compound)
529                for (coeff, compound) in reactants
530                if compound not in fixed_y0
531            )
532            new_products = tuple(
533                (coeff, compound)
534                for (coeff, compound) in products
535                if compound not in fixed_y0
536            )
537
538            for fixed_compound in fixed_y0:
539                for coeff, compound in reactants:
540                    if fixed_compound == compound:
541                        new_k[i] *= fixed_y0[fixed_compound] ** coeff
542
543            new_reactions.append((new_reactants, new_products, is_half_equilibrium))
544
545    new_reactions = tuple(r for r in rx.core._unparse_reactions(new_reactions))
546    new_is_half_equilibrium = scheme.is_half_equilibrium
547
548    new_A = []
549    new_B = []
550    new_compounds = []
551    for compound, row_A, row_B in zip(
552        scheme.compounds,
553        scheme.A,
554        scheme.B,
555    ):
556        if compound not in fixed_y0:
557            new_compounds.append(compound)
558            new_A.append(row_A)
559            new_B.append(row_B)
560
561    new_compounds = tuple(new_compounds)
562    new_A = tuple(new_A)
563    new_B = tuple(new_B)
564
565    return (
566        rx.core.Scheme(
567            compounds=new_compounds,
568            reactions=new_reactions,
569            is_half_equilibrium=new_is_half_equilibrium,
570            A=new_A,
571            B=new_B,
572        ),
573        new_k,
574    )

321def get_freeenergies(
322    compounds: dict,
323    bias: float = 0.0,
324    environment: str | None = None,
325    method: str = "standard",
326    qrrho: bool | tuple[bool, bool] = True,
327    temperature: float = 298.15,
328    pressure: float = constants.atm,
329):
330    """Obtain free energies for compounds at a given temperature and pressure.
331
332    Parameters
333    ----------
334    compounds : dict-like
335        A descriptor of the compounds.
336        Mostly likely, this comes from a parsed input file.
337        See `overreact.io.parse_model`.
338    bias : array-like, optional
339        Energy to be added to free energies.
340    environment : str or None, optional
341        Choose between "gas" and a solvent. This is chosen for you by default,
342        based on the names of each compound. If given, all compounds will
343        have the same behavior.
344    method : str, optional
345        This is a placeholder for future functionality.
346        There are plans to implement more sophisticated methods for calculating
347        entropies such as in
348        [*Phys. Chem. Chem. Phys.*, **2019**, 21, 18920-18929](https://doi.org/10.1039/C9CP03226F)
349        and
350        [*J. Chem. Theory Comput.* **2019**, 15, 5, 3204-3214](https://doi.org/10.1021/acs.jctc.9b00214).
351        Head over to the
352        [discussions](https://github.com/geem-lab/overreact/discussions) if
353        you're interested and would like to contribute.
354        Leave this as "standard" for now.
355    qrrho : bool or tuple-like, optional
356        Apply both the quasi-rigid rotor harmonic oscillator (QRRHO)
357        approximations of M. Head-Gordon and others (enthalpy correction, see
358        [*J. Phys. Chem. C* **2015**, 119, 4, 1840-1850](http://dx.doi.org/10.1021/jp509921r))
359        and S. Grimme (entropy correction, see
360        [*Theory. Chem. Eur. J.*, **2012**, 18: 9955-9964](https://doi.org/10.1002/chem.201200497))
361        on top of the classical RRHO.
362    temperature : array-like, optional
363        Absolute temperature in Kelvin.
364    pressure : array-like, optional
365        Reference gas pressure.
366
367    Returns
368    -------
369    array-like
370
371    Examples
372    --------
373    >>> import overreact as rx
374    >>> from overreact import _constants as constants
375    >>> model = rx.parse_model("data/ethane/B97-3c/model.k")
376    >>> freeenergies = get_freeenergies(model.compounds, qrrho=(False, True))
377    >>> (freeenergies - freeenergies.min()) / constants.kcal
378    array([0. , 2.62281461])
379    >>> freeenergies = get_freeenergies(model.compounds)
380    >>> (freeenergies - freeenergies.min()) / constants.kcal
381    array([0. , 2.62862818])
382
383    You can consider all compounds as solvated if you want:
384
385    >>> sol_freeenergies = get_freeenergies(model.compounds, environment="solvent")
386    >>> (sol_freeenergies - freeenergies) / constants.kcal
387    array([1.89432845, 1.89432845])
388
389    You can set a simple energy bias, either as a constant or compound-wise:
390
391    >>> get_freeenergies(model.compounds, bias=1.0) - freeenergies
392    array([1., 1.])
393    >>> get_freeenergies(model.compounds, bias=-1.0) - freeenergies
394    array([-1., -1.])
395    >>> get_freeenergies(model.compounds, bias=[1.0, -1.0]) - freeenergies
396    array([ 1., -1.])
397    """
398    qrrho_enthalpy, qrrho_entropy = _check_qrrho(qrrho)
399    enthalpies = get_enthalpies(
400        compounds,
401        qrrho=qrrho_enthalpy,
402        temperature=temperature,
403    )
404    entropies = get_entropies(
405        compounds,
406        environment=environment,
407        method=method,
408        qrrho=qrrho_entropy,
409        temperature=temperature,
410        pressure=pressure,
411    )
412    # TODO(schneiderfelipe): log the contribution of bias
413    return enthalpies - temperature * entropies + np.asarray(bias)

39def get_internal_energies(
40    compounds: dict,
41    qrrho: bool = True,
42    temperature: float = 298.15,
43):
44    """Obtain internal energies for compounds at a given temperature.
45
46    Parameters
47    ----------
48    compounds : dict-like
49        A descriptor of the compounds.
50        Mostly likely, this comes from a parsed input file.
51        See `overreact.io.parse_model`.
52    qrrho : bool, optional
53        Apply the quasi-rigid rotor harmonic oscillator (QRRHO) approximation of
54        M. Head-Gordon and others (see
55        [*J. Phys. Chem. C* **2015**, 119, 4, 1840-1850](http://dx.doi.org/10.1021/jp509921r))
56        on top of the classical RRHO.
57    temperature : array-like, optional
58        Absolute temperature in Kelvin.
59
60    Returns
61    -------
62    array-like
63
64    Examples
65    --------
66    >>> import overreact as rx
67    >>> from overreact import _constants as constants
68    >>> model = rx.parse_model("data/ethane/B97-3c/model.k")
69    >>> internal_energies = get_internal_energies(model.compounds)
70    >>> (internal_energies - internal_energies.min()) / constants.kcal
71    array([0. , 2.20053981])
72
73    """
74    compounds = rx.io._check_compounds(compounds)
75    internal_energies = []
76    for name in compounds:
77        logger.info(f"calculate internal energy: {name}")
78
79        # TODO(schneiderfelipe): inertia might benefit from caching
80        moments, _, _ = coords.inertia(
81            compounds[name].atommasses,
82            compounds[name].atomcoords,
83        )
84
85        internal_energy = rx.thermo.calc_internal_energy(
86            energy=compounds[name].energy,
87            degeneracy=compounds[name].mult,
88            moments=moments,
89            vibfreqs=compounds[name].vibfreqs,
90            qrrho=qrrho,
91            temperature=temperature,
92        )
93        internal_energies.append(internal_energy)
94    return np.array(internal_energies)

416def get_k(
417    scheme: Scheme,
418    compounds: dict | None = None,
419    bias: float = 0.0,
420    tunneling: str = "eckart",
421    qrrho: bool | tuple[bool, bool] = True,
422    scale: str = "l mol-1 s-1",
423    temperature: float = 298.15,
424    pressure: float = constants.atm,
425    delta_freeenergies: float | None = None,
426    molecularity: float | None = None,
427    volume: float | None = None,
428) -> float:
429    r"""Obtain reaction rate constants for a given reaction scheme.
430
431    Parameters
432    ----------
433    scheme : Scheme
434        A descriptor of the reaction scheme.
435        Mostly likely, this comes from a parsed input file.
436        See `overreact.io.parse_model`.
437    compounds : dict-like, optional
438        A descriptor of the compounds.
439        Mostly likely, this comes from a parsed input file.
440        See `overreact.io.parse_model`.
441    bias : array-like, optional
442        Energy to be added to free energies.
443    tunneling : str or None, optional
444        Choose between "eckart", "wigner" or None (or "none").
445    qrrho : bool or tuple-like, optional
446        Apply both the quasi-rigid rotor harmonic oscillator (QRRHO)
447        approximations of M. Head-Gordonand others (enthalpy correction, see
448        [*J. Phys. Chem. C* **2015**, 119, 4, 1840-1850](http://dx.doi.org/10.1021/jp509921r))
449        and S. Grimme (entropy correction, see
450        [*Theory. Chem. Eur. J.*, **2012**, 18: 9955-9964](https://doi.org/10.1002/chem.201200497))
451        on top of the classical RRHO.
452    scale : str, optional
453        Reaction rate units. Possible values are "cm3 mol-1 s-1",
454        "l mol-1 s-1", "m3 mol-1 s-1", "cm3 particle-1 s-1", "mmHg-1 s-1" and
455        "atm-1 s-1".
456    temperature : array-like, optional
457        Absolute temperature in Kelvin.
458    pressure : array-like, optional
459        Reference gas pressure.
460    delta_freeenergies : array-like, optional
461        Use this instead of obtaining delta free energies from the compounds.
462    molecularity : array-like, optional
463        Reaction order, i.e., number of molecules that come together to react.
464        If set, this is used to calculate `delta_moles` for
465        `overreact.thermo.equilibrium_constant`, which effectively calculates a solution
466        equilibrium constant between reactants and the transition state for
467        gas phase data. You should set this to `None` if your free energies
468        were already adjusted for solution Gibbs free energies.
469    volume : float, optional
470        Molar volume.
471
472    Returns
473    -------
474    array-like
475
476    Notes
477    -----
478    Some symbols are accepted as alternatives in "scale": "M-1", "ml" and
479    "torr-1" are understood as "l mol-1", "cm3" and "mmHg-1", respectively.
480
481    Examples
482    --------
483    Below is an example of an estimate for the rate of methyl rotation in
484    ethane (a trivial attempt to reproduce
485    [*Science*, **2006**, 313, 5795, 1951-1955](https://doi.org/10.1126/science.1132178)).
486    **How many turns it does per second?**
487
488    >>> model = rx.parse_model("data/ethane/B97-3c/model.k")
489    >>> get_k(model.scheme, model.compounds)
490    array([8.16e+10])
491    >>> get_k(model.scheme, model.compounds, qrrho=(False, True))
492    array([8.24968117e+10])
493    >>> get_k(model.scheme, model.compounds, qrrho=False)
494    array([8.26909266e+10])
495    >>> get_k(model.scheme, model.compounds, tunneling="wigner")
496    array([7.99e+10])
497    >>> get_k(model.scheme, model.compounds, tunneling=None)
498    array([7.35e+10])
499
500    The calculated value is off by less than 2% from the experimental value
501    (:math:`\frac{1}{12 \times 10^{-12}} \text{s}^{-1} = 8.33 \times 10^{10} \text{s}^{-1}`).
502    We use Eckart tunneling by default, but see the effect of changing it
503    above.
504
505    The units of the returned reaction rate constants can be selected for
506    non-unimolecular processes. The following is an attempt to reproduce
507    [*J Atmos Chem*, **1996** 23, 37-49](https://doi.org/10.1007/BF00058703) for
508    the reaction of proton-withdrawal by a chloride radical from the methane
509    molecule
510    :math:`\ce{CH4 + \cdot Cl -> [H3C\cdots H\cdots Cl]^\ddagger -> H3C\cdot + HCl}`:
511
512    >>> model = rx.parse_model("data/tanaka1996/UMP2/cc-pVTZ/model.jk")
513    >>> get_k(model.scheme, model.compounds, temperature=300,
514    ...       scale="cm3 particle-1 s-1")
515    array([9.60e-14])
516
517    (By the way, according to the Jet Propulsion Laboratory,
518    [Publication No. 19-5](https://jpldataeval.jpl.nasa.gov/),
519    the experimental reaction rate constant for this reaction is
520    :math:`1.0 \times 10^{-13} \text{cm}^3 \text{particle}^{-1} \text{s}^{-1}`.)
521
522    The returned units are "M-1 s-1" by default:
523
524    >>> get_k(model.scheme, model.compounds) \
525    ... == get_k(model.scheme, model.compounds, scale="l mol-1 s-1")
526    array([ True])
527
528    You can also turn the tunneling correction off by using the string "none":
529
530    >>> get_k(model.scheme, model.compounds, tunneling="none") \
531    ... == get_k(model.scheme, model.compounds, tunneling=None)
532    array([ True])
533
534    You can set a simple energy bias, either as a constant or compound-wise:
535
536    >>> get_k(model.scheme, model.compounds, bias=1.0 * constants.kcal,
537    ...       temperature=300.0, scale="cm3 particle-1 s-1")
538    array([5.14e-13])
539    >>> get_k(model.scheme, model.compounds,
540    ...       bias=np.array([0.0, 0.0, -1.4, 0.0, 0.0]) * constants.kcal,
541    ...       temperature=300.0, scale="cm3 particle-1 s-1")
542    array([1.1e-12])
543    """
544    qrrho_enthalpy, qrrho_entropy = _check_qrrho(qrrho)
545    scheme = rx.core._check_scheme(scheme)
546    if compounds is not None:
547        compounds = rx.io._check_compounds(compounds)
548    if delta_freeenergies is None:
549        assert compounds is not None, "compounds could not be inferred"
550        freeenergies = get_freeenergies(
551            compounds,
552            bias=bias,
553            qrrho=(qrrho_enthalpy, qrrho_entropy),
554            # NOTE(schneiderfelipe): ensure we get rate constants in M-1 s-1.
555            # TODO(schneiderfelipe): this strategy will have to change
556            # somewhat when we improve solvation entropy models.
557            environment="solvent",
558            temperature=temperature,
559            pressure=pressure,
560        )
561
562        # TODO(schneiderfelipe): log the contribution of reaction symmetry
563        delta_freeenergies = rx.get_delta(
564            scheme.B,
565            freeenergies,
566        ) - temperature * rx.get_reaction_entropies(
567            scheme.B,
568            temperature=temperature,
569            pressure=pressure,
570        )
571
572    if molecularity is None:
573        molecularity = rx.thermo.get_molecularity(scheme.A)
574
575    # NOTE(schneiderfelipe): passing molecularity here to rates.eyring messes up
576    # rate constant units (by a factor of M-1 s-1 to atm-1 s-1), so we leave it as is.
577    k = rates.eyring(
578        delta_freeenergies,
579        temperature=temperature,
580        pressure=pressure,
581        volume=volume,
582    )
583
584    # make reaction rate constants for equilibria as close as possible to one
585    i = 0
586    while i < len(scheme.is_half_equilibrium):
587        if scheme.is_half_equilibrium[i]:
588            pair = k[i : i + 2]
589            _K = pair[0] / pair[1]
590
591            denom = pair.min()
592            if denom == 0.0:
593                logger.warning(
594                    "found half-equilibrium reaction with zero rate constant: skipping equilibrium normalization",
595                )
596                denom = 1.0
597
598            k[i : i + 2] = pair / denom
599            assert np.isclose(_K, k[i] / k[i + 1]), (
600                f"reaction rate constants {k[i]} and {k[i + 1]} for "
601                "equilibria could not be made to match the expected "
602                f"equilibrium constant value {_K}"
603            )
604
605            # loop over pairs of equilibria
606            i += 1
607        i += 1
608
609    logger.info(
610        "(classical) reaction rate constants: "
611        f"{', '.join([f'{v:7.3g}' for v in k])} atm⁻ⁿ⁺¹·s⁻¹",
612    )
613    if tunneling not in {"none", None}:
614        if compounds is not None:
615            kappa = get_kappa(
616                scheme,
617                compounds,
618                method=tunneling,
619                qrrho=qrrho_enthalpy,
620                temperature=temperature,
621            )
622            k *= kappa
623        else:
624            # TODO(schneiderfelipe): when get_kappa accept deltas, this will
625            # be probably unnecessary.
626            logger.warning(
627                "assuming unitary tunneling coefficients due to incomplete compound data",
628            )
629        logger.info(
630            "(tunneling) reaction rate constants: "
631            f"{', '.join([f'{v:7.3g}' for v in k])} atm⁻ⁿ⁺¹·s⁻¹",
632        )
633
634    # TODO(schneiderfelipe): ensure diffusional limit for reactions in
635    # solvation using Collins-Kimball theory. This includes half-equilibria.
636    return rates.convert_rate_constant(
637        k,
638        new_scale=scale,
639        # NOTE(schneiderfelipe): all the code above should always generate
640        # rate constants in M-1 s-1, so that we convert them here.
641        old_scale="l mol-1 s-1",
642        molecularity=molecularity,
643        temperature=temperature,
644        pressure=pressure,
645    )