diff --git a/examples/sctep.py b/examples/sctep.py new file mode 100644 index 0000000..2eb7bdb --- /dev/null +++ b/examples/sctep.py @@ -0,0 +1,56 @@ +#! /usr/bin/env python + +import matplotlib.pyplot as plt +import numpy as np + +from radicalpy.data import Triplet +from radicalpy.experiments import steady_state_mary +from radicalpy.kinetics import Haberkorn, HaberkornFree +from radicalpy.simulation import Basis, LiouvilleSimulation, State + + +def main( + Bs=np.arange(0, 2500, 10), + D=-6.2, + E=35, + J=499.55, + k0=1, + ks=1.1e9, + kd=2.8e9, +): + m = Triplet() + sim = LiouvilleSimulation(molecules=[m, m], basis=Basis.ZEEMAN) + rhos, Phi_s = steady_state_mary( + sim, + obs=State.TP_SINGLET, + Bs=Bs, + D=D, + E=E, + J=J, + theta=np.pi / 4, + phi=0, + kinetics=[Haberkorn(ks, State.TP_SINGLET), HaberkornFree(kd)], + ) + rhos *= k0 + Phi_s *= k0 * ks + + MFE = ((np.abs(Phi_s) - np.abs(Phi_s[0])) / np.abs(Phi_s[0])) * 100 + # print(H) + 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="k") + plt.plot(Bs / J, MFE, linewidth=3, color="tab:red") + ax.set_xlabel("g$μ_B$$B_0$ / J", size=18) + ax.set_ylabel("MFE (%)", size=18) + ax.axvline(x=1.5, color="k", linestyle="--") + ax.axvline(x=3, color="k", linestyle="--") + plt.tick_params(labelsize=14) + fig.set_size_inches(8, 5) + plt.show() + + +if __name__ == "__main__": + main() diff --git a/pyproject.toml b/pyproject.toml index 05f9fbc..baf9f2d 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -32,16 +32,19 @@ dependencies = {file = ["requirements.txt"]} [tool.setuptools.package-data] "*" = ["*.json"] -[tool.black] -line-length = 88 - [tool.pydocstyle] convention = "google" +[tool.black] +line-length = 88 + [tool.pylint] max-line-length = 88 disable = ["C0103"] +[tool.isort] +line_length = 88 + [tool.mypy] plugins = ["numpy.typing.mypy_plugin"] disallow_untyped_defs = true diff --git a/radicalpy/data.py b/radicalpy/data.py index a036acb..26d1d6d 100644 --- a/radicalpy/data.py +++ b/radicalpy/data.py @@ -515,3 +515,10 @@ def effective_hyperfine(self) -> float: spns_np = np.array(list(map(multiplicity_to_spin, multiplicities))) hfcs_np = np.array([h.isotropic for h in hfcs]) return np.sqrt((4 / 3) * sum((hfcs_np**2 * spns_np) * (spns_np + 1))) + + +class Triplet(Molecule): + def __init__(self): + gamma = Isotope("E").gamma_mT + triplet = Nucleus(magnetogyric_ratio=gamma, multiplicity=3, hfc=0.0) + super().__init__(name="Triplet", nuclei=[], radical=triplet) diff --git a/radicalpy/experiments.py b/radicalpy/experiments.py new file mode 100644 index 0000000..a3a56cc --- /dev/null +++ b/radicalpy/experiments.py @@ -0,0 +1,33 @@ +#!/usr/bin/env python + +import numpy as np +from tqdm import tqdm + +from .simulation import HilbertIncoherentProcessBase, LiouvilleSimulation, State + + +def steady_state_mary( + sim: LiouvilleSimulation, + obs: State, + Bs: np.ndarray, + D: float, + E: float, + J: float, + theta: float, + phi: float, + kinetics: list[HilbertIncoherentProcessBase] = [], + # relaxations: list[HilbertIncoherentProcessBase] = [], +) -> np.ndarray: + HZFS = sim.zero_field_splitting_hamiltonian(D, E) + HJ = sim.exchange_hamiltonian(-J, prod_coeff=1) + rhos = np.zeros(shape=(len(Bs), sim.hamiltonian_size)) + Q = sim.projection_operator(obs) + for i, B in enumerate(tqdm(Bs)): + HZ = sim.zeeman_hamiltonian(B, theta=theta, phi=phi) + H = HZ + HZFS + HJ + H = sim.convert(H) + sim.apply_liouville_hamiltonian_modifiers(H, kinetics) # + relaxations) + rhos[i] = np.linalg.solve(H, Q.flatten()) + + Phi_s = rhos @ Q.flatten() + return rhos, Phi_s diff --git a/radicalpy/kinetics.py b/radicalpy/kinetics.py index 4252a4b..ea02f8f 100644 --- a/radicalpy/kinetics.py +++ b/radicalpy/kinetics.py @@ -1,13 +1,9 @@ -from math import prod +"""Kinetics classes.""" import numpy as np -from .simulation import ( - HilbertIncoherentProcessBase, - LiouvilleIncoherentProcessBase, - LiouvilleSimulation, - State, -) +from .simulation import (HilbertIncoherentProcessBase, LiouvilleIncoherentProcessBase, + LiouvilleSimulation, State) class HilbertKineticsBase(HilbertIncoherentProcessBase): @@ -25,10 +21,6 @@ def _name(self): name = super()._name() return f"Kinetics: {name}" - @staticmethod - def _convert(Q: np.ndarray) -> np.ndarray: - return np.kron(Q, np.eye(len(Q))) + np.kron(np.eye(len(Q)), Q) - class Exponential(HilbertKineticsBase): """Exponential model kinetics operator. @@ -79,7 +71,15 @@ class Haberkorn(LiouvilleKineticsBase): def __init__(self, rate_constant: float, target: State): super().__init__(rate_constant) self.target = target - if target not in {State.SINGLET, State.TRIPLET}: + if target not in { + State.SINGLET, + State.TP_SINGLET, + State.TRIPLET, + State.TRIPLET_MINUS, + State.TRIPLET_PLUS, + State.TRIPLET_PLUS_MINUS, + State.TRIPLET_ZERO, + }: raise ValueError( "Haberkorn kinetics supports only SINGLET and TRIPLET targets" ) @@ -87,7 +87,7 @@ def __init__(self, rate_constant: float, target: State): def init(self, sim: LiouvilleSimulation): """See `radicalpy.simulation.HilbertIncoherentProcessBase.init`.""" Q = sim.projection_operator(self.target) - self.subH = 0.5 * self.rate * self._convert(Q) + self.subH = 0.5 * self.rate * sim._convert(Q) def __repr__(self): lines = [ @@ -112,8 +112,7 @@ class HaberkornFree(LiouvilleKineticsBase): def init(self, sim: LiouvilleSimulation): """See `radicalpy.simulation.HilbertIncoherentProcessBase.init`.""" - size = prod(p.multiplicity for p in sim.particles) ** 2 - self.subH = self.rate * np.eye(size) + self.subH = self.rate * np.eye(sim.hamiltonian_size) class JonesHore(LiouvilleKineticsBase): @@ -143,8 +142,8 @@ def init(self, sim: LiouvilleSimulation): QS = sim.projection_operator(State.SINGLET) QT = sim.projection_operator(State.TRIPLET) self.subH = ( - 0.5 * self.singlet_rate * self._convert(QS) - + 0.5 * self.triplet_rate * self._convert(QT) + 0.5 * self.singlet_rate * sim._convert(QS) + + 0.5 * self.triplet_rate * sim._convert(QT) + (0.5 * (self.singlet_rate + self.triplet_rate)) * (np.kron(QS, QT) + np.kron(QT, QS)) ) diff --git a/radicalpy/simulation.py b/radicalpy/simulation.py index 2b305e9..e432ffc 100644 --- a/radicalpy/simulation.py +++ b/radicalpy/simulation.py @@ -21,6 +21,7 @@ class State(enum.Enum): TRIPLET_PLUS = "T_+" TRIPLET_PLUS_MINUS = "T_\\pm" TRIPLET_MINUS = "T_-" + TP_SINGLET = "TP_S" class Basis(enum.Enum): @@ -104,6 +105,10 @@ def nuclei(self): def particles(self): return self.radicals + self.nuclei + @property + def hamiltonian_size(self): + return np.prod([p.multiplicity for p in self.particles]) + def __repr__(self) -> str: return "\n".join( [ @@ -206,10 +211,10 @@ def spin_operator(self, idx: int, axis: str) -> np.ndarray: assert axis in "xyzpmu" sigma = self.pauli(self.particles[idx].multiplicity)[axis] - eye_before = np.eye(prod(p.multiplicity for p in self.particles[:idx])) - eye_after = np.eye(prod(p.multiplicity for p in self.particles[idx + 1 :])) - - spinop = np.kron(np.kron(eye_before, sigma), eye_after) + before_size = prod(p.multiplicity for p in self.particles[:idx]) + after_size = prod(p.multiplicity for p in self.particles[idx + 1 :]) + spinop = np.kron(np.eye(before_size), sigma) + spinop = np.kron(spinop, np.eye(after_size)) if self.basis == Basis.ST: return self.ST_basis(spinop) else: @@ -305,9 +310,22 @@ def projection_operator(self, state: State): State.TRIPLET_PLUS_MINUS: (2 * SAz**2 + SAz) * (2 * SBz**2 + SBz) + (2 * SAz**2 - SAz) * (2 * SBz**2 - SBz), State.EQUILIBRIUM: 1.05459e-34 / (1.38e-23 * 298), + State.TP_SINGLET: self.tp_singlet_projop(SAx, SAy, SAz, SBx, SBy, SBz), } + return result[state] + def tp_singlet_projop(self, SAx, SAy, SAz, SBx, SBy, SBz): + # For radical triplet pair (RTP) + SAsquared = SAx @ SAx + SAy @ SAy + SAz @ SAz + SBsquared = SBx @ SBx + SBy @ SBy + SBz @ SBz + Ssquared = SAsquared + SBsquared + 2 * (SAx @ SBx + SAy @ SBy + SAz @ SBz) # + return ( + (1 / 12) + * (Ssquared - (6 * np.eye(len(SAx)))) + @ (Ssquared - (2 * np.eye(len(SAx)))) + ) + def zeeman_hamiltonian( self, B0: float, theta: Optional[float] = None, phi: Optional[float] = None ) -> np.ndarray: @@ -455,7 +473,7 @@ def hyperfine_hamiltonian(self, hfc_anisotropy: bool = False) -> np.ndarray: ) ) - def exchange_hamiltonian(self, J: float) -> np.ndarray: + def exchange_hamiltonian(self, J: float, prod_coeff: float = 2) -> np.ndarray: """Construct the exchange Hamiltonian. Construct the exchange (J-coupling) Hamiltonian based on the @@ -470,6 +488,10 @@ def exchange_hamiltonian(self, J: float) -> np.ndarray: J (float): Exchange coupling constant. + prod_coeff (float): Coefficient of the product operator + (default, radical-pair convention uses 2.0, + spintronics convention uses 1.0). + Returns: np.ndarray: @@ -480,7 +502,7 @@ def exchange_hamiltonian(self, J: float) -> np.ndarray: """ Jcoupling = J * self.radicals[0].gamma_mT SASB = self.product_operator(0, 1) - return Jcoupling * (2 * SASB + 0.5 * np.eye(SASB.shape[0])) + return Jcoupling * (prod_coeff * SASB + 0.5 * np.eye(SASB.shape[0])) def dipolar_hamiltonian(self, D: float | np.ndarray) -> np.ndarray: """Construct the Dipolar Hamiltonian. @@ -571,6 +593,20 @@ def dipolar_hamiltonian_3d(self, dipolar_tensor: np.ndarray) -> np.ndarray: ) ) + def zero_field_splitting_hamiltonian(self, D, E) -> np.ndarray: + """Construct the Zero Field Splitting (ZFS) Hamiltonian.""" + Dmod = D * -self.radicals[0].gamma_mT + Emod = E * -self.radicals[0].gamma_mT + result = complex(0.0) + for idx, p in enumerate(self.particles): + Sx = self.spin_operator(idx, "x") + Sy = self.spin_operator(idx, "y") + Sz = self.spin_operator(idx, "z") + Ssquared = Sx @ Sx + Sy @ Sy + Sz @ Sz + result += Dmod * (Sz @ Sz - (1 / 3) * Ssquared) + result += Emod * ((Sx @ Sx) - (Sy @ Sy)) + return result + def total_hamiltonian( self, B0: float, @@ -670,8 +706,8 @@ def product_probability(self, obs: State, rhos: np.ndarray) -> np.ndarray: """Calculate the probability of the observable from the densities.""" if obs == State.EQUILIBRIUM: raise ValueError("Observable state should not be EQUILIBRIUM") - obs = self.observable_projection_operator(obs) - return np.real(np.trace(obs @ rhos, axis1=-2, axis2=-1)) + Q = self.observable_projection_operator(obs) + return np.real(np.trace(Q @ rhos, axis1=-2, axis2=-1)) @staticmethod def product_yield(product_probability, time, k): @@ -938,6 +974,10 @@ def anisotropy( def convert(H: np.ndarray) -> np.ndarray: return H + @staticmethod + def _convert(Q: np.ndarray) -> np.ndarray: + return Q + def initial_density_matrix(self, state: State, H: np.ndarray) -> np.ndarray: """Create an initial desity matrix. @@ -1043,7 +1083,17 @@ class LiouvilleSimulation(HilbertSimulation): def convert(H: np.ndarray) -> np.ndarray: """Convert the Hamiltonian from Hilbert to Liouville space.""" eye = np.eye(len(H)) - return 1j * (np.kron(H, eye) - np.kron(eye, H.T)) + tmp = np.kron(H, eye) - np.kron(eye, H.T) + return 1j * tmp + + @staticmethod + def _convert(Q: np.ndarray) -> np.ndarray: + eye = np.eye(len(Q)) + return np.kron(Q, eye) + np.kron(eye, Q) + + @property + def hamiltonian_size(self): + return super().hamiltonian_size ** 2 @staticmethod def _square_liouville_rhos(rhos): diff --git a/tests/test_simulation.py b/tests/test_simulation.py index 2a0d5f6..4ba2eb5 100644 --- a/tests/test_simulation.py +++ b/tests/test_simulation.py @@ -7,6 +7,7 @@ import matplotlib.pyplot as plt import numpy as np + import radicalpy as rp from radicalpy import estimations, kinetics, relaxation from radicalpy.simulation import Basis @@ -301,6 +302,8 @@ def test_dipolar_interaction_1d(self): def test_initial_density_matrix(self): H = self.sim.total_hamiltonian(PARAMS["B"][0], PARAMS["J"], PARAMS["D"]) for state in rp.simulation.State: + if state == rp.simulation.State.TP_SINGLET: + continue rho0 = self.sim.initial_density_matrix(state, H) rpstate = state2radpy(state) rho0_true = radpy.Hilbert_initial(rpstate, len(self.sim.particles), H) @@ -318,9 +321,13 @@ def test_time_evolution(self): H = self.sim.total_hamiltonian(PARAMS["B"][0], PARAMS["J"], PARAMS["D"]) Kexp = kinetics.Exponential(k) for init_state in rp.simulation.State: + if init_state == rp.simulation.State.TP_SINGLET: + continue for obs_state in rp.simulation.State: if obs_state == rp.simulation.State.EQUILIBRIUM: continue + if obs_state == rp.simulation.State.TP_SINGLET: + continue evol_true = radpy.TimeEvolution( len(self.sim.particles), state2radpy(init_state), @@ -500,6 +507,8 @@ def setUp(self): def test_initial_density_matrix(self): H = self.sim.total_hamiltonian(PARAMS["B"][0], PARAMS["J"], PARAMS["D"]) for state in rp.simulation.State: + if state == rp.simulation.State.TP_SINGLET: + continue rho0 = self.sim.initial_density_matrix(state, H) rpstate = state2radpy(state) rho0_true = radpy.Liouville_initial(rpstate, len(self.sim.particles), H)