Ensure standard state corrections are correct (#64)

geem-lab · Nov 25, 2020 · d568fb9 · d568fb9
1 parent 7ba6508
commit d568fb9
Show file tree

Hide file tree

Showing 18 changed files with 122 additions and 43 deletions.
diff --git a/data b/data
diff --git a/docs/_static/drc.png b/docs/_static/drc.png
diff --git a/docs/_static/first-order.png b/docs/_static/first-order.png
diff --git a/docs/_static/michaelis-menten-dydt.png b/docs/_static/michaelis-menten-dydt.png
diff --git a/docs/_static/michaelis-menten-tof.png b/docs/_static/michaelis-menten-tof.png
diff --git a/docs/_static/michaelis-menten.png b/docs/_static/michaelis-menten.png
diff --git a/docs/_static/simple-first-order.png b/docs/_static/simple-first-order.png
diff --git a/docsrc/notebooks/#2 A mini language for chemical reactions.ipynb b/docsrc/notebooks/#2 A mini language for chemical reactions.ipynb
@@ -59,6 +59,7 @@
    "metadata": {},
    "outputs": [
     {
+     "output_type": "execute_result",
      "data": {
       "text/plain": [
        "(['A', 'B', 'TS*', 'S', 'E', 'ES', 'ES*', 'P'],\n",
@@ -86,9 +87,8 @@
        "  [0.0, 0.0, 0.0, 0.0, 0.0, 1.0]])"
       ]
      },
-     "execution_count": 4,
      "metadata": {},
-     "output_type": "execute_result"
+     "execution_count": 4
     }
    ],
    "source": [
@@ -112,26 +112,23 @@
    "metadata": {},
    "outputs": [
     {
-     "data": {
-      "text/plain": [
-       "(array([ 8. , -5. , -5. ,  0. ,  6.5, -4.5]),\n",
-       " array([1.36760068e-06, 4.62408100e+03, 4.62408100e+03, 1.00000000e+00,\n",
-       "        1.71975641e-05, 1.98848761e+03]),\n",
-       " array([8.49613440e+06, 2.87268165e+16, 2.87268165e+16, 6.21243799e+12,\n",
-       "        1.06838801e+08, 1.23533560e+16]))"
-      ]
-     },
-     "execution_count": 5,
-     "metadata": {},
-     "output_type": "execute_result"
+     "output_type": "error",
+     "ename": "AttributeError",
+     "evalue": "module 'scipy.constants' has no attribute 'kcal'",
+     "traceback": [
+      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
+      "\u001b[0;31mAttributeError\u001b[0m                            Traceback (most recent call last)",
+      "\u001b[0;32m<ipython-input-5-adcc75c3508e>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[1;32m      3\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      4\u001b[0m \u001b[0mdelta_gibbs\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0m_thermo\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mget_delta\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mscheme\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mB\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m[\u001b[0m\u001b[0;36m0.\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;36m1.\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;36m10.\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;36m5.\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;36m1.\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;36m1.\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;36m7.5\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;36m2.\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m----> 5\u001b[0;31m \u001b[0mK\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mexp\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m-\u001b[0m \u001b[0mdelta_gibbs\u001b[0m \u001b[0;34m/\u001b[0m \u001b[0;34m(\u001b[0m\u001b[0mconstants\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mR\u001b[0m \u001b[0;34m*\u001b[0m \u001b[0mT\u001b[0m \u001b[0;34m/\u001b[0m \u001b[0mconstants\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mkcal\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m      6\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      7\u001b[0m \u001b[0;31m# This is Eyring's equation\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
+      "\u001b[0;31mAttributeError\u001b[0m: module 'scipy.constants' has no attribute 'kcal'"
+     ]
     }
    ],
    "source": [
     "T = 298.15\n",
     "kappa = 1.\n",
     "\n",
     "delta_gibbs = _thermo.get_delta(scheme.B, [0., 1., 10., 5., 1., 1., 7.5, 2.])\n",
-    "K = np.exp(- delta_gibbs / (constants.R * T / (constants.kilo * constants.calorie)))\n",
+    "K = np.exp(- delta_gibbs / (constants.R * T / constants.kcal))\n",
     "\n",
     "# This is Eyring's equation\n",
     "k = kappa * (constants.k * T / constants.h) * K\n",
@@ -219,9 +216,9 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.8.5"
+   "version": "3.8.5-final"
   }
  },
  "nbformat": 4,
  "nbformat_minor": 2
-}
+}
diff --git a/docsrc/tutorials/1-rate_constant.rst b/docsrc/tutorials/1-rate_constant.rst
@@ -5,9 +5,9 @@ There are functions for calculating reaction rate constants from energy deltas.
 But first we need an energy barrier:
 
 >>> import numpy as np
