Allow the user to choose a simulation time

docsrc/notebooks/#2 A mini language for chemical reactions.ipynb

@@ -151,7 +151,7 @@
    "source": [
     "# TODO: scipy 1.4+ will support passing args from solve_ivp, but for now let's do something rather mundane here\n",
     "def dydt(t, y):\n",
-    "    A = np.asanyarray(scheme.A)\n",
+    "    A = np.asarray(scheme.A)\n",
     "    r = k * np.prod(np.power(y, np.where(A > 0, 0, -A).T), axis=1)\n",
     "    return np.dot(scheme.A, r)"
    ]
@@ -216,7 +216,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.8.6-final"
+   "version": "3.8.5-final"
   }
  },
  "nbformat": 4,

docsrc/tutorials/2-equilibrium-constants.rst

@@ -20,7 +20,7 @@ So let's check it:
 >>> scheme.compounds, scheme.reactions
 (('Cd2p', 'MeNH2', '[Cd(MeNH2)4]2p'),
  ('Cd2p + 4 MeNH2 -> [Cd(MeNH2)4]2p', '[Cd(MeNH2)4]2p -> Cd2p + 4 MeNH2'))
->>> dydt = simulate.get_dydt(scheme, [K, 1.])
+>>> dydt = simulate.get_dydt(scheme, np.array([K, 1.]))
 >>> y, r = simulate.get_y(dydt, y0=[0., 0., 1.])
 >>> y(y.t_max)
 array([0.01608807, 0.06435228, 0.98391193])

docsrc/tutorials/simulation/0-basic-reaction-simulation.rst

@@ -42,8 +42,9 @@ The above model translates to the following in overreact:
 overreact helps us to define the equations and solve the initial value problem.
 First, let's define the system of ordinary differential equations:
 
+>>> import numpy as np
 >>> from overreact import simulate
->>> dydt = simulate.get_dydt(scheme, [kf])
+>>> dydt = simulate.get_dydt(scheme, np.array([kf]))
 
 The returned object above is a function of concentrations and time that defines
 a set of ordinary differential equations in time:
@@ -57,8 +58,6 @@ We are going to simulate 10 seconds, starting with an initial concentration of
 reaction rate constant).
 
 >>> y, r = simulate.get_y(dydt, y0=[1., 0.])
-
->>> import numpy as np
 >>> t = np.linspace(y.t_min, 5.0)
 
 The simulation data is stored in `t` (points in time) and `y` (concentrations).
@@ -84,4 +83,4 @@ We can see that the reaction went to full completion by checking the final
 concentrations:
 
 >>> y(y.t_max)
-array([0.0000, 1.0000])
+array([0.00, 1.00])

docsrc/tutorials/simulation/1-basic-reaction-simulation.rst

@@ -72,8 +72,9 @@ values from the table:
 We need one reaction rate constant per reaction, in a
 list with the same ordering as the ``scheme.reactions`` variable obtained above.
 
+>>> import numpy as np
 >>> k_f = (k_r + k_cat) / K_M  # using the equation above
->>> k = [k_f, k_r, k_cat]
+>>> k = np.array([k_f, k_r, k_cat])
 
 Solving the initial value problem
 ---------------------------------
@@ -83,7 +84,6 @@ units, as long as they match with the reaction rate constants) and solve the
 initial value problem. Below we set the substrate to one molar and the enzyme
 to 5% of it:
 
->>> import numpy as np
 >>> y0 = np.array([0.05, 1.00, 0.00, 0.00, 0.00])
 
 One return value that ``core.parse_reactions`` has given us was the :math:`A`

overreact/_thermo/__init__.py

@@ -510,7 +510,7 @@ def get_molecularity(transform):
     >>> get_molecularity([[0.], [0.]])
     array([1])
     """
-    res = np.sum(np.asanyarray(transform) < 0, axis=0)
+    res = np.sum(np.asarray(transform) < 0, axis=0)
     return np.where(res > 0, res, 1)
 
 
@@ -554,7 +554,7 @@ def get_delta(transform, property):
     ...            [ 1,  3]], [-5, 12])
     array([17, 46])
     """
-    return np.asanyarray(transform).T @ np.asanyarray(property)
+    return np.asarray(transform).T @ np.asarray(property)
 
 
 # TODO(schneiderfelipe): further test the usage of delta_moles with values
