Tau calculation working

examples/semiclassical_micelles.py

@@ -8,31 +8,27 @@
 
 import radicalpy as rp
 from radicalpy.data import Molecule
+from radicalpy.estimations import (autocorrelation, autocorrelation_fit,
+                                   exchange_interaction_in_solution_MC, k_STD)
 from radicalpy.experiments import semiclassical_mary
+from radicalpy.kinetics import Haberkorn
+from radicalpy.relaxation import SingletTripletDephasing
 from radicalpy.simulation import SemiclassicalSimulation, State
+from radicalpy.utils import Bhalf_fit, read_trajectory_files
 
 
-def main(data_path="./examples/data/md_fad_trp_aot"):
-    flavin = Molecule.all_nuclei("flavin_anion")
-    trp = Molecule.all_nuclei("tryptophan_cation")
-    sim = SemiclassicalSimulation([flavin, trp], basis="Zeeman")
-
-    all_data = rp.utils.read_trajectory_files(data_path)
-
-    time = np.linspace(0, len(all_data), len(all_data)) * 5e-12 * 1e9
-    j = rp.estimations.exchange_interaction_in_solution_MC(all_data[:, 1], J0=5)
-
+def plot_exchange_interaction_in_solution(ts, trajectory_data, j):
     fig = plt.figure(1)
     ax = fig.add_axes([0, 0, 1, 1])
     ax.set_facecolor("none")
     ax.grid(False)
     plt.axis("on")
     plt.rc("axes", edgecolor="black")
     color = "tab:red"
-    plt.plot(time, all_data[:, 1] * 1e9, color=color)
+    plt.plot(ts, trajectory_data[:, 1] * 1e9, color=color)
     ax2 = ax.twinx()
     color2 = "tab:blue"
-    plt.plot(time, -j, color=color2)
+    plt.plot(ts, -j, color=color2)
     ax.set_xlabel("Time (ns)", size=24)
     ax.set_ylabel("Radical pair separation (nm)", size=24, color=color)
     ax2.set_ylabel("Exchange interaction (mT)", size=24, color=color2)
@@ -42,17 +38,8 @@ def main(data_path="./examples/data/md_fad_trp_aot"):
     fig.set_size_inches(7, 5)
     plt.show()
 
-    # Calculate the autocorrelation, tau_c, and k_STD
-    acf_j = rp.utils.autocorrelation(j, factor=1)
-    zero_point_crossing_j = np.where(np.diff(np.sign(acf_j)))[0][0]
-    t_j_max = max(time[:zero_point_crossing_j]) * 1e-9
-    t_j = np.linspace(5e-12, t_j_max, zero_point_crossing_j)
-
-    acf_j_fit = rp.estimations.autocorrelation_fit(t_j, j, 5e-12, t_j_max)
-    acf_j_fit["tau_c"]
-    kstd = rp.estimations.k_STD(j, acf_j_fit["tau_c"])
-    k_STD = np.mean(kstd)  # singlet-triplet dephasing rate
 
+def plot_autocorrelation_fit(t_j, acf_j, acf_j_fit, zero_point_crossing_j):
     fig = plt.figure(2)
     ax = fig.add_axes([0, 0, 1, 1])
     ax.set_facecolor("none")
@@ -68,30 +55,54 @@ def main(data_path="./examples/data/md_fad_trp_aot"):
     fig.set_size_inches(7, 5)
     plt.show()
 
-    kq = 5e6  # triplet excited state quenching rate
-    krec = 8e6  # recombination rate
-    kesc = 5e5  # escape rate
 
-    time = np.arange(0, 10e-6, 10e-9)
+def main(data_path="./examples/data/md_fad_trp_aot"):
+    flavin = Molecule.all_nuclei("flavin_anion")
+    trp = Molecule.all_nuclei("tryptophan_cation")
+    sim = SemiclassicalSimulation([flavin, trp], basis="Zeeman")
+    print(f"{sim.molecules[0].semiclassical_tau=}")
+    print(f"{sim.molecules[1].semiclassical_tau=}")
+
+    trajectory_data = read_trajectory_files(data_path)
+    ts = np.linspace(0, len(trajectory_data), len(trajectory_data)) * 5e-12 * 1e9
+    j = exchange_interaction_in_solution_MC(trajectory_data[:, 1], J0=5)
+
+    # plot_exchange_interaction_in_solution(ts, trajectory_data, j)
+
+    # Calculate the autocorrelation, tau_c, and k_STD
+    acf_j = autocorrelation(j, factor=1)
+    zero_point_crossing_j = np.where(np.diff(np.sign(acf_j)))[0][0]
+    t_j_max = max(ts[:zero_point_crossing_j]) * 1e-9
+    t_j = np.linspace(5e-12, t_j_max, zero_point_crossing_j)
+
+    acf_j_fit = autocorrelation_fit(t_j, j, 5e-12, t_j_max)
+    acf_j_fit["tau_c"]
+    kstd = rp.estimations.k_STD(j, acf_j_fit["tau_c"])
+    # k_STD = np.mean(kstd)  # singlet-triplet dephasing rate ##################################
+
+    # plot_autocorrelation_fit(t_j, acf_j, acf_j_fit, zero_point_crossing_j)
+
+    triplet_excited_state_quenching_rate = 5e6
+    recombination_rate = 8e6
+    free_radical_escape_rate = 5e5
+
+    ts = np.arange(0, 10e-6, 10e-9)
     Bs = np.arange(0, 50, 1)
 
