add: drivers for Marcudim-scan and -parametrization

pysisyphus/drivers/marcusdim.py

@@ -1,8 +1,11 @@
 # [1] https://doi.org/10.26434/chemrxiv-2022-253hc-v2
 #     Identifying the Marcus dimension of electron transfer from
 #     ab initio calculations
-#     Šrut, Lear, Krewald, 2022
-
+#     Šrut, Lear, Krewald, 2022 on chemrxiv
+# [2] https://doi.org/10.1039/D3SC01402A
+#     Identifying the Marcus dimension of electron transfer from
+#     ab initio calculations
+#     Šrut, Lear, Krewald, 2023, actual published version
 
 import datetime
 from enum import Enum

pysisyphus/drivers/marcusdim_param.py

@@ -0,0 +1,151 @@
+# [1] https://doi.org/10.1039/D3SC01402A
+#     Identifying the Marcus dimension of electron transfer from
+#     ab initio calculations
+#     Šrut, Lear, Krewald, 2023
+
+
+from dataclasses import dataclass
+
+import matplotlib.pyplot as plt
+import sympy as sym
+
+from pysisyphus.constants import BOHR2ANG, NU2AU
+
+
+@dataclass
+class MarcusModel:
+    reorg_en: float  # Lambda, Hartree
+    dG: float  # ΔG, barrier height, in Hartree
+    coupling: float  # Vab, in Hartree
+    R: float  # Distance of adiabatic minimum to top of barrier in Bohr
+    f: float  # Force constant in Hartree / Bohr**2
+    d: float  # Separation of diabatic states in Bohr
+
+    def as_wavenums_and_ang_tuple(self):
+        return (
+            self.reorg_en / NU2AU,
+            self.dG / NU2AU,
+            self.coupling / NU2AU,
+            self.R * BOHR2ANG,
+            self.f / NU2AU / BOHR2ANG**2,
+            self.d * BOHR2ANG,
+        )
+
+    def G_diabatic(self, x):
+        Ga = self.f * x**2
+        Gb = self.f * (x - self.d) ** 2
+        return Ga, Gb
+
+    def plot_diabatic(self, x, show=False):
+        Ga, Gb = self.G_diabatic(x)
+        fig, ax = plt.subplots()
+        for state in (Ga, Gb):
+            ax.plot(x, state)
+        if show:
+            plt.show()
+        return fig, ax
+
+    def pretty(self):
+        reorg_en, dG, coupling, *_, d = self.as_wavenums_and_ang_tuple()
+        reorg_en = f"{reorg_en:.0f} cm⁻¹"
+        dG = f"{dG:.0f} cm⁻¹"
+        _2coupling = f"{2*coupling:.0f} cm⁻¹"
+        d = f"{d:.3f} Å"
+        return f"MarcusModel(λ={reorg_en}, ΔG={dG}, 2Vab={_2coupling}, d={d})"
+
+
+def solve_marcus(R, Vab, dG=None, en_exc_min=None):
+    """Fit Marcus-Hush model to given parameters R, Vab, dG and en_exc_min.
+    Either dG or en_exc_min is mandatory.
+
+    All argument should be given in atomic units (Hartree and Bohr).
+
+    Quartic extension is described here:
+        https://www.science.org/doi/epdf/10.1126/science.278.5339.846
+    This is not yet implemented!"""
+
+    assert (dG is not None) or (
+        en_exc_min is not None
+    ), "Either 'dG' or 'en_exc_min' must be given!"
+    f, d = sym.symbols("f d", real=True)
+    # Equation 1 is the same for both parametrizations; either eq. (5b) or (6b)
+    eq1 = (
+        d / 2
+        - 1 / 2 * ((d**2 * f - sym.sqrt(d**4 * f**2 - 4 * Vab**2)) / (d * f))
+        - R
+    )
+
+    def inner(eq2):
+        # Solve equation system using sympy
+        results = sym.solve([eq1, eq2], d, f, dict=True)
+        assert len(results) == 1
+        res = results[0]
+        fval = res[f]
+        dval = res[d]
+        reorg_en = fval * dval**2
+        dG = (dval**2 * fval - 2 * Vab) ** 2 / (4 * dval**2 * fval)
+        return MarcusModel(
+            reorg_en=float(reorg_en),
+            dG=float(dG),
+            coupling=float(Vab),
+            R=float(R),
+            f=float(fval),
+            d=float(dval),
+        )
+
+    results = dict()
+    # Use parametrization A
+    if dG is not None:
+        # Eq. (5c) in SI of [1]
+        eq2_a = (d**2 * f - 2 * Vab) ** 2 / (4 * d**2 * f) - dG
+        results["a"] = inner(eq2_a)
+    # Use parametrization B
+    if en_exc_min is not None:
+        # Eq. (6c) in SI of [1]
+        eq2_b = d**2 * f - en_exc_min
+        results["b"] = inner(eq2_b)
+    # If dG and en_exc_min are given we can use both parametrizations.
+    return results
+
+
+def solve_marcus_wavenums_and_ang(R, Vab, dG=None, en_exc_min=None):
+    """Wrapper for solve_marcus that accepts wavenumbers and Ångstrom."""
+    kwargs = {
+        "R": R / BOHR2ANG,
+        "Vab": Vab * NU2AU,
+        "dG": dG * NU2AU if dG is not None else None,
+        "en_exc_min": en_exc_min * NU2AU if en_exc_min is not None else None,
+    }
+    return solve_marcus(**kwargs)
+
+
+def find_minima(arr):
+    assert (arr.ndim == 1) and (len(arr) >= 3)
+    return [i for i in range(1, arr.size) if arr[i - 1] > arr[i] < arr[i + 1]]
+
+
+def param_marcus(coordinate, energies, scheme="B"):
+    assert energies.ndim == 2
+
+    # Excitation energy at adiabatic minimum
+    min_inds = find_minima(energies[:, 0])  # Search minima in lower state
+    if len(min_inds) == 2:
+        ind0, ind1 = min_inds
+        adia_min_ind = ind0 if energies[ind0, 0] < energies[ind1, 0] else ind1
+        barr_ind = energies[ind0 : ind1 + 1].argmax() + ind0
+    elif len(min_inds) == 1:
+        adia_min_ind = barr_ind = min_inds[0]
+    else:
+        raise Exception("How did I get here?!")
+
+    adia_min_ens = energies[adia_min_ind]
+    en_exc_min = adia_min_ens[1] - adia_min_ens[0]
+    barr_ens = energies[barr_ind]
+    en_exc_barr = barr_ens[1] - barr_ens[0]
+
+    # Electronic coupling
+    V_ab = en_exc_barr / 2
+    # Distance R between adiabatic minimum and top of barrier
+    R = coordinate[barr_ind] - coordinate[adia_min_ind]
+    # Barrier height ΔG
+    dG = energies[barr_ind, 0]