@@ -614,13 +614,13 @@ def equilibrium_constant(
     >>> equilibrium_constant(64187.263215698644, delta_moles=-2, temperature=745.0)
     0.118e-6
     """
-    temperature = np.asanyarray(temperature)
+    temperature = np.asarray(temperature)
 
     if volume is None:
         volume = molar_volume(temperature=temperature, pressure=pressure)
 
     equilibrium_constant = (
-        np.exp(-np.asanyarray(delta_freeenergy) / (constants.R * temperature))
+        np.exp(-np.asarray(delta_freeenergy) / (constants.R * temperature))
         * (volume) ** -delta_moles
     )
     logger.info(f"equilibrium constant = {equilibrium_constant}")
@@ -704,7 +704,7 @@ def change_reference_state(
     if old_reference is None:
         if volume is None:
             volume = molar_volume(
-                temperature=np.asanyarray(temperature), pressure=pressure
+                temperature=np.asarray(temperature), pressure=pressure
             )
         old_reference = 1.0 / volume
     return sign * constants.R * np.log(new_reference / old_reference)
@@ -745,5 +745,5 @@ def get_reaction_entropies(transform):
     >>> get_reaction_entropies(scheme.A)
     array([-5.763, 5.763])
     """
-    sym = _factorial(np.abs(np.asanyarray(transform)))
+    sym = _factorial(np.abs(np.asarray(transform)))
     return np.sum(np.sign(transform) * change_reference_state(sym, 1), axis=0)

overreact/_thermo/_gas.py

@@ -95,7 +95,7 @@ def calc_elec_energy(energy=0.0, degeneracy=1, temperature=298.15):
     321.00
 
     """
-    min_energy = np.asanyarray(energy).min()
+    min_energy = np.asarray(energy).min()
     if np.isclose(temperature, 0.0):
         logger.warning("assuming ground state as electronic energy at zero temperature")
         return min_energy
@@ -176,7 +176,7 @@ def calc_elec_entropy(energy=0.0, degeneracy=1, temperature=298.15):
         logger.warning("assuming electronic entropy zero at zero temperature")
         return 0.0
 
-    min_energy = np.asanyarray(energy).min()
+    min_energy = np.asarray(energy).min()
     energy = energy - min_energy
 
     q_elec_terms = degeneracy * np.exp(-energy / (constants.R * temperature))
@@ -880,6 +880,6 @@ def molar_volume(temperature=298.15, pressure=constants.atm):
     >>> molar_volume() / (constants.angstrom ** 3 * constants.N_A)
     40625.758632362515
     """
-    molar_volume = constants.R * np.asanyarray(temperature) / np.asanyarray(pressure)
+    molar_volume = constants.R * np.asarray(temperature) / np.asarray(pressure)
     logger.debug(f"molar volume = {molar_volume} ų")
     return molar_volume

overreact/api.py

@@ -19,15 +19,13 @@
 from overreact.core import get_transition_states
 from overreact.core import is_transition_state
 from overreact.core import parse_reactions
-from overreact.datasets import data_path
 from overreact.io import parse_compounds
 from overreact.io import parse_model
 from overreact.simulate import get_dydt
 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
-from overreact.misc import cache
 
 logger = logging.getLogger(__name__)
 
@@ -264,10 +262,9 @@ def get_freeenergies(
     enthalpies = get_enthalpies(compounds, qrrho=qrrho, temperature=temperature)
     entropies = get_entropies(compounds, qrrho=qrrho, temperature=temperature)
     # TODO(schneiderfelipe): log the contribution of bias
-    return enthalpies - temperature * entropies + _np.asanyarray(bias)
+    return enthalpies - temperature * entropies + _np.asarray(bias)
 
 
-# @cache
 def get_k(
     scheme,
     compounds=None,

overreact/cli.py

@@ -74,6 +74,7 @@ def __init__(
         temperature=298.15,
         bias=0.0,
         method="Radau",
+        max_time=5 * 60,
         rtol=1e-5,
         atol=1e-11,
         box_style=box.SIMPLE,
@@ -86,6 +87,7 @@ def __init__(
         self.temperature = temperature
         self.bias = bias
         self.method = method
+        self.max_time = max_time
         self.rtol = rtol
         self.atol = atol
         self.box_style = box_style
@@ -474,7 +476,14 @@ def _yield_kinetics(self):
                 # TODO(schneiderfelipe): the following is inefficient but probably OK
                 y0[self.model.scheme.compounds.index(name)] = quantity
 
-            y, r = api.get_y(dydt, y0=y0)
+            y, r = api.get_y(
+                dydt,
+                y0=y0,
+                method=self.method,
+                rtol=self.rtol,
+                atol=self.atol,
+                max_time=self.max_time,
+            )
             conc_table = Table(
                 Column("no", justify="right"),
                 Column("compound", justify="left"),
@@ -492,9 +501,13 @@ def _yield_kinetics(self):
                 )
             yield conc_table
 
-            # TODO(schneiderfelipe): we can get a max time now based on the
-            # changes through time: stop when the graph gets boring.
-            t = np.linspace(y.t_min, y.t_max)
+            t_max = y.t_max
+            factor = y(y.t_max).max()
+            reference = y(y.t_max) / factor
+            while np.allclose(y(t_max) / factor, reference, atol=1e-2):
+                t_max = 0.95 * t_max
+
+            t = np.linspace(y.t_min, t_max, num=100)
             if self.plot:
                 for i, name in enumerate(self.model.scheme.compounds):
                     if not core.is_transition_state(name):
@@ -554,6 +567,7 @@ def main():
         choices=["BDF", "LSODA", "Radau"],
         default="Radau",
     )
+    parser.add_argument("--max-time", type=float, default=5 * 60)
     parser.add_argument("--rtol", type=float, default=1e-5)
     parser.add_argument("--atol", type=float, default=1e-11)
     parser.add_argument(
@@ -596,6 +610,7 @@ def main():
 - Temperature    = {args.temperature} K
 - Pressure       = {args.pressure} Pa
 - Integrator     = {args.method}
+- Max. Time      = {args.max_time}
 - Rel. Tol.      = {args.rtol}
 - Abs. Tol.      = {args.atol}
 - Bias           = {args.bias / constants.kcal} kcal/mol
@@ -622,6 +637,7 @@ def main():
         temperature=args.temperature,
         bias=args.bias,
         method=args.method,
+        max_time=args.max_time,
         rtol=args.rtol,
         atol=args.atol,
     )

overreact/coords.py

@@ -728,8 +728,8 @@ def _get_proper_axes(
     if rotor_class[1] == "atomic" or len(atomcoords) == 1:
         return list()
 
-    axes = np.asanyarray(axes)
-    atomcoords = np.asanyarray(atomcoords)
+    axes = np.asarray(axes)
+    atomcoords = np.asarray(atomcoords)
     orders = _guess_orders(groups, rotor_class)
 
     found_axes = list()
@@ -944,8 +944,8 @@ def _get_improper_axes(
     if rotor_class[1] == "atomic" or len(atomcoords) == 1:
         return list()
 
-    axes = np.asanyarray(axes)
-    atomcoords = np.asanyarray(atomcoords)
+    axes = np.asarray(axes)
+    atomcoords = np.asarray(atomcoords)
 
     if proper_axes is None:
         proper_axes = _get_proper_axes(
@@ -1089,8 +1089,8 @@ def _get_mirror_planes(
     if rotor_class[1] == "atomic" or len(atomcoords) == 1:
         return list()
 
-    axes = np.asanyarray(axes)
-    atomcoords = np.asanyarray(atomcoords)
+    axes = np.asarray(axes)
+    atomcoords = np.asarray(atomcoords)
 
     if proper_axes is None:
         proper_axes = _get_proper_axes(
@@ -1192,7 +1192,7 @@ def _has_inversion_center(atomcoords, groups, rtol=0.0, atol=1.0e-2, slack=1.888
     """
     rtol, atol = slack * rtol, slack * atol
 
-    atomcoords = np.asanyarray(atomcoords)
+    atomcoords = np.asarray(atomcoords)
     return all(
         _is_symmetric(atomcoords[group], _operation("i"), rtol=rtol, atol=atol)
         for group in groups[::-1]
@@ -1332,7 +1332,7 @@ def _operation(name, order=2, axis=None):
     elif name == "e":
         return np.eye(3)
     elif name in {"c", "σ", "sigma", "s"}:  # normalize axis
-        axis = np.asanyarray(axis)
+        axis = np.asarray(axis)
         axis = axis / np.linalg.norm(axis)
 
         if name in {"c", "s"}:
@@ -1424,7 +1424,7 @@ def _classify_rotor(moments, rtol=0.0, atol=1.0e-2, slack=0.870):
 
     if np.isclose(moments[2], 0.0, rtol=inner_slack * rtol, atol=inner_slack * atol):
         return "spheric", "atomic"
-    moments = np.asanyarray(moments) / moments[2]
+    moments = np.asarray(moments) / moments[2]
 
     # basic tests for tops
     is_oblate = np.isclose(
@@ -1645,7 +1645,7 @@ def calc_hessian(atommasses, atomcoords, vibfreqs, vibdisps):
             -0.01313383, -0.00678954,  0.25045664,  0.29443748]])
     """
     dof = 3 * len(atommasses)
-    L_cart = np.asanyarray(vibdisps).reshape((len(vibfreqs), dof)).T
+    L_cart = np.asarray(vibdisps).reshape((len(vibfreqs), dof)).T
     # correct until here
     L_cart = np.linalg.qr(L_cart, mode="complete")[0]
 
@@ -1657,7 +1657,7 @@ def calc_hessian(atommasses, atomcoords, vibfreqs, vibdisps):
     assert np.allclose(M @ D @ L, L_cart)
 
     # correct from here
-    nu = np.asanyarray(vibfreqs) * constants.c / constants.centi
+    nu = np.asarray(vibfreqs) * constants.c / constants.centi
     eigenvalues = (
         (2.0 * np.pi * nu) ** 2
         * (constants.atomic_mass * constants.bohr ** 2)
@@ -1701,7 +1701,7 @@ def calc_vibfreqs(hessian, atommasses):
     atommasses_sqrt = np.sqrt([mass for mass in atommasses for _ in range(3)])
 
     # mass-weighted Hessian
-    hessian = np.asanyarray(hessian) / np.outer(atommasses_sqrt, atommasses_sqrt)
+    hessian = np.asarray(hessian) / np.outer(atommasses_sqrt, atommasses_sqrt)
 
     eigenvalues = np.linalg.eigvals(hessian)
     # TODO(schneiderfelipe): the following probably misses some linear
@@ -1772,7 +1772,7 @@ def eckart_transform(atommasses, atomcoords):
              2.94635849e-01, -3.12654446e-01,  5.71869440e-01,
              5.73721626e-01, -1.13019078e-01,  3.00863871e-01]])
     """
-    atommasses = np.asanyarray(atommasses)
+    atommasses = np.asarray(atommasses)
     natom = len(atommasses)
     dof = 3 * natom
 

overreact/core.py

@@ -111,7 +111,7 @@ def get_transition_states(A, B, is_half_equilibrium):
     (None, None, 3)
 
     """
-    tau = np.asanyarray(B) - np.asanyarray(A) > 0  # transition state matrix
+    tau = np.asarray(B) - np.asarray(A) > 0  # transition state matrix
     return tuple(
         x if not is_half_equilibrium[i] and tau[:, i].any() else None
         for i, x in enumerate(np.argmax(tau, axis=0))

overreact/rates.py

@@ -100,17 +100,17 @@ def smoluchowski(
     >>> smoluchowski(radii, viscosity=8.91e-4) / constants.liter  # doctest: +SKIP
     3.6e9
     """
-    radii = np.asanyarray(radii)
-    temperature = np.asanyarray(temperature)
-    pressure = np.asanyarray(pressure)  # TODO(schneiderfelipe): do we need this?
+    radii = np.asarray(radii)
+    temperature = np.asarray(temperature)
+    pressure = np.asarray(pressure)  # TODO(schneiderfelipe): do we need this?
 
     if mutual_diff_coef is None:
         if callable(viscosity):
             viscosity = viscosity(temperature)
         elif isinstance(viscosity, str):
             viscosity = liquid_viscosity(viscosity, temperature, pressure)
         mutual_diff_coef = (
-            constants.k * temperature / (6.0 * np.pi * np.asanyarray(viscosity))
+            constants.k * temperature / (6.0 * np.pi * np.asarray(viscosity))
         ) * np.sum(1.0 / radii)
 
     if reactive_radius is None:
@@ -309,8 +309,8 @@ def eyring(
     >>> eyring(18.86 * constants.kcal)  # unimolecular, s-1
     0.093
     """
-    temperature = np.asanyarray(temperature)
-    delta_freeenergy = np.asanyarray(delta_freeenergy)
+    temperature = np.asarray(temperature)
+    delta_freeenergy = np.asarray(delta_freeenergy)
     delta_moles = 1 - molecularity
     return (
         _thermo.equilibrium_constant(