From 20d08848c0d30ac7bbb9218f4185356dcfceb1d1 Mon Sep 17 00:00:00 2001 From: Johannes Steinmetzer Date: Mon, 21 Aug 2023 13:33:24 +0200 Subject: [PATCH] add: drivers for Marcudim-scan and -parametrization --- pysisyphus/drivers/marcusdim.py | 7 +- pysisyphus/drivers/marcusdim_param.py | 151 ++++++++++++++++++++++ pysisyphus/drivers/marcusdim_scan.py | 128 ++++++++++++++++++ tests/test_marcusdim/test_param_marcus.py | 41 ++++++ 4 files changed, 325 insertions(+), 2 deletions(-) create mode 100644 pysisyphus/drivers/marcusdim_param.py create mode 100644 pysisyphus/drivers/marcusdim_scan.py create mode 100644 tests/test_marcusdim/test_param_marcus.py diff --git a/pysisyphus/drivers/marcusdim.py b/pysisyphus/drivers/marcusdim.py index bf4f2a313..f39325a20 100644 --- a/pysisyphus/drivers/marcusdim.py +++ b/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 diff --git a/pysisyphus/drivers/marcusdim_param.py b/pysisyphus/drivers/marcusdim_param.py new file mode 100644 index 000000000..9c68019c2 --- /dev/null +++ b/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] diff --git a/pysisyphus/drivers/marcusdim_scan.py b/pysisyphus/drivers/marcusdim_scan.py new file mode 100644 index 000000000..a40b13f4d --- /dev/null +++ b/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 diff --git a/tests/test_marcusdim/test_param_marcus.py b/tests/test_marcusdim/test_param_marcus.py new file mode 100644 index 000000000..f74ede456 --- /dev/null +++ b/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()