Skip to content

Commit

Permalink
Merge pull request #142 from Spin-Chemistry-Labs/139-implement-sctep
Browse files Browse the repository at this point in the history
139 implement sctep
  • Loading branch information
vatai authored Jan 24, 2024
2 parents c96150d + b52b33d commit 39facdc
Show file tree
Hide file tree
Showing 7 changed files with 186 additions and 29 deletions.
56 changes: 56 additions & 0 deletions examples/sctep.py
Original file line number Diff line number Diff line change
@@ -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()
9 changes: 6 additions & 3 deletions pyproject.toml
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
7 changes: 7 additions & 0 deletions radicalpy/data.py
Original file line number Diff line number Diff line change
Expand Up @@ -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)
33 changes: 33 additions & 0 deletions radicalpy/experiments.py
Original file line number Diff line number Diff line change
@@ -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
33 changes: 16 additions & 17 deletions radicalpy/kinetics.py
Original file line number Diff line number Diff line change
@@ -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):
Expand All @@ -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.
Expand Down Expand Up @@ -79,15 +71,23 @@ 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"
)

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 = [
Expand All @@ -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):
Expand Down Expand Up @@ -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))
)
Expand Down
68 changes: 59 additions & 9 deletions radicalpy/simulation.py
Original file line number Diff line number Diff line change
Expand Up @@ -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):
Expand Down Expand Up @@ -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(
[
Expand Down Expand Up @@ -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:
Expand Down Expand Up @@ -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:
Expand Down Expand Up @@ -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
Expand All @@ -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:
Expand All @@ -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.
Expand Down Expand Up @@ -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,
Expand Down Expand Up @@ -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):
Expand Down Expand Up @@ -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.
Expand Down Expand Up @@ -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):
Expand Down
9 changes: 9 additions & 0 deletions tests/test_simulation.py
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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)
Expand All @@ -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),
Expand Down Expand Up @@ -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)
Expand Down

0 comments on commit 39facdc

Please sign in to comment.