pysisyphus/drivers/marcusdim_scan.py

@@ -0,0 +1,128 @@
+import math
+import sys
+import warnings
+
+import numpy as np
+
+from pysisyphus.helpers import rms
+
+
+def scan_dir(
+    x0,
+    direction,
+    get_property,
+    step_size=0.05,
+    add_steps=10,
+    max_steps=500,
+    min_steps=5,
+    rms_grad_thresh=1e-2,
+):
+    assert add_steps >= 0, f"{add_steps=:} must be >= 0!"
+    assert max_steps > 0, f"{max_steps=:} must be positive!"
+
+    grad = np.nan
+    step = step_size * direction
+    stop_in = add_steps
+    converged = False
+
+    xcur = x0 + step
+    prop_prev = None
+    grad_rms_prev = None
+
+    all_factors = np.arange(max_steps) + 1
+    all_coords = np.empty((max_steps, *x0.shape))
+    all_energies = list()
+    all_properties = np.empty(max_steps)
+    for i in range(max_steps):
+        all_coords[i] = xcur
+        # Calculate & store property
+        energies, prop = get_property(xcur)
+        all_properties[i] = prop
+        all_energies.append(energies)
+
+        # Determine gradient from finite differences
+        if prop_prev is not None:
+            grad = (prop - prop_prev) / step_size
+        print(f"{i=:03d}, property={prop: >12.4f}, {grad=: >12.4f}")
+        sys.stdout.flush()
+
+        rms_grad = rms(grad)
+        # Break when the gradient increases. This is probably enough and the
+        # following rms(grad) check is not needed.
+        if (
+            (i >= min_steps)
+            and (grad_rms_prev is not None)
+            and rms_grad > grad_rms_prev
+        ):
+            print("Property gradient increased. Breaking")
+            break
+
+        # Also check convergence via rms; this check is skipped once convergence
+        # is indicated.
+        if not converged and (converged := np.abs(rms_grad) <= rms_grad_thresh):
+            print("Converged!")
+
+        # After convergence we can carry out some additional steps, if requested.
+        if add_steps and converged:
+            stop_in -= 1
+
+        # Break directly when we don't want to do any additional steps
+        if converged and add_steps == 0:
+            break
+        elif add_steps and stop_in < 0:
+            print("Did additional steps")
+            break
+
+        # Update variables
+        xcur = xcur + step
+        prop_prev = prop
+        grad_rms_prev = rms_grad
+    else:
+        raise Exception("Scan did not converge!")
+    all_energies = np.array(all_energies)
+    # Truncate arrays and drop empty part. This will also drop the last calculation
+    # that lead to the break from the loop.
+    return (
+        all_factors[:i] * step_size,
+        all_coords[:i],
+        all_energies[:i],
+        all_properties[:i],
+    )
+
+
+def scan(x0, direction, get_properties, **kwargs):
+    dir_norm = np.linalg.norm(direction)
+    if not math.isclose(dir_norm, 1.0):
+        warnings.warn(f"norm(direction)={dir_norm:.6f} is not 1.0! Renormalizing.")
+        direction = direction / dir_norm
+    # Carry out calculation on initial geometry.
+    ens0, prop0 = get_properties(x0)
+
+    def get_property_changes(xi):
+        """Get property changes w.r.t. initial geometry."""
+        ens, prop = get_properties(xi)
+        return ens, prop  # - prop0
+
+    print("Positive direction")
+    pos_dir = direction
+    pos_facts, pos_coords, pos_ens, pos_props = scan_dir(
+        x0, pos_dir, get_property_changes, **kwargs
+    )
+    print()
+
+    print("Negative direction")
+    neg_dir = -1.0 * direction
+    neg_facts, neg_coords, neg_ens, neg_props = scan_dir(
+        x0, neg_dir, get_property_changes, **kwargs
+    )
+
+    def concat(neg, init, pos):
+        return np.concatenate((neg[::-1], [init], pos))
+
+    all_facts = concat(-neg_facts, 0.0, pos_facts)
+    all_coords = concat(neg_coords, x0, pos_coords)
+    all_energies = concat(neg_ens, ens0, pos_ens)
+    # When we consider the difference w.r.t. initial geometry then
+    # the property is always 0.0.
+    all_properties = concat(neg_props, prop0, pos_props)
+    return all_facts, all_coords, all_energies, all_properties