->>> from scipy.constants import kilo, calorie
+>>> from overreact import constants
 >>> # TODO(schneiderfelipe): transform all "delta" into "delta_freeenergy"
->>> delta = np.array([17.26, 18.86]) * kilo * calorie
+>>> delta = np.array([17.26, 18.86]) * constants.kcal
 >>> delta  # J/mol
 array([72215.8, 78910.2])
 

diff --git a/docsrc/tutorials/2-quantum_tunnel.rst b/docsrc/tutorials/2-quantum_tunnel.rst
@@ -5,7 +5,7 @@ Transmission coefficients (:math:`\kappa`) take quantum tunneling effects into
 account in reaction rate constants.
 
 >>> from overreact import api, tunnel
->>> from scipy import constants
+>>> from overreact import constants
 >>> tunnel.wigner(266.5144)
 1.069
 >>> kappa = tunnel.wigner(1000.0)
@@ -26,7 +26,7 @@ array([-15.,   0.,  15.])
 Now comes Eyring's equation
 
 >>> T = 298.15
->>> K = np.exp(- delta_gibbs / (constants.R * T / (constants.kilo * constants.calorie)))
+>>> K = np.exp(- delta_gibbs / (constants.R * T / constants.kcal))
 >>> k = kappa * (constants.k * T / constants.h) * K
 >>> K, k
 (array([9.887e+10, 1.000e+00, 1.011e-11]),

diff --git a/docsrc/tutorials/4-curtin-hammett.rst b/docsrc/tutorials/4-curtin-hammett.rst
@@ -3,8 +3,8 @@ Curtin-Hammett Principle
 
 Let's do a toy example that addresses the Curtin-Hammett principle:
 
->>> from scipy.constants import kilo, calorie
 >>> from overreact import core, api, rates, simulate
+>>> from overreact import constants
 >>> scheme = core.parse_reactions("""
 ...     I1 <=> I2
 ...     I1 -> T1‡ -> P1
@@ -20,7 +20,7 @@ The selectivity is defined as the ratio between both products:
 >>> scheme.compounds
 ['I1', 'I2', 'T1‡', 'P1', 'T2‡', 'P2']
 >>> freeenergy = [0.0, -0.5, 10.0, 1.0, 11.0, 1.0]
->>> k = rates.eyring(kilo * calorie * api.get_delta(scheme.B, freeenergy))
+>>> k = rates.eyring(constants.kcal * api.get_delta(scheme.B, freeenergy))
 >>> k
 array([1.445e+13, 6.212e+12, 2.905e+05, 2.310e+04])
 

diff --git a/docsrc/tutorials/degree-of-rate-control/degree-of-rate-control.rst b/docsrc/tutorials/degree-of-rate-control/degree-of-rate-control.rst
@@ -20,8 +20,8 @@ where ``*`` denotes a free surface site.
 
 >>> import numpy as np
 >>> from scipy.misc import derivative
->>> from scipy.constants import kilo, calorie
->>> k = rates.eyring(kilo * calorie * api.get_delta(scheme.B, [0, 0, -2.5, 0, -2.5, 15.0, 0.0, -2.5]))
+>>> from overreact import constants
+>>> k = rates.eyring(constants.kcal * api.get_delta(scheme.B, [0, 0, -2.5, 0, -2.5, 15.0, 0.0, -2.5]))
 >>> k
 array([4.2245e+14, 6.2124e+12, 4.2245e+14, 6.2124e+12,
        1.3588e-02, 4.2245e+14, 6.2124e+12])

diff --git a/overreact/api.py b/overreact/api.py
@@ -26,6 +26,7 @@
 from overreact.simulate import get_y
 from overreact._thermo import get_delta
 from overreact._thermo import get_reaction_entropies
+from overreact._thermo import change_reference_state
 
 logger = logging.getLogger(__name__)
 
@@ -409,7 +410,6 @@ def get_k(
             pressure=pressure,
         )
 
-        # TODO(schneiderfelipe): test this implementation of reaction symmetry
         # TODO(schneiderfelipe): log the contribution of reaction symmetry
         delta_freeenergies = get_delta(
             scheme.B, freeenergies

diff --git a/overreact/cli.py b/overreact/cli.py
@@ -183,11 +183,11 @@ def _yield_scheme(self):
             reactants, _, products = re.split(r"\s*(->|<=>|<-)\s*", reaction)
             # TODO(schneiderfelipe): should we use "No" instead of None for
             # "half-equilib.?"?
-            row = [f"{i:d}", reactants, None, products, None]
+            row = [f"{i:d}", reactants, None, products, "NO"]
             if transition_states[i] is not None:
                 row[2] = scheme.compounds[transition_states[i]]
             elif scheme.is_half_equilibrium[i]:
-                row[4] = True
+                row[4] = "YES"
             parsed_table.add_row(*row)
         yield parsed_table
 
@@ -229,8 +229,10 @@ def _yield_compounds(self):
             box=self.box_style,
         )
         for i, (name, data) in enumerate(self.model.compounds.items()):
-            path_text = Text(data.logfile)
-            path_text.highlight_regex(r"[^\/]+$", "bright_blue")
+            path_text = None
+            if data.logfile is not None:
+                path_text = Text(data.logfile)
+                path_text.highlight_regex(r"[^\/]+$", "bright_blue")
             logfiles_table.add_row(
                 f"{i:d}",
                 name,
@@ -240,10 +242,12 @@ def _yield_compounds(self):
                 path_text,
             )
 
-            vibfreqs_text = Text(
-                ", ".join([f"{vibfreq:+7.1f}" for vibfreq in data.vibfreqs[:3]])
-            )
-            vibfreqs_text.highlight_regex(r"-\d+\.\d", "bright_yellow")
+            vibfreqs_text = None
+            if data.vibfreqs is not None:
+                vibfreqs_text = Text(
+                    ", ".join([f"{vibfreq:+7.1f}" for vibfreq in data.vibfreqs[:3]])
+                )
+                vibfreqs_text.highlight_regex(r"-\d+\.\d", "bright_yellow")
             compounds_table.add_row(
                 f"{i:d}",
                 name,
@@ -281,6 +285,16 @@ def _yield_thermochemistry(self):
             self.model.compounds, qrrho=self.qrrho, temperature=self.temperature
         )
         freeenergies = enthalpies - self.temperature * entropies
+        assert np.allclose(
+            freeenergies,
+            api.get_freeenergies(
+                self.model.compounds,
+                # bias=bias,
+                qrrho=self.qrrho,
+                temperature=self.temperature,
+                # pressure=pressure,
+            ),
+        )
 
         compounds_table = Table(
             Column("no", justify="right"),
@@ -309,25 +323,33 @@ def _yield_thermochemistry(self):
         delta_energies = api.get_delta(scheme.A, energies)
         delta_internal_energies = api.get_delta(scheme.A, internal_energies)
         delta_enthalpies = api.get_delta(scheme.A, enthalpies)
-        # TODO(schneiderfelipe): test this implementation of reaction symmetry
         # TODO(schneiderfelipe): log the contribution of reaction symmetry
         delta_entropies = api.get_delta(
             scheme.A, entropies
         ) + api.get_reaction_entropies(scheme.A)
         delta_freeenergies = delta_enthalpies - self.temperature * delta_entropies
+        assert np.allclose(
+            delta_freeenergies,
+            api.get_delta(scheme.A, freeenergies)
+            - self.temperature * api.get_reaction_entropies(scheme.A),
+        )
 
         delta_activation_mass = api.get_delta(scheme.B, molecular_masses)
         delta_activation_energies = api.get_delta(scheme.B, energies)
         delta_activation_internal_energies = api.get_delta(scheme.B, internal_energies)
         delta_activation_enthalpies = api.get_delta(scheme.B, enthalpies)
-        # TODO(schneiderfelipe): test this implementation of reaction symmetry
         # TODO(schneiderfelipe): log the contribution of reaction symmetry
         delta_activation_entropies = api.get_delta(
             scheme.B, entropies
         ) + api.get_reaction_entropies(scheme.B)
         delta_activation_freeenergies = (
             delta_activation_enthalpies - self.temperature * delta_activation_entropies
         )
+        assert np.allclose(
+            delta_activation_freeenergies,
+            api.get_delta(scheme.B, freeenergies)
+            - self.temperature * api.get_reaction_entropies(scheme.B),
+        )
 
         circ_table = Table(
             Column("no", justify="right"),
@@ -450,9 +472,9 @@ def _yield_kinetics(self):
             box=self.box_style,
         )
         for i, reaction in enumerate(self.model.scheme.reactions):
-            row = [f"{i:d}", reaction, None] + [f"{k[scale][i]:.3g}" for scale in k]
+            row = [f"{i:d}", reaction, "NO"] + [f"{k[scale][i]:.3g}" for scale in k]
             if self.model.scheme.is_half_equilibrium[i]:
-                row[2] = True
+                row[2] = "YES"
 
             kinetics_table.add_row(*row)
         yield kinetics_table

diff --git a/overreact/core.py b/overreact/core.py
@@ -769,7 +769,7 @@ def _parse_side(side):
     '2 *A*1* + 40 B1 + chlorophyll'
 
     """
-    for token in re.split(r"\s*\+\s*", side):
+    for token in re.split(r"\s+\+\s+", side):
         token = re.match(
             r"\s*(?P<coefficient>\d+)?\s*(?P<compound>[^\s]+)\s*", token
         ).groupdict(1)

diff --git a/tests/test_rates.py b/tests/test_rates.py
@@ -66,19 +66,19 @@ def test_basic_example_for_chemical_kinetics():
     delta_freeenergy -= (
         temperatures
         * _thermo.change_reference_state(temperature=temperatures)
-        / (constants.kilo * constants.calorie)
+        / constants.kcal
     )  # 1 atm to 1 M
     assert delta_freeenergy == pytest.approx([6.9, 8.4, 8.4, 9.9], 8e-3)
 
     delta_freeenergy += (
         temperatures
         * _thermo.change_reference_state(4, 1, sign=-1, temperature=temperatures)
-        / (constants.kilo * constants.calorie)
+        / constants.kcal
     )  # 4-fold symmetry TS
     assert delta_freeenergy == pytest.approx([6.3, 7.6, 7.6, 8.8], 9e-3)
 
     k = rates.eyring(
-        delta_freeenergy * constants.kilo * constants.calorie,
+        delta_freeenergy * constants.kcal,
         temperature=temperatures,
         molecularity=2,
     )

diff --git a/tests/test_regressions.py b/tests/test_regressions.py
@@ -0,0 +1,52 @@
+#!/usr/bin/env python3
+
+"""Regressions against experimental/reference values."""
+
+# TODO(schneiderfelipe): transfer all comparisons with experimental/reference
+# values to this file.
+
+import numpy as np
+import pytest
+
+from overreact import api
+from overreact import constants
+
+
+def test_solubility_of_acetic_acid():
+    """Reproduce literature data for AcOH(g) <=> AcOH(aq).
+
+    Data is as cited in doi:10.1021/jp810292n and is experimental except when
+    otherwise indicated in the comments.
+    """
+    model = api.parse_model("data/acetate/model.k")
+    temperature = 298.15
+    pK = 4.76  # doi:10.1063/1.1416902
+
+    enthalpies = api.get_enthalpies(model.compounds, temperature=temperature)
+    entropies = api.get_entropies(model.compounds, temperature=temperature)
+
+    delta_enthalpies = api.get_delta(model.scheme.A, enthalpies)
+    # TODO(schneiderfelipe): log the contribution of reaction symmetry
+    delta_entropies = api.get_delta(
+        model.scheme.A, entropies
+    ) + api.get_reaction_entropies(model.scheme.A)
+    delta_freeenergies = delta_enthalpies - temperature * delta_entropies
+
+    assert delta_freeenergies / constants.kcal == pytest.approx(
+        [6.50, -6.50, -7.94, 7.94, -74.77, 74.77], 3e-4
+    )
+
+    concentration_correction = temperature * api.change_reference_state(
+        sign=-1.0,
+        temperature=temperature,
+    )
+    # doi:10.1021/jp810292n
+    assert delta_freeenergies[2] / constants.kcal == pytest.approx(
+        -6.70 + concentration_correction / constants.kcal, 8e-2
+    )
+    assert delta_freeenergies[0] == pytest.approx(
+        -constants.R * temperature * np.log(10 ** -pK), 7e-4
+    )
+
+    k = api.get_k(model.scheme, model.compounds, temperature=temperature)
+    assert -np.log10(k[0] / k[1]) == pytest.approx(pK, 5e-2)
diff --git a/tests/test_thermo_gas.py b/tests/test_thermo_gas.py
@@ -1283,3 +1283,11 @@ def test_can_calculate_reaction_entropies():
         )
         * np.array([1, -1])
     )
+
+    assert _thermo.get_reaction_entropies(
+        core.parse_reactions("2A(g) + 3B(g) + 4C(g) <=> 5D(g) + 6E(g) + 7F(g").A
+    ) == pytest.approx(
+        _thermo.get_reaction_entropies(
+            core.parse_reactions("2A(w) + 3B(w) + 4C(w) <=> 5D(w) + 6E(w) + 7F(w").A
+        )
+    )