From dd4a1004494b4fe8666a2df19859088893117f27 Mon Sep 17 00:00:00 2001 From: Emil VATAI Date: Thu, 23 Nov 2023 18:41:43 +0900 Subject: [PATCH] Add zero field splitting Hamiltonian --- examples/sctep.py | 31 ++++++++++++++++++++++++++++--- radicalpy/simulation.py | 31 +++++++++++++++++-------------- 2 files changed, 45 insertions(+), 17 deletions(-) diff --git a/examples/sctep.py b/examples/sctep.py index 404d26b..05cf229 100644 --- a/examples/sctep.py +++ b/examples/sctep.py @@ -1,11 +1,36 @@ #! /usr/bin/env python -from radicalpy.data import Molecule +from radicalpy.data import Isotope, Molecule, Nucleus +from radicalpy.experiments import steady_state_mary +from radicalpy.simulation import Basis, LiouvilleSimulation def main(): - m = Molecule(nuclei=[]) - print(m) + gamma = Isotope("E").gamma_mT + D = -6.2 * gamma + E = 35 * gamma + triplet = Nucleus( + magnetogyric_ratio=gamma, + multiplicity=3, + hfc=0.0, + name="Triplet", + ) + m = Molecule(radical=triplet) + sim = LiouvilleSimulation(molecules=[m, m], basis=Basis.ZEEMAN) + # print(sim) + steady_state_mary( + sim, + obs_state=State.TP_SINGLET, + B=Bs, + D=D, + E=E, + J=0, + kinetics=[ + rp.kinetics.Haberkorn(krec, State.TP_SINGLET), + rp.kinetics.HaberkornFree(kesc), + ], + ) + print(H) if __name__ == "__main__": diff --git a/radicalpy/simulation.py b/radicalpy/simulation.py index 3dec627..01c432f 100644 --- a/radicalpy/simulation.py +++ b/radicalpy/simulation.py @@ -571,22 +571,25 @@ def dipolar_hamiltonian_3d(self, dipolar_tensor: np.ndarray) -> np.ndarray: ) ) - def zero_field_splitting_hamiltonian(self, E, D) -> np.ndarray: + def zero_field_splitting_hamiltonian(self, D, E) -> np.ndarray: """Construct the Zero Field Splitting (ZFS) Hamiltoninan.""" - for idx in range(self.particles): - spinops = [] - coeffs = [E - D / 3, -E - D / 3, 2 * D / 3] - sum( - coeff * np.linalg.matrix_power(self.spin_operator(idx, ax), 2) - for idx, (coeff, ax) in itertools.product( - self.particles, zip(coeffs, "xyz") - ) - ) + coeffs = [E - D / 3, -E - D / 3, 2 * D / 3] + result = complex(0.0) + other = complex(0.0) + for idx, p in enumerate(self.particles): + tmp = complex(0.0) + for coeff, ax in zip(coeffs, "xyz"): + sqr = self.spin_operator(idx, ax) @ self.spin_operator(idx, ax) + result += coeff * sqr + tmp += coeff * sqr Sx = self.spin_operator(idx, "x") @ self.spin_operator(idx, "x") - # Sx, Sy, Sz = spinops(pos, 2, spin=3) - # Ssquared = Sx @ Sx + Sy @ Sy + Sz @ Sz - # return D * (Sz @ Sz - (1 / 3) * Ssquared) + E * ((Sx @ Sx) - (Sy @ Sy)) - return 0 + Sy = self.spin_operator(idx, "y") @ self.spin_operator(idx, "y") + Sz = self.spin_operator(idx, "z") @ self.spin_operator(idx, "z") + Ssquared = Sx @ Sx + Sy @ Sy + Sz @ Sz + other = D * (Sz @ Sz - (1 / 3) * Ssquared) + E * ((Sx @ Sx) - (Sy @ Sy)) + delta = np.linalg.norm(other - sqr) + print(f"{delta=}") + return result def total_hamiltonian( self,