tests/test_marcusdim/test_param_marcus.py

@@ -0,0 +1,41 @@
+import math
+
+import numpy as np
+
+from pysisyphus.drivers.marcusdim_param import solve_marcus_wavenums_and_ang
+
+
+def assert_parametrization(model, lambda_ref, dG_ref, _2Vab_ref, d_ref):
+    *energies, _, _, d = model.as_wavenums_and_ang_tuple()
+    np.testing.assert_allclose(
+        energies, (lambda_ref, dG_ref, _2Vab_ref / 2.0), atol=0.5
+    )
+    assert math.isclose(d, d_ref, abs_tol=5e-4)
+
+
+def test_mdnb():
+    # Ab initio scan data from Table 1 in the manuscript for m-DNB·⁻, second to last row.
+    dG_ref = 1233
+    en_exc_min_ref = 7434
+    Vab_ref = 2752 / 2.0
+    res = solve_marcus_wavenums_and_ang(
+        R=0.06, Vab=Vab_ref, dG=dG_ref, en_exc_min=en_exc_min_ref
+    )
+    ma = res["a"]
+    mb = res["b"]
+
+    print("Parametrization A")
+    print("\t", ma.pretty())
+    print("Parametrization B")
+    print("\t", mb.pretty())
+
+    assert_parametrization(ma, 9651, dG_ref, 2 * Vab_ref, 0.125)
+    assert_parametrization(mb, en_exc_min_ref, 737, 2 * Vab_ref, 0.129)
+
+    # import matplotlib.pyplot as plt
+    # xs = np.linspace(-0.10, 0.2)
+    # figa, _ = ma.plot_diabatic(xs, show=False)
+    # figa.suptitle("Model A")
+    # figb, _ = mb.plot_diabatic(xs, show=False)
+    # figb.suptitle("Model B")
+    # plt.show()