-    kinetics = [
-        rp.kinetics.Haberkorn(krec, State.SINGLET),
-        rp.kinetics.FreeRadical(kesc),
-        rp.kinetics.ElectronTransfer(kq),
-    ]
-    relaxations = [rp.relaxation.SingletTripletDephasing(kstd)]
     results = semiclassical_mary(
+        sim=sim,
+        num_samples=10,
         init_state=State.TRIPLET,
-        obs_state=State.TRIPLET,
-        time=time,
+        # obs_state=State.TRIPLET,
+        time=ts,
         B0=Bs,
         D=0,
         J=0,
-        kinetics=kinetics,
-        relaxations=relaxations,
+        triplet_excited_state_quenching_rate=triplet_excited_state_quenching_rate,
+        free_radical_escape_rate=free_radical_escape_rate,
+        kinetics=[Haberkorn(recombination_rate, State.SINGLET)],
+        relaxations=[SingletTripletDephasing(kstd)],
     )
-    results = sim.MARY()
 
     # Calculate time evolution of the B1/2
     bhalf_time = np.zeros((len(results)))
@@ -105,16 +116,16 @@ def main(data_path="./examples/data/md_fad_trp_aot"):
             fit_time[:, i],
             fit_error_time[:, i],
             R2_time[i],
-        ) = rp.utils.Bhalf_fit(B, results["MARY"])
+        ) = Bhalf_fit(Bs, results["MARY"])
 
     # Plotting
     factor = 1e6
 
     plt.figure(3)
