geem-lab / overreact Goto Github PK

I encountered some numerical problems when attempting to solve a large system (observe the initial concentration, 1e-5):

$ overreact data/perez-soto2020/RI/BLYP-D4/def2-TZVP/model.k \
    "Benzaldehyde(dcm): 1e-5" "NButylamine(dcm): 1e-5" --plot

gives rise to (observe the explosion in the benzaldehyde concentration and the negative hemiaminal concentration, among other discrepancies):

  no   compound            t = 0 s   t = 1e+11 s
 ────────────────────────────────────────────────
   0   Benzaldehyde(dcm)     0.000     17223.494
   1   NButylamine(dcm)      0.000         0.039
   2   A_N(dcm)              0.000         0.022
   3   A_N_N(dcm)            0.000         0.000
   4   Water(dcm)            0.000        -0.001
   5   A_N_W(dcm)            0.000        -0.000
   6   A_N_N_W(dcm)          0.000        -0.000
   7   A_N_W_W(dcm)          0.000         0.000
   8   TS1_#(dcm)            0.000         0.000
   9   Hemiaminal(dcm)       0.000        -8.834
  10   TS2_#(dcm)            0.000         0.000
  11   I_W(dcm)              0.000        -0.000
  12   TS1N_#(dcm)           0.000         0.000
  13   Int_N(dcm)            0.000         0.000
  14   TS2N_#(dcm)           0.000         0.000
  15   I_N_W(dcm)            0.000        -0.000
  16   TS1W_#(dcm)           0.000         0.000
  17   Int_W(dcm)            0.000        -0.000
  18   TS2W_#(dcm)           0.000         0.000
  19   I_W_W(dcm)            0.000         0.000
  20   TS1WW_#(dcm)          0.000         0.000
  21   Int_W_W(dcm)          0.000         0.000
  22   TS2WW_#(dcm)          0.000         0.000
  23   I_W_W_W(dcm)          0.000        -0.000
  24   Imine(dcm)            0.000         1.273
  25   I_N_W_W(dcm)          0.000         0.000

and the following profile:

When using 1e-6 as concentration, reaction time becomes sane (0.0001 s versus 1e+11 s), plus the graph becomes flat (no reaction, that is; final concentrations are precisely the same as in the beginning).

When using 1e-4 as concentration, I get an overflow:

/home/schneider/.local/lib/python3.8/site-packages/scipy/integrate/_ivp/common.py:336: RuntimeWarning: overflow encountered in multiply
  new_factor = NUM_JAC_FACTOR_INCREASE * factor[ind]
/home/schneider/.local/lib/python3.8/site-packages/scipy/integrate/_ivp/common.py:358: RuntimeWarning: overflow encountered in multiply
  factor[max_diff < NUM_JAC_DIFF_SMALL * scale] *= NUM_JAC_FACTOR_INCREASE

Discussion

The above is a problem with the Jacobian, which we currently numerically approximate (EDIT: more precisely, the initial rate is probably poorly estimated). That is not ideal, so I propose a couple of things to mitigate this problem:

• Allow the user to choose among ODE solvers (see method in scipy.integrate.solve_ivp).
• Allow the user to select an initial time step.
• Calculate an analytical Jacobian and use it.
  • Print its eigenvalues at the beginning of the output.

Bear in mind that this is probably related to #65.

EDIT: I think we can consider two other things:

• the relative tolerance (rtol, a good starting point is 1e-6)
• the absolute tolerance (atol, a good starting point is 1e-11)

EDIT: Other things to consider:

• Ensure the diffusion limit is enforced in all cases for reactions in solution (Collins-Kimball theory). This gives an upper bound to the reaction rates, which is the reason this is considered here.

Reactions and equilibria having several undistinguishable compounds on either side should have energy quantities adjusted by a so-called translational symmetry number (Fernández-Ramos (2007)).

As such, a reaction such as n A -> B should have an added n! symmetry number, which accounts for A molecules change position. This is important for activation free energies of bimolecular reactions.

