Improve CLI, plots, and some tests (#73)

* Test more liquid viscosities

* Check and manage some TODOs and docs

* Enlarge the default max simulation time

* Base max. time and plot on active species

* Select plot points according to simulation profile

* Allow selection of tunneling correction in the CLI

* Use an integrated keyword for --plot

* Improve CLI help page

Co-authored-by: Felipe S. S. Schneider <[email protected]>

.flake8

@@ -3,7 +3,7 @@ max-line-length = 88
 max-doc-length = 88
 enable-extensions = H106,H203,H204,H205,H210,H904
 extend-ignore = E203
-max-complexity = 17
+max-complexity = 18
 doctests = True
 docstring-convention = numpy
 include-in-doctest = "*.rst"

overreact/_thermo/__init__.py

@@ -15,8 +15,6 @@
 
 logger = logging.getLogger(__name__)
 
-# TODO(schneiderfelipe): implement frequency scaling factors.
-
 
 def calc_trans_entropy(
     atommasses,
@@ -83,7 +81,7 @@ def calc_trans_entropy(
     # certainly wrong (https://physics.stackexchange.com/a/400431/77366).
     # See https://physics.stackexchange.com/a/468649/77366 and
     # https://physics.stackexchange.com/a/335828/77366 for further details on what we
-    # should do (disclaimer: no Sackur-Tetrode!).
+    # should do (disclaimer: no Sackur-Tetrode at 0 K!).
     if np.isclose(temperature, 0.0):
         logger.warning("assuming translational entropy zero at zero temperature")
         return 0.0
@@ -111,8 +109,8 @@ def calc_trans_entropy(
     return translational_entropy
 
 
-# TODO(schneiderfelipe): "energy" has potentially two meanings here. Solve for
-# the whole package.
+# TODO(schneiderfelipe): "energy" has potentially two meanings here. Correct
+# naming for the whole package?
 def calc_internal_energy(
     energy=0.0,
     degeneracy=1,
@@ -178,8 +176,8 @@ def calc_internal_energy(
     return internal_energy
 
 
-# TODO(schneiderfelipe): "energy" has potentially two meanings here. Solve for
-# the whole package.
+# TODO(schneiderfelipe): "energy" has potentially two meanings here. Correct
+# naming for the whole package?
 def calc_enthalpy(
     energy=0.0,
     degeneracy=1,
@@ -251,11 +249,8 @@ def calc_enthalpy(
     return enthalpy
 
 
-# TODO(schneiderfelipe): this function should probably go to _gas.py, but it
-# still calls a function from _thermo._solv. As such, we need to separate
-# things and make a transfer.
-# TODO(schneiderfelipe): "energy" has potentially two meanings here. Solve for
-# the whole package.
+# TODO(schneiderfelipe): "energy" has potentially two meanings here. Correct
+# naming for the whole package?
 def calc_entropy(
     atommasses,
     atomnos=None,

overreact/_thermo/_gas.py

@@ -170,8 +170,6 @@ def calc_elec_entropy(energy=0.0, degeneracy=1, temperature=298.15):
     13.175
 
     """
-    # TODO(schneiderfelipe): This is probably an ugly hack for zero temperature and
-    # probably wrong.
     if np.isclose(temperature, 0.0):
         logger.warning("assuming electronic entropy zero at zero temperature")
         return 0.0
@@ -361,8 +359,6 @@ def calc_rot_entropy(
     ...                  atomcoords=data.atomcoords)
     47.1
     """
-    # TODO(schneiderfelipe): This is probably an ugly hack for zero temperature and
-    # probably wrong.
     if np.isclose(temperature, 0.0):
         logger.warning("assuming rotational entropy zero at zero temperature")
         return 0.0
@@ -548,8 +544,6 @@ def calc_vib_entropy(vibfreqs=None, qrrho=True, temperature=298.15):
     0.0
 
     """
-    # TODO(schneiderfelipe): This is probably an ugly hack for zero temperature and
-    # probably wrong.
     if np.isclose(temperature, 0.0):
         logger.warning("assuming vibrational entropy zero at zero temperature")
         return 0.0
@@ -766,7 +760,7 @@ def _check_vibfreqs(vibfreqs=None, cutoff=-50.0):
     return np.abs(vibfreqs[vibfreqs > cutoff])
 
 
-# TODO(schneiderfelipe): B_av chosen as 1.0e-44 / (atomic_mass * angstrom**2)
+# B_av was chosen as 1.0e-44 / (atomic_mass * angstrom**2)
 def _vibrational_moment(vibfreqs=None, B_av=602.2140762081121):
     """Calculate moments of inertia for a free rotors with the same frequencies.
 
@@ -795,19 +789,19 @@ def _vibrational_moment(vibfreqs=None, B_av=602.2140762081121):
     """
     # TODO(schneiderfelipe): should we receive vibrational temperatures and
     # avoid calling it twice when calling calc_vib_entropy?
-    # TODO(schneiderfelipe): I sincerely have absolutely no clue as to why we
-    # should divide by np.pi and not by 2 * np.pi. The original paper says
-    # nothing about it (and the equation 4 there is probably wrong), but the
-    # expression below reproduces both figure 2 of the original paper and the
-    # results returned by ORCA.
+
+    # I sincerely have absolutely no clue as to why we should divide by np.pi
+    # and not by 2 * np.pi. The original paper says nothing about it (and the
+    # equation 4 there is probably wrong), but the expression below reproduces
+    # both Figure 2 of the original paper and the results returned by ORCA.
     moments = constants.hbar ** 2 / (
         2.0 * constants.k * _vibrational_temperature(vibfreqs) / np.pi
     )
     moments = moments / (constants.atomic_mass * constants.angstrom ** 2)
     return moments * B_av / (moments + B_av)
 
 
-# TODO(schneiderfelipe): omega chosen as (k * 298.15 / (2.0 * h)) * centi / c
+# omega was chosen as (k * 298.15 / (2.0 * h)) * centi / c
 def _head_gordon_damping(vibfreqs, omega=103.61231288246945, alpha=4):
     """Calculate the Head-Gordon damping function.
 

overreact/_thermo/_solv.py

@@ -18,17 +18,17 @@
 # associated with this (see doi:10.1021/jp9716997 and references therein, e.g.,
 # doi:10.1016/0009-2614(96)00349-1, doi:10.1021/cr60304a002,
 # doi:10.1002/jcc.540100504). As such, including this in the total free energy
-# might overcount energy contributions. As an alternative, I might:
-# a. want to remove this contribution,
-# b. want to keep this and remove the contribution from the cavity enthalpy,
+# might overcount energy contributions. As an alternative, I might want to
+# a. remove this contribution,
+# b. keep this and remove the contribution from the cavity enthalpy,
 #    either by
 #    i.  doing a similar calculation as here and removing from the enthalpy or
 #    ii. by actually implementing the original method of C-PCM and excluding
 #        it from the enthalpy.
-# c. want to implement something closer to the original and do one of the
+# c. implement something closer to the original and do one of the
 #    cited in b. above.
 #
-# TODO(schneiderfelipe): see doi:10.1021/cr60304a002.
+# See doi:10.1021/cr60304a002.
 def calc_cav_entropy(
     atomnos,
     atomcoords,
@@ -88,7 +88,7 @@ def calc_cav_entropy(
 
     def func(temperature, solvent):
         # TODO(schneiderfelipe): allow passing a "solvent" object everywhere so
-        # that we don't repeat ourselves.
+        # that we don't repeat ourselves?
         _, _, ratio = coords._garza(
             vdw_volume,
             solvent,
@@ -242,13 +242,8 @@ def molar_free_volume(
         r_M, r_cav = np.cbrt(vdw_volume), np.cbrt(cav_volume)
         molar_free_volume = (r_cav - r_M) ** 3 * constants.angstrom ** 3 * constants.N_A
     elif method == "garza":
-        # TODO(schneiderfelipe): test for the following solvents at the
-        # following temperatures:
-        # - water: 274 K -- 373 K
-        # - pentane: 144 K -- 308 K
-        # - hexane: 178 K -- 340 K
-        # - heptane: 183 K -- 370 K
-        # - octane: 217 K -- 398 K
+        # TODO(schneiderfelipe): test for the following solvents: water,
+        # pentane, hexane, heptane and octane.
         cav_volume, N_cav, _ = coords._garza(
             coords.get_molecular_volume(atomnos, atomcoords),
             environment,

overreact/api.py

@@ -278,7 +278,7 @@ def get_k(
     molecularity=None,
     volume=None,
 ):
-    """Obtain reaction rate constants.
+    r"""Obtain reaction rate constants.
 
     Parameters
     ----------
@@ -287,7 +287,7 @@ def get_k(
     bias : array-like, optional
         Energy to be added to free energies.
     tunneling : str or None, optional
-        Choose between "eckart", "wigner" or None.
+        Choose between "eckart", "wigner" or None (or "none").
     qrrho : bool, optional
         Apply both the quasi-rigid rotor harmonic oscilator (QRRHO)
         approximations of M. Head-Gordon (enthalpy correction, see
@@ -337,6 +337,9 @@ def get_k(
     >>> model = parse_model("data/tanaka1996/UMP2/6-311G(2df,2pd)/model.jk")
     >>> get_k(model.scheme, model.compounds, tunneling=None)
     array([14820222.69476697])
+    >>> get_k(model.scheme, model.compounds, tunneling="none") \
+    ... == get_k(model.scheme, model.compounds, tunneling=None)
+    array([ True])
     >>> get_k(model.scheme, model.compounds, tunneling=None, scale="atm-1 s-1")
     array([605762.44228638])
     >>> get_k(model.scheme, model.compounds, tunneling=None,
@@ -384,7 +387,7 @@ def get_k(
         "(classical) reaction rate constants: "
         f"{', '.join([f'{v:7.3g}' for v in k])} atm⁻ⁿ⁺¹·s⁻¹"
     )
-    if tunneling is not None:
+    if tunneling not in {"none", None}:
         if compounds is not None:
             kappa = get_kappa(
                 scheme, compounds, method=tunneling, temperature=temperature
@@ -403,7 +406,7 @@ def get_k(
         )
 
     # TODO(schneiderfelipe): ensure diffusional limit for reactions in
-    # solvation using Collins-Kimball theory.
+    # solvation using Collins-Kimball theory. This includes half-equilibria.
     return rates.convert_rate_constant(
         k, scale, molecularity=molecularity, temperature=temperature, pressure=pressure
     )
@@ -421,8 +424,8 @@ def get_kappa(scheme, compounds, method="eckart", temperature=298.15):
     ----------
     scheme : Scheme
     compounds : dict-like
-    tunneling : str or None, optional
-        Choose between "eckart" or "wigner".
+    method : str or None, optional
+        Choose between "eckart", "wigner" or None (or "none").
     temperature : array-like, optional
         Absolute temperature in Kelvin.
 
@@ -442,8 +445,13 @@ def get_kappa(scheme, compounds, method="eckart", temperature=298.15):
     >>> kappa = get_kappa(model.scheme, model.compounds)
     >>> kappa
     array([1.10949425])
+    >>> get_kappa(model.scheme, model.compounds, method="none")
+    array([1.0])
     >>> kappa * get_k(model.scheme, model.compounds, tunneling=None)
     array([8.e+10])
+    >>> get_kappa(model.scheme, model.compounds, method="none") \
+    ... == get_kappa(model.scheme, model.compounds, method=None)
+    array([ True])
 
     Beware that the values below have not been validated yet:
 
@@ -486,6 +494,8 @@ def get_kappa(scheme, compounds, method="eckart", temperature=298.15):
                         kappa = tunnel.wigner(vibfreq, temperature=temperature)
             elif method == "wigner":
                 kappa = tunnel.wigner(vibfreq, temperature=temperature)
+            elif method in {"none", None}:
+                kappa = 1.0
             else:
                 raise ValueError(f"unavailable method: '{method}'")