-    for i in range(2, len(time), 35):
-        plt.plot(time[i] * factor, bhalf_time[i], "ro", linewidth=3)
+    for i in range(2, len(ts), 35):
+        plt.plot(ts[i] * factor, bhalf_time[i], "ro", linewidth=3)
         plt.errorbar(
-            time[i] * factor,
+            ts[i] * factor,
             bhalf_time[i],
             fit_error_time[1, i],
             color="k",
@@ -123,15 +134,15 @@ def main(data_path="./examples/data/md_fad_trp_aot"):
     plt.xlabel("Time ($\mu s$)", size=18)
     plt.ylabel("$B_{1/2}$ (mT)", size=18)
     plt.tick_params(labelsize=14)
-    fig.set_size_inches(10, 5)
+    plt.gcf().set_size_inches(10, 5)
     plt.show()
 
     fig = plt.figure(figsize=plt.figaspect(1.0))
     ax = fig.add_subplot(projection="3d")
     cmap = plt.cm.ScalarMappable(cmap=plt.get_cmap("viridis"))
     ax.set_facecolor("none")
     ax.grid(False)
-    X, Y = np.meshgrid(Bs, time)
+    X, Y = np.meshgrid(Bs, ts)
     ax.plot_surface(
         X,
         Y * factor,

radicalpy/data.py

@@ -1,7 +1,7 @@
 #! /usr/bin/env python
 import json
 from functools import singledispatchmethod
-from typing import Optional
+from typing import Iterator, Optional, Tuple
 
 import numpy as np
 from importlib_resources import files
@@ -569,6 +569,40 @@ def effective_hyperfine(self) -> float:
         hfcs_np = np.array([h.isotropic for h in hfcs])
         return np.sqrt((4 / 3) * sum((hfcs_np**2 * spns_np) * (spns_np + 1)))
 
+    @property
+    def semiclassical_tau(self) -> float:
+        """Calculate the `tau` coefficient used in `Molecule.semiclassical_random_hfc_theta_phi_rng`.
+
+        Examples:
+            >>> m = Molecule.fromdb("flavin_anion", nuclei=["N14"])
+            >>> m.semiclassical_tau
+        """
+        tmp = sum(
+            n.spin_quantum_number * (n.spin_quantum_number + 1) * n.hfc.isotropic**2
+            for n in self.nuclei
+        )
+        return (tmp / 6) ** -0.5
+
+    def semiclassical_random_hfc_theta_phi_rng(
+        self, num_samples: int, I_max: float, fI_max: float
+    ) -> Iterator[Tuple[float, float, float]]:
+        hfc_rng = np.random.default_rng(42)
+        ang_rng = np.random.default_rng(43)
+        rng = np.random.default_rng(42)
+        tau = self.semiclassical_tau
+        for j in range(num_samples):
+            while True:
+                I_r = rng.random() * I_max
+                fI_r = rng.random() * fI_max
+                f = (
+                    I_r**2
+                    * ((tau**2 / (4 * np.pi)) ** 1.5)
+                    * np.exp(-1 / 4 * I_r**2 * tau**2)
+                )
+                if fI_r < f:
+                    yield I_r
+                break
+
 
 class Triplet(Molecule):
     def __init__(self):

radicalpy/experiments.py

@@ -1,9 +1,11 @@
 #!/usr/bin/env python
 
 import numpy as np
+from numpy.typing import NDArray
 from tqdm import tqdm
 
-from .simulation import HilbertIncoherentProcessBase, LiouvilleSimulation, State
+from .simulation import (HilbertIncoherentProcessBase, LiouvilleSimulation,
+                         SemiclassicalSimulation, State)
 from .utils import mary_lorentzian, modulated_signal, reference_signal
 
 
@@ -46,7 +48,7 @@ def modulated_mary_brute_force(
 def steady_state_mary(
     sim: LiouvilleSimulation,
     obs: State,
-    Bs: np.ndarray,
+    Bs: NDArray[float],
     D: float,
     E: float,
     J: float,
@@ -71,32 +73,24 @@ def steady_state_mary(
 
 
 def semiclassical_mary(
-    sim, num_samples, init_state, obs_state, time, B0, D, J, kinetics, relaxations
+    sim: SemiclassicalSimulation,
+    num_samples: int,
+    init_state: State,
+    time: NDArray[float],
+    B0,
+    D,
+    J,
+    triplet_excited_state_quenching_rate,
+    free_radical_escape_rate,
+    kinetics,
+    relaxations,
 ):
+    print(f"{sim.molecules[0].semiclassical_tau=}")
+    print(f"{sim.molecules[1].semiclassical_tau=}")
+    exit()
     for i, B0 in enumerate(tqdm(B)):
-        for j in range(0, sample_number):
-            FAD_omegax = (
-                FAD_gamma * FAD_I[j] * np.sin(FAD_theta[j]) * np.cos(FAD_phi[j])
-            )
-            Trp_omegax = (
-                Trp_gamma * Trp_I[j] * np.sin(Trp_theta[j]) * np.cos(Trp_phi[j])
-            )
-            hamiltonian_x = FAD_omegax * FAD_Sx + Trp_omegax * Trp_Sx
-
-            FAD_omegay = (
-                FAD_gamma * FAD_I[j] * np.sin(FAD_theta[j]) * np.sin(FAD_phi[j])
-            )
-            Trp_omegay = (
-                Trp_gamma * Trp_I[j] * np.sin(Trp_theta[j]) * np.sin(Trp_phi[j])
-            )
-            hamiltonian_y = FAD_omegay * FAD_Sy + Trp_omegay * Trp_Sy
-
-            FAD_omegaz = FAD_gamma * (B0 + FAD_I[j] * np.cos(FAD_theta[j]))
-            Trp_omegaz = Trp_gamma * (B0 + Trp_I[j] * np.cos(Trp_theta[j]))
-            hamiltonian_z = FAD_omegaz * FAD_Sz + Trp_omegaz * Trp_Sz
-
-            HZ = hamiltonian_x + hamiltonian_y + hamiltonian_z
-            HZL = rp.simulation.LiouvilleSimulation.convert(HZ)
+        gen = sim.semiclassical_gen(num_samples)
+        for HZL in gen:
             sim.apply_liouville_hamiltonian_modifiers(HZL, kinetics + relaxation)
             L = HZL
 
@@ -106,14 +100,12 @@ def semiclassical_mary(
             triplet_initial_population = 1  # triplet excited state
 
             initial_temp = np.reshape(initial / 3, (M, 1))
-            initial_density = np.reshape(np.zeros(16), (M, 1))
-            density = initial_density
+            density = np.reshape(np.zeros(16), (M, 1))
 
             for k in range(0, len(time)):
                 FR_density = density
                 population = np.trace(FR_density)
-                free_radical = FR_initial_population + population * kesc * dt
-                rho = population + free_radical
+                rho = population + (FR_initial_population + population * kesc * dt)
                 trace[j, k] = np.real(rho)
                 density = (
                     propagator * density
@@ -123,7 +115,7 @@ def semiclassical_mary(
                     -kq * dt
                 )
 
-        average = np.ones(sample_number) * 0.01
+        average = np.ones(num_samples) * 0.01
         decay = average @ trace
         if i == 0:
             decay0 = np.real(decay)

radicalpy/simulation.py

@@ -3,10 +3,11 @@
 import enum
 import itertools
 from math import prod
-from typing import Optional
+from typing import Iterator, Optional
 
 import numpy as np
 import scipy as sp
+from numpy.typing import NDArray
 from tqdm import tqdm
 
 from . import utils
@@ -1205,6 +1206,9 @@ def adjust_hamiltonian(self, H: np.ndarray):
 
 
 class SemiclassicalSimulation(LiouvilleSimulation):
+    def semiclassical_gen(num_samples) -> Iterator[NDArray[float]]:
+        pass
+
     @property
     def nuclei(self):
         return []