Observe that the above describes a correction to free energies associated with the reaction, not any compound in particular.

There are many different models out there. Computational chemistry packages usually deal with the ideal-gas limit only, but there are models for hindered rotors and models for the harmonic limit and crystals.

I mostly want to predict the experimental enthalpies and entropies of molecules in the gas phase (this issue) and solution (see #26). As such, quasi-harmonic harmonic approximations are useful. Excellent models for supramolecular chemistry (mean absolute deviation of 2.1 kcal/mol) are given by Grimme (2012), together with a quasi-harmonic approximation to entropies. Li (2015) introduced a robust and theoretically sound quasi-harmonic estimate to enthalpies, with 1.8 kcal/mol of error for a particularly well-modeled system. Ribeiro (2011) offers some interesting numerical examples to test.

As suggested by G. Caramori, it would be nice to have some different units being shown at the same time for reaction rate constants when doing $ overreact model.k. For this to work, we need to:

• Show a single column in s-1 for unimolecular reactions.
• Show a couple of columns for bimolecular reactions with other units (at least M s-1, cm3 particle-1 s-1, and atm s-1).
• Generalize the above and ensure similar output for termolecular reactions.

As a low priority feature, after implementing #5, it is interesting to implement the contributions to heat capacity.

• Check dependabot & GitHub Actions

Adding some references here: https://doi.org/10.5194/acp-21-18351-2021

*I'm able to provide help, if wanted

We could benefit from a high-level interface, defined as an interface that deals with logfile and model objects directly. All functions would be higher-level versions of internal procedures. This is related to the discussion in #28.

• As such, the user would be able to do something similar to:

from overreact import api, io, simulate
import matplotlib.pyplot as plt

model = io.parse_model("my-model.jk")
for temperature in [200.0, 300.0, 400.0]:
    enthalpies = api.get_enthalpies(model.compounds, temperature=temperature)
    entropies = api.get_entropies(model.compounds, temperature=temperature)
    freeenenergies = enthalpies - temperature * entropies

    delta_freeenergies = api.get_delta(model.scheme.B, freeenergies)
    k = api.get_k(delta_freeenergies)
    kappa = api.get_kappa(delta_freeenergies, model.scheme)
    k = kappa * k

    y0 = [1.0, 0.0]
    dydt = simulate.get_dydt(scheme, k)
    t, y = simulate.get_y(dydt, y0=y0, method="Radau")
    plt.plot(t, y[1], label=f"P @{temperature}K")

plt.legend()
plt.xlabel("Time (s)")
plt.ylabel("Concentration (M)")
plt.show()

Furthermore, other related things are proposed:

• Calculate the degree of rate control with a factory function that receives the free energies as first parameter and returns one reaction rate:

...
def factory(freeenergies, scheme, y0, t):
    delta_freeenergies = api.get_delta(scheme.B, freeenergies)
    k = api.get_k(delta_freeenergies)
    kappa = api.get_kappa(delta_freeenergies, scheme)
    k = kappa * k

    dydt = simulate.get_dydt(scheme, k)
    t, y = simulate.get_y(dydt, y0, [0.0, t], method="Radau")
    return dydt(t[-1], y.T)

t = 60.0
y0 = [1.0, 0.0]
drc = get_drc(factory, args=(model.scheme, y0, t))

# degree of rate control of the formation of 1 (P) with respect to all compounds
print(drc[1, :])

The (square) array-like object drc would contain the derivatives of equation 1 of Campbell (2020) for the rate of formation of all compounds with respect to the free energies of all compounds.

I need to finish some steps in using the GNU General Public License:

• Add a COPYING LICENSE file with a copy of the GNU GPL MIT.
• Put a license notice in each file (no need for MIT).
• Give each file the proper copyright notices. Make sure to clearly identify which versions of the license users can use (no need for MIT).
• Make the program display a startup notice.

Support for type hints is a plus to the future.

Hey @schneiderfelipe, the build action for Python 3.10 still failing on pytest.

from . import vq, hierarchy
../../../.cache/pypoetry/virtualenvs/overreact-ykHAvF72-py3.10/lib/python3.10/site-packages/scipy/cluster/hierarchy.py:135: in <module>
    from . import _hierarchy, _optimal_leaf_ordering
E   ImportError: /home/runner/.cache/pypoetry/virtualenvs/overreact-ykHAvF72-py3.10/lib/python3.10/site-packages/scipy/cluster/_optimal_leaf_ordering.cpython-310-x86_64-linux-gnu.so: undefined symbol: _PyGen_Send
=============================== warnings summary ===============================
../../../.cache/pypoetry/virtualenvs/overreact-ykHAvF72-py3.10/lib/python3.10/site-packages/coverage/inorout.py:473
  /home/runner/.cache/pypoetry/virtualenvs/overreact-ykHAvF72-py3.10/lib/python3.10/site-packages/coverage/inorout.py:473: CoverageWarning: --include is ignored because --source is set (include-ignored)
    self.warn("--include is ignored because --source is set", slug="include-ignored")

We currently access logfile dictionary objects through the ["..."], but it would be friendly (and compatible with cclib) if we could use the dot syntax as well, e.g., data["energy"] being equivalent to data.energy.

The whole thing.

Semantic versioning is something we have to inform in the documentation (and README.rst).

I believe Rich can contribute a lot to the current verbose logging system.

• Overreact.jl
• Graphical User Interface (Reactive.jl) → Maybe I should implement it online as an API

@schneiderfelipe, please assign me on this issue

Deadline: 25/12

*don't worry about the deadline, @schneiderfelipe, 😂😂 it is only for me. You will be able to check it out after of your summer break.

https://doi.org/10.1021/jacs.1c07351

I currently have no configuration in the repo. This could be an issue when going public.

The user experience is essential for the widespread of the project. As such, considerable time should be devoted to its planning, and proper documentation should match an excellent interface. Currently, I envision three separate processes for using overreact in a project:

1. Create a model
2. Validate
3. Analyze

All steps take place after calculations have been effectuated. In reality, those steps can happen in a nonlinear fashion, e.g., you might calculate structures, then attempt a model and validate it, only to realize that the model needs adjustments, in which case you go back to your calculations. The analysis is the last step, where conclusions are taken, and quantitative results are gathered.

Creating a model

Model creation consists of specifying reactions and species, together with data sources from computational chemistry logfiles. A dictionary object encodes all information necessary for a model. This object can be stored in the JSON format for convenience (extension .jk). Additionally, we offer an input format for compiling JSON-encoded models files (extension .k).

Below are the usual ways of constructing a model file:

1. Using orxn -c to compile a .jk file from a .k file (might read logfiles).
2. Specifying reactions and species in Python and writing the resulting object to a .jk file (might read logfiles).

Below are the usual ways of constructing dictionary model objects:

1. Reading a .jk file to build the resulting object.
2. Reading a .k file directly (might read logfiles) to build the resulting object.
3. Specifying reactions and species in Python to build the resulting object (might read logfiles).

The most efficient way is to compile a .jk file from a .k file, which ensures logfiles are only read once. You can always come back and improve your .k file and recompile.

Validating

This is a broad concept but consists of checking how good is the model in predicting something you're interested in. As such, it consists of property calculation and comparison with experimental evidence. Some commonly compared properties are:

• Non-simulated properties:
  • Rate constants
  • Equilibrium constants
  • Arrhenius plots
  • Apparent activation energy
  • Boltzmann-averaged temperature-resolved spectroscopic properties
  • Kinetic isotope effects
• End-of-simulation properties:
  • Reaction yields
  • Product selectivity
• Simulated properties:
  • Boltzmann-averaged time-resolved spectroscopic properties
  • Reaction progress kinetic analysis (see below)

Reaction progress kinetic analysis (Blackmond (2005)) consists of a framework to investigate reaction mechanisms, in which graphs are used to understand chemical reactions:

• Reaction concentration profiles
  • Substrate concentration versus time ("integral method")
• Reaction rate profiles
  • Reaction rate versus time ("differential method")
• "Graphical rate equations" (see here as well)
  • Zeroth order
  • First-order
  • Second-order
• Turnover frequency
  • Turnover frequency versus substrate concentration
• Michaelis-Menten plots
  • Rate versus substrate concentration
  • Rate over one substrate concentration versus another substrate concentration
• Lineweaver-Burk plots
  • Reciprocal rate versus reciprocal substrate

The only things overreact needs to expose in order to allow reaction progress kinetic analysis are some properties over time (concentrations, rates, and turnover numbers and frequencies). The graphing part can be very easily done in matplotlib.

Analyzing models

In my view, the analysis consists of explaining the phenomena. This means determining the most important effects and measuring them. Determining the most important reaction steps is part of this:

• Rate-determining state analysis (Kozuch (2011))
• Degree of rate control (Campbell (2017))

The functionality of overreact ends here. After the most important steps are known, specific analysis with particular transition states can be employed somewhere else, such as the reaction force analysis (Yepes (2013)), resonance NBOs (Glendening (2019)) and the activation strain model (Vermeeren (2020)).

Support a direct bias adjustment of computed free energies for fitting a target (e.g., experimental) concentration profile (see, e.g., Pérez-Soto (2020)).

There is extensive evidence suggesting that computational errors are intrinsically systematic (EDIT: in most cases). As such, correcting energies through a learned bias (i.e., add optimum value to all calculated energies) seems justified. Furthermore, if different computational methods were used in different reaction steps, it seems logical to have a single fitted bias for each of them (EDIT: this is of very low priority at this point, so that I will skip it completely; that would require rewriting a substantial amount of code).

That allows one to correct a single kinetic model for experimental evidence and propagate the correction for a wide range of related kinetic models (e.g., similar catalysts, similar substrates, etc.).

Furthermore, this may improve posterior results of the degree of rate control (DRC) analysis. In particular, a DRC-based sensitivity analysis (see, e.g., Meskine (2009)) can be used to check which reaction steps require more accurate levels of theory.

Fixing the concentration of a particular compound can be useful for, e.g., simulating under specific pH conditions or setting pseudo order conditions (w.r.t. solvent, for instance). For the case of pH-dependent reactions, that is also necessary for solvent kinetic isotope effects (SKIE).

Below I sketch a simple method for simulating a fixed concentration of compound A:

• Calculate all reaction rate constants.
• For each direct reaction in our the model where A takes part as a reactant:
  • Multiply the reaction rate constant by the fixed concentration of A (at the power of its stoichiometric coefficient) and store it.
  • Remove A entirely from the list of reactants.

The above will generate reaction rates 1. compatible with a fixed concentration of A and 2. where A is absent. I believe this suffices, but numerical experiments must be done to access whether this is fruitful. Observe that I still need to know the energetic contribution of A to all reaction rate constants.

The method proposed here should produce the same results as, e.g., Scorsin (2020). See Sicilio (1961) for a systematic overview of errors related to the experimental treatment of pseudo first order reactions.

It is a common use case, most of the concentrations in a complex reaction network will be zero. So assume them zero if some are missing.

Often we want to sum concentrations among different conformational structures or solvated ones. Instead of requiring post-processing, one would only need to add some unique IDs shared among some species in the input.

Something like

using piccolomo/plotext.

Rich is a "Python library for rich text and beautiful formatting in the terminal." In particular, it can also render pretty tables, which will 1. make overreact look awesome, and 2. simplify a lot the current codebase regarding the interface.

❯ overreact model.k
Traceback (most recent call last):
  File "/home/schneider/.local/bin/overreact", line 5, in <module>
    from overreact._cli import main
  File "/home/schneider/.local/lib/python3.9/site-packages/overreact/_cli.py", line 15, in <module>
    import matplotlib.pyplot as plt
ModuleNotFoundError: No module named 'matplotlib'

The offending line is the following:

The following is a brief, cascading outline of all contributions to thermochemistry:

• G = H - T * S
  • H = U + kB * T
    • U = E(el) + E(ZPE) + E(vib) + E(rot) + E(trans)
  • S = S(el) + S(vib) + S(rot) + S(trans)

Above I use the following definitions:

• E(el) is the total electronic energy [equals E(kin-el) + E(nuc-el) + E(el-el) + E(nuc-nuc)]
• S(el) is the electronic entropy
• E(ZPE) is the zero temperature vibrational energy
• E(vib) is the finite temperature correction to E(ZPE) due to populated excited vibrational states
• S(vib) is the vibrational entropy
• E(rot) is the rotational thermal energy
• S(rot) is the rotational entropy
• E(trans) is the translational thermal energy
• S(trans) is the translational entropy

All values incorporated into G, H, and S will be calculated using the following models:

• RRHO
• QRRHO

For some information on how RRHO is implemented in some codes, see here, here and here. All this should be organized in a proper module. Ideas for names: thermo, thermostat, mechstat, rovib, rovibs.

Collateral features

This scheme will allow us to

• calculate all values above for any temperature
• calculate all values above for any isotope substitution

Initially, I wanted to exclusively use cclib for reading logfiles. Unfortunately, this was not possible and I patched the cclib output when reading ORCA logfiles. As such, cclib might become an optional dependency for reading logfiles other than ORCA generated ones.

When doing this, I might need to test some logfiles from at least one open source package (NWChem or Psi4 are fine).

What I propose is something similar to what I did to the data directory.

Furthermore, I want to simplify the documentation:

• Use less text, straight to the point
• Document API for only classes and functions I intend to support (choose wisely!) (EDIT: this will live in the code only for now)
• Transform all code tutorials into Jupyter Notebooks
• Transform all code examples into Jupyter Notebooks
• Test Jupyter Notebooks by using a py.test plugin such as nbval (EDIT: postponed)

People should be able to do the following as soon as they hit the documentation page:

• Learn how overreact can be useful to them
• See a code example (~5 lines of code max.)
• Learn how to install overreact
• Jump straight to tutorials in Jupyter Notebooks (GitHub-hosted subfolder)
• Jump straight to the paper (coming soon!)
• Jump straight to the license

Imagine what happens when we get

$scheme
 [...]
 // Harmless, huh?
 B -> E# -> H -> L# -> M:RNH2
 [...]
$compounds
 [...]
 // Actually eats your brain out 🤯
 M:RNH2: 0w/M/product_0.out
 [...]

Pretty bad.

Currently, overreact -h is not that useful. We need

• Properly describe all options, including their units if applicable
• Tell about the program and guide the user through online documentation, citation, license, etc.

That last point can be set into a single description string and wired so that both the output header and the help page show the same preset information.

Option: run the GUI online (?) → Add references

Related to #83. From https://sphinx-themes.org/#themes:

• Alabaster: too boring
• Read the Docs: too common
• Furo
• Book
• PyData: too informal
• Press: too corporatey
• Celery: boring²
• Insegel: corporatey²
• Insipid: ...insipid
• Material: too Androidy/Googley
• Material: ugly³
• Basicstrap: boring
• Bernard: makes me puke
• Better: ugly, boring, oh boy
• Bootstrap: makes me cry
• Cloud: God, no
• Documatt: too universitey
• Groundwork: good, but dark and too ruby/rails/jekyll
• Hachibee: too kawaii
• Karma (decent)
• Kotti: too flat
• Maisie: nice but spacing is bad
• Murray: too latexy/waxy/scholarly
• PD: too facebookey
• PDJ: too serious
• Python: too... Python!
• Readable: better Alabaster, but still bad
• Redactor: simply bad
• renku: doesn't even compile
• Sandstone: too Firefoxey
• Sizzle: too corporatey
• Solar: please no
• Stanford (nice, an improved RTDs)
• Topos: almost cried
• Yummy: we don't sell cosmetics
• ZeroVM: too 1999
• tibas.tt: too animey
• agogo: NO
• bizstyle: NO²
• classic: ugy, 1996!
• haiku: clean but boring
• nature: can't stand
• pyramid: sad gray
• scrolls: God, woody background! NO
• sphinxdoc: default is too default
• traditional: surprisingly works, but boring

Round two

Veredict

• Furo* got 2 + 0 + 1 + 1 + 1 + 2 = 7 points
• Book got 3 + 1 + 2 + 1 + 1 + 1 = 9 points
• Karma got 0 + 0 + 0 + 0 + 0 + 0 = 0 points
• Stanford got 1 + 1 + 2 + 1 + 1 + 1 = 7 points

Being a dictator I opted for Stanford because it is cleaner, has more contrast and less clutter.

This is essential for quickly checking results and is useful when working with equilibria, specificity, and selectivity.

Some links will be broken when we publish the webpage since they reflect the repository state before transferring to @geem-lab.

Correct standard states (dG°(conc)) are crucial for thermochemistry and chemical kinetics. It has been emphasized to be particularly important in the context of the degree of rate control analysis, too (see, e.g., Campbell (2017)).

• Inform about those corrections in the output (requires -vv)

Tests should, ideally, include experimental comparisons with:

• Bimolecular reactions in both gas phase and solution.
• Gas-liquid transfer processes.

Using -vv or higher triggers weird logging from matplotlib regarding fonts:

...
    findfont: score(<Font 'cmr10' (cmr10.ttf) normal normal 400
        normal>) = 10.05
    findfont: score(<Font 'Trebuchet MS' (Trebuchet_MS_Bold.ttf)
        normal normal 700 normal>) = 10.335
    findfont: score(<Font 'Latin Modern Sans Quotation'
        (lmsansquot8-oblique.otf) oblique normal 700 normal>) = 11.335
    findfont: score(<Font 'Nakula' (nakula.ttf) normal normal 400
        normal>) = 10.05
    findfont: score(<Font 'KacstDecorative' (KacstDecorative.ttf)
        normal normal 400 normal>) = 10.05
    findfont: score(<Font 'Latin Modern Mono Prop Light'
        (lmmonoproplt10-regular.otf) normal normal light normal>) =
        10.24
    findfont: score(<Font 'TeX Gyre Termes' (texgyretermes-italic.otf)
        italic normal 400 normal>) = 11.05
    findfont: Matching :family=sans-serif:style=normal:variant=normal:
        weight=normal:stretch=normal:size=10.0 to DejaVu Sans
        ('/usr/share/matplotlib/mpl-data/fonts/ttf/DejaVuSans.ttf')
        with score of 0.050000.

This should be suppressed.

As suggested by G. Caramori, it is better if we inform the user about units in all logging messages, as units are different from the normal output. As such,

$ overreact model.k -v

will properly show units in all instances.

In the acetic acid example, I get the equivalent of ~0.3 kcal/mol of error when using final concentrations or direct/reverse reaction rate constants for the calculation of the acid dissociation constant, when compared to using the Gibbs free energies.

This is a very small error, but it might be worth checking where it is coming from.

Say,

$scheme
 A(w) -> B(w)
 [...]

and then

$ overreact model.k "A:1"  # should have been "A(w):1"

Recipe for trouble 🙅🏽

As pointed out by G. Caramori, we log two different values for rotational internal energies and entropies:

@ calculate internal energy: NH3
    translational energy = 3718.43554433505
    electronic energy = -148242230.4843498
    rotational energy = 3667.0130846884526
    rotational energy = 0.05436732296723629
    vibrational energy = 87465.47422773033
    internal energy = -148147379.56149304

...

@ calculate entropy: NH3
    point group = C3v
    symmetry number = 3
    translational entropy = 144.1035562838303
    electronic entropy = 0.0
    rotational entropy = 48.087550244988776
    rotational entropy = 0.0002310639595995566
    vibrational entropy = 0.4263105601512432
    entropy = 192.61741708897034

That has to do with how we calculate the QRRHO model (the smaller value is the quasi-harmonic component of the vibrational energies and entropies). It would be interesting to either remove it from the log or to state its meaning explicitly.

This has to do with #36, too, as the units have to be shown as well.