diff --git a/examples/semiclassical_micelles.py b/examples/semiclassical_micelles.py index c067779..88035f4 100644 --- a/examples/semiclassical_micelles.py +++ b/examples/semiclassical_micelles.py @@ -8,20 +8,16 @@ 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") @@ -29,10 +25,10 @@ def main(data_path="./examples/data/md_fad_trp_aot"): 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,7 +134,7 @@ 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)) @@ -131,7 +142,7 @@ def main(data_path="./examples/data/md_fad_trp_aot"): 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, diff --git a/radicalpy/data.py b/radicalpy/data.py index b84dcb7..9a50aa8 100644 --- a/radicalpy/data.py +++ b/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): diff --git a/radicalpy/experiments.py b/radicalpy/experiments.py index 46adff6..3782d67 100644 --- a/radicalpy/experiments.py +++ b/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) diff --git a/radicalpy/simulation.py b/radicalpy/simulation.py index 23f8201..3a7066c 100644 --- a/radicalpy/simulation.py +++ b/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 []