From 4f65a4a2fc141ca0b001dfda675f56c6b2e15e2b Mon Sep 17 00:00:00 2001 From: Emil VATAI Date: Wed, 25 Oct 2023 08:15:39 +0900 Subject: [PATCH 01/30] Start sctep implementation --- examples/sctep.py | 11 +++++++++++ radicalpy/simulation.py | 14 ++++++++++---- 2 files changed, 21 insertions(+), 4 deletions(-) create mode 100644 examples/sctep.py diff --git a/examples/sctep.py b/examples/sctep.py new file mode 100644 index 0000000..e55f655 --- /dev/null +++ b/examples/sctep.py @@ -0,0 +1,11 @@ +#! /usr/bin/env python + +import radicalpy as rp + + +def main(): + print(rp) + + +if __name__ == "__main__": + main() diff --git a/radicalpy/simulation.py b/radicalpy/simulation.py index 2b305e9..7d56bd5 100644 --- a/radicalpy/simulation.py +++ b/radicalpy/simulation.py @@ -206,10 +206,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: @@ -571,6 +571,12 @@ def dipolar_hamiltonian_3d(self, dipolar_tensor: np.ndarray) -> np.ndarray: ) ) + def zero_field_splitting_hamiltonian(self) -> np.ndarray: + """Construct the Zero Field Splitting (ZFS) Hamiltoninan.""" + 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)) + def total_hamiltonian( self, B0: float, From 8c9496f5bcab08060daba90031db4f4d3685ea7e Mon Sep 17 00:00:00 2001 From: Emil VATAI Date: Wed, 25 Oct 2023 10:17:01 +0900 Subject: [PATCH 02/30] WIP --- examples/sctep.py | 5 +++-- radicalpy/simulation.py | 19 +++++++++++++++---- 2 files changed, 18 insertions(+), 6 deletions(-) diff --git a/examples/sctep.py b/examples/sctep.py index e55f655..404d26b 100644 --- a/examples/sctep.py +++ b/examples/sctep.py @@ -1,10 +1,11 @@ #! /usr/bin/env python -import radicalpy as rp +from radicalpy.data import Molecule def main(): - print(rp) + m = Molecule(nuclei=[]) + print(m) if __name__ == "__main__": diff --git a/radicalpy/simulation.py b/radicalpy/simulation.py index 7d56bd5..3dec627 100644 --- a/radicalpy/simulation.py +++ b/radicalpy/simulation.py @@ -571,11 +571,22 @@ def dipolar_hamiltonian_3d(self, dipolar_tensor: np.ndarray) -> np.ndarray: ) ) - def zero_field_splitting_hamiltonian(self) -> np.ndarray: + def zero_field_splitting_hamiltonian(self, E, D) -> np.ndarray: """Construct the Zero Field Splitting (ZFS) Hamiltoninan.""" - 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)) + 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") + ) + ) + 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 def total_hamiltonian( self, From dd4a1004494b4fe8666a2df19859088893117f27 Mon Sep 17 00:00:00 2001 From: Emil VATAI Date: Thu, 23 Nov 2023 18:41:43 +0900 Subject: [PATCH 03/30] 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, From dc3e99b66db3511624fb7eee2b165a8b40ddf883 Mon Sep 17 00:00:00 2001 From: Emil VATAI Date: Fri, 24 Nov 2023 03:27:01 +0900 Subject: [PATCH 04/30] Add experiments.py --- radicalpy/experiments.py | 2 ++ 1 file changed, 2 insertions(+) create mode 100644 radicalpy/experiments.py diff --git a/radicalpy/experiments.py b/radicalpy/experiments.py new file mode 100644 index 0000000..9449e4d --- /dev/null +++ b/radicalpy/experiments.py @@ -0,0 +1,2 @@ +def steady_state_mary(): + pass From 12653c5a237cec62c13662c0d5a1842c55cfe661 Mon Sep 17 00:00:00 2001 From: Emil VATAI Date: Thu, 30 Nov 2023 19:18:23 +0900 Subject: [PATCH 05/30] Add Triplet class --- radicalpy/data.py | 7 +++++++ 1 file changed, 7 insertions(+) diff --git a/radicalpy/data.py b/radicalpy/data.py index a036acb..f7e9774 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) + super().__init__(name="Triplet", nuclei=[], radical=triplet) From 614736115cb90fd775f8f0744bbe15801226739b Mon Sep 17 00:00:00 2001 From: Emil VATAI Date: Thu, 7 Dec 2023 19:05:10 +0900 Subject: [PATCH 06/30] Add missing `.` --- radicalpy/data.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/radicalpy/data.py b/radicalpy/data.py index f7e9774..26d1d6d 100644 --- a/radicalpy/data.py +++ b/radicalpy/data.py @@ -520,5 +520,5 @@ def effective_hyperfine(self) -> float: class Triplet(Molecule): def __init__(self): gamma = Isotope("E").gamma_mT - triplet = Nucleus(magnetogyric_ratio=gamma, multiplicity=3, hfc=0) + triplet = Nucleus(magnetogyric_ratio=gamma, multiplicity=3, hfc=0.0) super().__init__(name="Triplet", nuclei=[], radical=triplet) From 01bbbe9a448820992c9195bc9857ae46bdce63b4 Mon Sep 17 00:00:00 2001 From: Emil VATAI Date: Thu, 7 Dec 2023 19:13:50 +0900 Subject: [PATCH 07/30] WIP --- examples/sctep.py | 42 ++++++++++++++++++---------------------- radicalpy/experiments.py | 8 ++++++-- 2 files changed, 25 insertions(+), 25 deletions(-) diff --git a/examples/sctep.py b/examples/sctep.py index 05cf229..ea9f8fb 100644 --- a/examples/sctep.py +++ b/examples/sctep.py @@ -1,6 +1,6 @@ #! /usr/bin/env python -from radicalpy.data import Isotope, Molecule, Nucleus +from radicalpy.data import Isotope, Molecule, Nucleus, Triplet from radicalpy.experiments import steady_state_mary from radicalpy.simulation import Basis, LiouvilleSimulation @@ -9,28 +9,24 @@ def main(): 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) + m = Triplet() + sim = LiouvilleSimulation(molecules=[m, m]) + print(sim) + # 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/experiments.py b/radicalpy/experiments.py index 9449e4d..25f481a 100644 --- a/radicalpy/experiments.py +++ b/radicalpy/experiments.py @@ -1,2 +1,6 @@ -def steady_state_mary(): - pass +#!/usr/bin/env python + + +def steady_state_mary(sim, D, E): + H = sim.zero_field_splitting_hamiltonian(D, E) + return H From 552edf905942441800913aff33a27bea2af304d1 Mon Sep 17 00:00:00 2001 From: Emil VATAI Date: Thu, 7 Dec 2023 19:51:04 +0900 Subject: [PATCH 08/30] Add commented code for new state --- examples/sctep.py | 25 +++++++++++-------------- radicalpy/simulation.py | 10 ++++++++++ 2 files changed, 21 insertions(+), 14 deletions(-) diff --git a/examples/sctep.py b/examples/sctep.py index ea9f8fb..4b5fd83 100644 --- a/examples/sctep.py +++ b/examples/sctep.py @@ -12,20 +12,17 @@ def main(): m = Triplet() sim = LiouvilleSimulation(molecules=[m, m]) print(sim) - # 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), - # ], - # ) + steady_state_mary( + sim, + obs_state=State.TP_SINGLET, + B=Bs, + D=D, + E=E, + kinetics=[ + rp.kinetics.Haberkorn(krec, State.TP_SINGLET), + rp.kinetics.HaberkornFree(kesc), + ], + ) # print(H) diff --git a/radicalpy/simulation.py b/radicalpy/simulation.py index 01c432f..2bcfb6a 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 = "RTP_S" class Basis(enum.Enum): @@ -306,6 +307,15 @@ def projection_operator(self, state: State): + (2 * SAz**2 - SAz) * (2 * SBz**2 - SBz), State.EQUILIBRIUM: 1.05459e-34 / (1.38e-23 * 298), } + # state.TP_SINGLET + # (1 / 12) * (Ssquared - (6 * np.eye(len(S1x)))) @ (Ssquared - (2 * np.eye(len(S1x)))) + # # For radical triplet pair (RTP) + # S1x, S1y, S1z = spinops(0, 2, spin=3) + # S2x, S2y, S2z = spinops(1, 2, spin=3) + # S1squared = S1x @ S1x + S1y @ S1y + S1z @ S1z + # S2squared = S2x @ S2x + S2y @ S2y + S2z @ S2z + # Ssquared = S1squared + S2squared + 2 * (S1x @ S2x + S1y @ S2y + S1z @ S2z) # + return result[state] def zeeman_hamiltonian( From fa5084bf319b5cb97956ab8b44727d9c16aa95d6 Mon Sep 17 00:00:00 2001 From: Emil VATAI Date: Thu, 11 Jan 2024 18:30:46 +0900 Subject: [PATCH 09/30] WIP --- examples/sctep.py | 30 ++++++++++++++++++++---------- radicalpy/experiments.py | 39 ++++++++++++++++++++++++++++++++++++--- radicalpy/kinetics.py | 23 ++++++++++++++++------- radicalpy/simulation.py | 22 +++++++++++++--------- 4 files changed, 85 insertions(+), 29 deletions(-) diff --git a/examples/sctep.py b/examples/sctep.py index 4b5fd83..ecf4477 100644 --- a/examples/sctep.py +++ b/examples/sctep.py @@ -1,27 +1,37 @@ #! /usr/bin/env python +import numpy as np from radicalpy.data import Isotope, Molecule, Nucleus, Triplet from radicalpy.experiments import steady_state_mary -from radicalpy.simulation import Basis, LiouvilleSimulation +from radicalpy.kinetics import Haberkorn, HaberkornFree +from radicalpy.simulation import Basis, HilbertSimulation, State def main(): - gamma = Isotope("E").gamma_mT - D = -6.2 * gamma - E = 35 * gamma + # gamma = Isotope("E").gamma_mT + Bs = np.arange(0, 2500, 1) + D = -6.2 # * gamma + E = 35 # * gamma + J = 499.55 + k0 = 1 + ks = 1.1e9 + kd = 2.8e9 m = Triplet() - sim = LiouvilleSimulation(molecules=[m, m]) + sim = HilbertSimulation(molecules=[m, m], basis=Basis.ZEEMAN) print(sim) + H = sim.zero_field_splitting_hamiltonian( + D=D, + E=E, + ) + print(H) steady_state_mary( sim, obs_state=State.TP_SINGLET, - B=Bs, + Bs=Bs, D=D, E=E, - kinetics=[ - rp.kinetics.Haberkorn(krec, State.TP_SINGLET), - rp.kinetics.HaberkornFree(kesc), - ], + J=J, + kinetics=[Haberkorn(ks, State.TP_SINGLET), HaberkornFree(kd)], ) # print(H) diff --git a/radicalpy/experiments.py b/radicalpy/experiments.py index 25f481a..406ae42 100644 --- a/radicalpy/experiments.py +++ b/radicalpy/experiments.py @@ -1,6 +1,39 @@ #!/usr/bin/env python +import numpy as np +from tqdm import tqdm -def steady_state_mary(sim, D, E): - H = sim.zero_field_splitting_hamiltonian(D, E) - return H +from .simulation import (HilbertIncoherentProcessBase, HilbertSimulation, + LiouvilleSimulation, State) + + +def steady_state_mary( + sim: HilbertSimulation, + obs_state: State, + Bs: np.ndarray, + D: float, + E: float, + J: float, + kinetics: list[HilbertIncoherentProcessBase] = [], + relaxations: list[HilbertIncoherentProcessBase] = [], +) -> np.ndarray: + theta = np.pi / 4 + phi = 0 + HZFS = sim.zero_field_splitting_hamiltonian(D, E) + HJ = sim.exchange_hamiltonian(J) + Phi_s = np.zeros((len(Bs)), dtype=complex) + Ps = sim.projection_operator(State.TP_SINGLET) + for i, B in enumerate(tqdm(Bs)): + HZ = sim.zeeman_hamiltonian_3d(B, theta, phi) + print(f"{HZ.shape=}") + print(f"{HZFS.shape=}") + print(f"{HJ.shape=}") + H = HZ + HZFS + HJ + sim.apply_liouville_hamiltonian_modifiers(H, kinetics + relaxations) + # HL = ry.Hilbert2Liouville(H) + # L = 1j * H + sum(kinetics) + rho = np.linalg.solve(H, Ps) # Density operator + + Phi_s[i] = np.matmul(Ps.T, rho) + + return rho, Phi_s diff --git a/radicalpy/kinetics.py b/radicalpy/kinetics.py index 4252a4b..a9cce4c 100644 --- a/radicalpy/kinetics.py +++ b/radicalpy/kinetics.py @@ -2,12 +2,9 @@ import numpy as np -from .simulation import ( - HilbertIncoherentProcessBase, - LiouvilleIncoherentProcessBase, - LiouvilleSimulation, - State, -) +from .simulation import (HilbertIncoherentProcessBase, + LiouvilleIncoherentProcessBase, LiouvilleSimulation, + State) class HilbertKineticsBase(HilbertIncoherentProcessBase): @@ -17,6 +14,10 @@ def _name(self): name = super()._name() return f"Kinetics: {name}" + @staticmethod + def _convert(Q: np.ndarray) -> np.ndarray: + return Q + class LiouvilleKineticsBase(LiouvilleIncoherentProcessBase): """Base class for kinetics superoperators (Liouville space).""" @@ -79,7 +80,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" ) diff --git a/radicalpy/simulation.py b/radicalpy/simulation.py index 2bcfb6a..39a677b 100644 --- a/radicalpy/simulation.py +++ b/radicalpy/simulation.py @@ -21,7 +21,7 @@ class State(enum.Enum): TRIPLET_PLUS = "T_+" TRIPLET_PLUS_MINUS = "T_\\pm" TRIPLET_MINUS = "T_-" - TP_SINGLET = "RTP_S" + TP_SINGLET = "TP_S" class Basis(enum.Enum): @@ -306,18 +306,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), } - # state.TP_SINGLET - # (1 / 12) * (Ssquared - (6 * np.eye(len(S1x)))) @ (Ssquared - (2 * np.eye(len(S1x)))) - # # For radical triplet pair (RTP) - # S1x, S1y, S1z = spinops(0, 2, spin=3) - # S2x, S2y, S2z = spinops(1, 2, spin=3) - # S1squared = S1x @ S1x + S1y @ S1y + S1z @ S1z - # S2squared = S2x @ S2x + S2y @ S2y + S2z @ S2z - # Ssquared = S1squared + S2squared + 2 * (S1x @ S2x + S1y @ S2y + S1z @ S2z) # 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: From 42683752c0854641a8f240c9be3c55de3d8fc46c Mon Sep 17 00:00:00 2001 From: Emil VATAI Date: Sun, 21 Jan 2024 12:36:36 +0900 Subject: [PATCH 10/30] Add isort config line length --- pyproject.toml | 9 ++++++--- 1 file changed, 6 insertions(+), 3 deletions(-) diff --git a/pyproject.toml b/pyproject.toml index 05f9fbc..19d205b 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 From f144ab32e422080aff39a1c3aea58e073d4233b9 Mon Sep 17 00:00:00 2001 From: Emil VATAI Date: Sun, 21 Jan 2024 12:41:51 +0900 Subject: [PATCH 11/30] Format conversion functions (while comparing them) --- radicalpy/kinetics.py | 4 +++- radicalpy/simulation.py | 3 ++- 2 files changed, 5 insertions(+), 2 deletions(-) diff --git a/radicalpy/kinetics.py b/radicalpy/kinetics.py index a9cce4c..16b5a7e 100644 --- a/radicalpy/kinetics.py +++ b/radicalpy/kinetics.py @@ -1,3 +1,4 @@ +"""Kinetics classes.""" from math import prod import numpy as np @@ -28,7 +29,8 @@ def _name(self): @staticmethod def _convert(Q: np.ndarray) -> np.ndarray: - return np.kron(Q, np.eye(len(Q))) + np.kron(np.eye(len(Q)), Q) + eye = np.eye(len(Q)) + return np.kron(Q, eye) + np.kron(eye, Q) class Exponential(HilbertKineticsBase): diff --git a/radicalpy/simulation.py b/radicalpy/simulation.py index 39a677b..b922bd2 100644 --- a/radicalpy/simulation.py +++ b/radicalpy/simulation.py @@ -1077,7 +1077,8 @@ 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 _square_liouville_rhos(rhos): From 79f555488ddd6ee4664bdd50c6fef9196a84e652 Mon Sep 17 00:00:00 2001 From: Emil VATAI Date: Sun, 21 Jan 2024 13:43:30 +0900 Subject: [PATCH 12/30] Move `_convert` method to simulation --- radicalpy/kinetics.py | 24 +++++++++--------------- radicalpy/simulation.py | 9 +++++++++ 2 files changed, 18 insertions(+), 15 deletions(-) diff --git a/radicalpy/kinetics.py b/radicalpy/kinetics.py index 16b5a7e..c364535 100644 --- a/radicalpy/kinetics.py +++ b/radicalpy/kinetics.py @@ -3,9 +3,12 @@ import numpy as np -from .simulation import (HilbertIncoherentProcessBase, - LiouvilleIncoherentProcessBase, LiouvilleSimulation, - State) +from .simulation import ( + HilbertIncoherentProcessBase, + LiouvilleIncoherentProcessBase, + LiouvilleSimulation, + State, +) class HilbertKineticsBase(HilbertIncoherentProcessBase): @@ -15,10 +18,6 @@ def _name(self): name = super()._name() return f"Kinetics: {name}" - @staticmethod - def _convert(Q: np.ndarray) -> np.ndarray: - return Q - class LiouvilleKineticsBase(LiouvilleIncoherentProcessBase): """Base class for kinetics superoperators (Liouville space).""" @@ -27,11 +26,6 @@ def _name(self): name = super()._name() return f"Kinetics: {name}" - @staticmethod - def _convert(Q: np.ndarray) -> np.ndarray: - eye = np.eye(len(Q)) - return np.kron(Q, eye) + np.kron(eye, Q) - class Exponential(HilbertKineticsBase): """Exponential model kinetics operator. @@ -98,7 +92,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 = [ @@ -154,8 +148,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 b922bd2..72c3006 100644 --- a/radicalpy/simulation.py +++ b/radicalpy/simulation.py @@ -972,6 +972,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. @@ -1080,6 +1084,11 @@ def convert(H: np.ndarray) -> np.ndarray: 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) + @staticmethod def _square_liouville_rhos(rhos): shape = rhos.shape From d866dcf2c2d08386eac1a655c1ff736bfba29ba5 Mon Sep 17 00:00:00 2001 From: Emil VATAI Date: Sun, 21 Jan 2024 14:01:42 +0900 Subject: [PATCH 13/30] Fix isort config in pyproject.toml --- pyproject.toml | 2 +- radicalpy/kinetics.py | 11 +++-------- 2 files changed, 4 insertions(+), 9 deletions(-) diff --git a/pyproject.toml b/pyproject.toml index 19d205b..baf9f2d 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -43,7 +43,7 @@ max-line-length = 88 disable = ["C0103"] [tool.isort] -line-length = 88 +line_length = 88 [tool.mypy] plugins = ["numpy.typing.mypy_plugin"] diff --git a/radicalpy/kinetics.py b/radicalpy/kinetics.py index c364535..e0ee3cd 100644 --- a/radicalpy/kinetics.py +++ b/radicalpy/kinetics.py @@ -1,14 +1,9 @@ """Kinetics classes.""" -from math import prod import numpy as np -from .simulation import ( - HilbertIncoherentProcessBase, - LiouvilleIncoherentProcessBase, - LiouvilleSimulation, - State, -) +from .simulation import (HilbertIncoherentProcessBase, LiouvilleIncoherentProcessBase, + LiouvilleSimulation, State) class HilbertKineticsBase(HilbertIncoherentProcessBase): @@ -117,7 +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 + size = np.prod(p.multiplicity for p in sim.particles) ** 2 self.subH = self.rate * np.eye(size) From 42604f5558d530d0f3c3bb31316e210b1b049843 Mon Sep 17 00:00:00 2001 From: Emil VATAI Date: Sun, 21 Jan 2024 17:15:45 +0900 Subject: [PATCH 14/30] Move Hamiltonian size calculation --- radicalpy/kinetics.py | 3 +-- radicalpy/simulation.py | 8 ++++++++ 2 files changed, 9 insertions(+), 2 deletions(-) diff --git a/radicalpy/kinetics.py b/radicalpy/kinetics.py index e0ee3cd..ea02f8f 100644 --- a/radicalpy/kinetics.py +++ b/radicalpy/kinetics.py @@ -112,8 +112,7 @@ class HaberkornFree(LiouvilleKineticsBase): def init(self, sim: LiouvilleSimulation): """See `radicalpy.simulation.HilbertIncoherentProcessBase.init`.""" - size = np.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): diff --git a/radicalpy/simulation.py b/radicalpy/simulation.py index 72c3006..f64af40 100644 --- a/radicalpy/simulation.py +++ b/radicalpy/simulation.py @@ -105,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( [ @@ -1089,6 +1093,10 @@ 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): shape = rhos.shape From 6221ee02ba589dffbb59de67e3d93af83c63a177 Mon Sep 17 00:00:00 2001 From: Emil VATAI Date: Sun, 21 Jan 2024 17:16:17 +0900 Subject: [PATCH 15/30] SCTEP looks good --- examples/sctep.py | 19 ++++++++++--------- radicalpy/experiments.py | 3 +++ 2 files changed, 13 insertions(+), 9 deletions(-) diff --git a/examples/sctep.py b/examples/sctep.py index ecf4477..43c0e10 100644 --- a/examples/sctep.py +++ b/examples/sctep.py @@ -1,21 +1,22 @@ #! /usr/bin/env python import numpy as np + from radicalpy.data import Isotope, Molecule, Nucleus, Triplet from radicalpy.experiments import steady_state_mary from radicalpy.kinetics import Haberkorn, HaberkornFree from radicalpy.simulation import Basis, HilbertSimulation, State -def main(): - # gamma = Isotope("E").gamma_mT - Bs = np.arange(0, 2500, 1) - D = -6.2 # * gamma - E = 35 # * gamma - J = 499.55 - k0 = 1 - ks = 1.1e9 - kd = 2.8e9 +def main( + Bs=np.arange(0, 2500, 1), + D=-6.2, + E=35, + J=499.55, + k0=1, + ks=1.1e9, + kd=2.8e9, +): m = Triplet() sim = HilbertSimulation(molecules=[m, m], basis=Basis.ZEEMAN) print(sim) diff --git a/radicalpy/experiments.py b/radicalpy/experiments.py index 406ae42..1ba7dc4 100644 --- a/radicalpy/experiments.py +++ b/radicalpy/experiments.py @@ -34,6 +34,9 @@ def steady_state_mary( # L = 1j * H + sum(kinetics) rho = np.linalg.solve(H, Ps) # Density operator + print(f"{type(rho)=} {rho.shape=}") + print(f"{type(Ps)=} {Ps.shape=}") + Phi_s[i] = np.matmul(Ps.T, rho) return rho, Phi_s From b68184b554c840de3cf890cfd73441dc2edfed4a Mon Sep 17 00:00:00 2001 From: Emil VATAI Date: Sun, 21 Jan 2024 18:39:54 +0900 Subject: [PATCH 16/30] It runs (but wrong output) --- examples/sctep.py | 25 +++++++++++++++++++++++-- radicalpy/experiments.py | 29 ++++++++++------------------- 2 files changed, 33 insertions(+), 21 deletions(-) diff --git a/examples/sctep.py b/examples/sctep.py index 43c0e10..0cd5178 100644 --- a/examples/sctep.py +++ b/examples/sctep.py @@ -1,5 +1,6 @@ #! /usr/bin/env python +import matplotlib.pyplot as plt import numpy as np from radicalpy.data import Isotope, Molecule, Nucleus, Triplet @@ -25,16 +26,36 @@ def main( E=E, ) print(H) - steady_state_mary( + rhos, Phi_s = steady_state_mary( sim, - obs_state=State.TP_SINGLET, + 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__": diff --git a/radicalpy/experiments.py b/radicalpy/experiments.py index 1ba7dc4..e44f93b 100644 --- a/radicalpy/experiments.py +++ b/radicalpy/experiments.py @@ -9,34 +9,25 @@ def steady_state_mary( sim: HilbertSimulation, - obs_state: State, + obs: State, Bs: np.ndarray, D: float, E: float, J: float, + theta: float, + phi: float, kinetics: list[HilbertIncoherentProcessBase] = [], - relaxations: list[HilbertIncoherentProcessBase] = [], + # relaxations: list[HilbertIncoherentProcessBase] = [], ) -> np.ndarray: - theta = np.pi / 4 - phi = 0 HZFS = sim.zero_field_splitting_hamiltonian(D, E) HJ = sim.exchange_hamiltonian(J) - Phi_s = np.zeros((len(Bs)), dtype=complex) - Ps = sim.projection_operator(State.TP_SINGLET) + rhos = np.zeros(shape=(len(Bs), *HJ.shape)) + Q = sim.projection_operator(obs) for i, B in enumerate(tqdm(Bs)): HZ = sim.zeeman_hamiltonian_3d(B, theta, phi) - print(f"{HZ.shape=}") - print(f"{HZFS.shape=}") - print(f"{HJ.shape=}") H = HZ + HZFS + HJ - sim.apply_liouville_hamiltonian_modifiers(H, kinetics + relaxations) - # HL = ry.Hilbert2Liouville(H) - # L = 1j * H + sum(kinetics) - rho = np.linalg.solve(H, Ps) # Density operator + sim.apply_liouville_hamiltonian_modifiers(H, kinetics) # + relaxations) + rhos[i] = np.linalg.solve(H, Q) + Phi_s = sim.product_probability(obs, rhos) - print(f"{type(rho)=} {rho.shape=}") - print(f"{type(Ps)=} {Ps.shape=}") - - Phi_s[i] = np.matmul(Ps.T, rho) - - return rho, Phi_s + return rhos, Phi_s From d3f01e4e76cf4c6fa5a715f858607ed0748552e1 Mon Sep 17 00:00:00 2001 From: Emil VATAI Date: Sun, 21 Jan 2024 22:09:33 +0900 Subject: [PATCH 17/30] Skip TP_SINGLET state in tests --- tests/test_simulation.py | 3 +++ 1 file changed, 3 insertions(+) diff --git a/tests/test_simulation.py b/tests/test_simulation.py index 2a0d5f6..2a88793 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 @@ -500,6 +501,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) From 8489835747b76f98aaa036593f7faea9015cb912 Mon Sep 17 00:00:00 2001 From: Emil VATAI Date: Sun, 21 Jan 2024 22:14:11 +0900 Subject: [PATCH 18/30] More TP_SINGLET skipping --- tests/test_simulation.py | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/tests/test_simulation.py b/tests/test_simulation.py index 2a88793..f9a3302 100644 --- a/tests/test_simulation.py +++ b/tests/test_simulation.py @@ -302,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) @@ -322,6 +324,8 @@ def test_time_evolution(self): for obs_state in rp.simulation.State: if obs_state == rp.simulation.State.EQUILIBRIUM: continue + if state == rp.simulation.State.TP_SINGLET: + continue evol_true = radpy.TimeEvolution( len(self.sim.particles), state2radpy(init_state), From 083da890df9010643c9b9d41f44fd376035753bc Mon Sep 17 00:00:00 2001 From: Emil VATAI Date: Sun, 21 Jan 2024 22:16:03 +0900 Subject: [PATCH 19/30] Fix typo --- tests/test_simulation.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/test_simulation.py b/tests/test_simulation.py index f9a3302..df5520e 100644 --- a/tests/test_simulation.py +++ b/tests/test_simulation.py @@ -324,7 +324,7 @@ def test_time_evolution(self): for obs_state in rp.simulation.State: if obs_state == rp.simulation.State.EQUILIBRIUM: continue - if state == rp.simulation.State.TP_SINGLET: + if obs_state == rp.simulation.State.TP_SINGLET: continue evol_true = radpy.TimeEvolution( len(self.sim.particles), From d86c670bda2c65d027470bfb07fa49fde9b1a981 Mon Sep 17 00:00:00 2001 From: Emil VATAI Date: Sun, 21 Jan 2024 22:17:58 +0900 Subject: [PATCH 20/30] More TP_SINGLET skipping --- tests/test_simulation.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/tests/test_simulation.py b/tests/test_simulation.py index df5520e..4ba2eb5 100644 --- a/tests/test_simulation.py +++ b/tests/test_simulation.py @@ -321,6 +321,8 @@ 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 From c6b2085b1792abd186e782e33c85965513183f1f Mon Sep 17 00:00:00 2001 From: Emil VATAI Date: Sun, 21 Jan 2024 23:54:04 +0900 Subject: [PATCH 21/30] Fix lint warning --- radicalpy/simulation.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/radicalpy/simulation.py b/radicalpy/simulation.py index f64af40..e37d698 100644 --- a/radicalpy/simulation.py +++ b/radicalpy/simulation.py @@ -708,8 +708,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): From 738b236d8ca12b214c5ef427e30005217e84b8b0 Mon Sep 17 00:00:00 2001 From: Emil VATAI Date: Mon, 22 Jan 2024 00:31:33 +0900 Subject: [PATCH 22/30] Disable modified ZFSH calculation --- radicalpy/simulation.py | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) diff --git a/radicalpy/simulation.py b/radicalpy/simulation.py index e37d698..9cf820b 100644 --- a/radicalpy/simulation.py +++ b/radicalpy/simulation.py @@ -598,15 +598,16 @@ def zero_field_splitting_hamiltonian(self, D, E) -> np.ndarray: 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 + # result += coeff * sqr tmp += coeff * sqr Sx = self.spin_operator(idx, "x") @ self.spin_operator(idx, "x") 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=}") + # other = D * (Sz @ Sz - (1 / 3) * Ssquared) + E * ((Sx @ Sx) - (Sy @ Sy)) + # delta = np.linalg.norm(other - sqr) + # print(f"{delta=}") + result += D * (Sz @ Sz - (1 / 3) * Ssquared) + E * ((Sx @ Sx) - (Sy @ Sy)) return result def total_hamiltonian( From 88a250351069bac9afa04122d3a29ea9cfad5656 Mon Sep 17 00:00:00 2001 From: Emil VATAI Date: Mon, 22 Jan 2024 00:31:08 +0900 Subject: [PATCH 23/30] Switch to Liouvillian --- examples/sctep.py | 12 +++--------- radicalpy/experiments.py | 15 +++++++++++---- 2 files changed, 14 insertions(+), 13 deletions(-) diff --git a/examples/sctep.py b/examples/sctep.py index 0cd5178..d129957 100644 --- a/examples/sctep.py +++ b/examples/sctep.py @@ -6,11 +6,11 @@ from radicalpy.data import Isotope, Molecule, Nucleus, Triplet from radicalpy.experiments import steady_state_mary from radicalpy.kinetics import Haberkorn, HaberkornFree -from radicalpy.simulation import Basis, HilbertSimulation, State +from radicalpy.simulation import Basis, LiouvilleSimulation, State def main( - Bs=np.arange(0, 2500, 1), + Bs=np.arange(0, 2500, 10), D=-6.2, E=35, J=499.55, @@ -19,13 +19,7 @@ def main( kd=2.8e9, ): m = Triplet() - sim = HilbertSimulation(molecules=[m, m], basis=Basis.ZEEMAN) - print(sim) - H = sim.zero_field_splitting_hamiltonian( - D=D, - E=E, - ) - print(H) + sim = LiouvilleSimulation(molecules=[m, m], basis=Basis.ZEEMAN) rhos, Phi_s = steady_state_mary( sim, obs=State.TP_SINGLET, diff --git a/radicalpy/experiments.py b/radicalpy/experiments.py index e44f93b..c719de5 100644 --- a/radicalpy/experiments.py +++ b/radicalpy/experiments.py @@ -8,7 +8,7 @@ def steady_state_mary( - sim: HilbertSimulation, + sim: LiouvilleSimulation, obs: State, Bs: np.ndarray, D: float, @@ -21,13 +21,20 @@ def steady_state_mary( ) -> np.ndarray: HZFS = sim.zero_field_splitting_hamiltonian(D, E) HJ = sim.exchange_hamiltonian(J) - rhos = np.zeros(shape=(len(Bs), *HJ.shape)) + 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_3d(B, theta, phi) H = HZ + HZFS + HJ + H = sim.convert(H) sim.apply_liouville_hamiltonian_modifiers(H, kinetics) # + relaxations) - rhos[i] = np.linalg.solve(H, Q) - Phi_s = sim.product_probability(obs, rhos) + rhos[i] = np.linalg.solve(H, Q.flatten()) + # Phi_s = sim.product_probability(obs, rhos) + print(f"{Q.shape=}") + print(f"{rhos.shape=}") + # Phi_s = np.sum(Q * rhos, axis=(-1, -2)) + Phi_s = Q.flatten() @ rhos.T + # print(f"{(Q * rhos).shape=}") + print(f"{Phi_s.shape=}") return rhos, Phi_s From 84b65e382a3ad6e4cc02cfcd158a93eb4db8e9de Mon Sep 17 00:00:00 2001 From: Emil VATAI Date: Mon, 22 Jan 2024 20:29:33 +0900 Subject: [PATCH 24/30] Remove unused imports --- examples/sctep.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/examples/sctep.py b/examples/sctep.py index d129957..2eb7bdb 100644 --- a/examples/sctep.py +++ b/examples/sctep.py @@ -3,7 +3,7 @@ import matplotlib.pyplot as plt import numpy as np -from radicalpy.data import Isotope, Molecule, Nucleus, Triplet +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 From 3c839d3d467216fe09197c0cd171133ee2895865 Mon Sep 17 00:00:00 2001 From: Emil VATAI Date: Mon, 22 Jan 2024 21:39:08 +0900 Subject: [PATCH 25/30] Check ZFS hamiltonian --- radicalpy/simulation.py | 16 +++------------- 1 file changed, 3 insertions(+), 13 deletions(-) diff --git a/radicalpy/simulation.py b/radicalpy/simulation.py index 9cf820b..101d3f8 100644 --- a/radicalpy/simulation.py +++ b/radicalpy/simulation.py @@ -591,22 +591,12 @@ 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) Hamiltoninan.""" - 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") - Sy = self.spin_operator(idx, "y") @ self.spin_operator(idx, "y") - Sz = self.spin_operator(idx, "z") @ self.spin_operator(idx, "z") + 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 - # other = D * (Sz @ Sz - (1 / 3) * Ssquared) + E * ((Sx @ Sx) - (Sy @ Sy)) - # delta = np.linalg.norm(other - sqr) - # print(f"{delta=}") result += D * (Sz @ Sz - (1 / 3) * Ssquared) + E * ((Sx @ Sx) - (Sy @ Sy)) return result From c99de107cb07c81e780e3ccf9d7e8aaa0773e750 Mon Sep 17 00:00:00 2001 From: Emil VATAI Date: Mon, 22 Jan 2024 22:14:21 +0900 Subject: [PATCH 26/30] Add `prod_coeff` to `exchange_hamiltonian` --- radicalpy/experiments.py | 4 ++-- radicalpy/simulation.py | 8 ++++++-- 2 files changed, 8 insertions(+), 4 deletions(-) diff --git a/radicalpy/experiments.py b/radicalpy/experiments.py index c719de5..aa7992c 100644 --- a/radicalpy/experiments.py +++ b/radicalpy/experiments.py @@ -19,8 +19,8 @@ def steady_state_mary( kinetics: list[HilbertIncoherentProcessBase] = [], # relaxations: list[HilbertIncoherentProcessBase] = [], ) -> np.ndarray: - HZFS = sim.zero_field_splitting_hamiltonian(D, E) - HJ = sim.exchange_hamiltonian(J) + 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)): diff --git a/radicalpy/simulation.py b/radicalpy/simulation.py index 101d3f8..e809e28 100644 --- a/radicalpy/simulation.py +++ b/radicalpy/simulation.py @@ -473,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 @@ -488,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: @@ -498,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. From 4a037d3421fa237fd05a6e9a85f076634fec845a Mon Sep 17 00:00:00 2001 From: Emil VATAI Date: Mon, 22 Jan 2024 22:15:14 +0900 Subject: [PATCH 27/30] Fix typo --- radicalpy/simulation.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/radicalpy/simulation.py b/radicalpy/simulation.py index e809e28..8c9d8d4 100644 --- a/radicalpy/simulation.py +++ b/radicalpy/simulation.py @@ -594,7 +594,7 @@ 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) Hamiltoninan.""" + """Construct the Zero Field Splitting (ZFS) Hamiltonian.""" result = complex(0.0) for idx, p in enumerate(self.particles): Sx = self.spin_operator(idx, "x") From 2a19910e5ca2e9b743b0f70cf5adcd647a80df5b Mon Sep 17 00:00:00 2001 From: Emil VATAI Date: Wed, 24 Jan 2024 22:21:00 +0900 Subject: [PATCH 28/30] Fix sctep.py --- radicalpy/experiments.py | 6 ++++-- radicalpy/simulation.py | 5 ++++- 2 files changed, 8 insertions(+), 3 deletions(-) diff --git a/radicalpy/experiments.py b/radicalpy/experiments.py index aa7992c..4f28cad 100644 --- a/radicalpy/experiments.py +++ b/radicalpy/experiments.py @@ -24,16 +24,18 @@ def steady_state_mary( 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_3d(B, theta, phi) + 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 = sim.product_probability(obs, rhos) print(f"{Q.shape=}") print(f"{rhos.shape=}") # Phi_s = np.sum(Q * rhos, axis=(-1, -2)) - Phi_s = Q.flatten() @ rhos.T + print(f"{Q.shape=} {rhos.shape=}") + Phi_s = rhos @ Q.flatten() # print(f"{(Q * rhos).shape=}") print(f"{Phi_s.shape=}") diff --git a/radicalpy/simulation.py b/radicalpy/simulation.py index 8c9d8d4..85adadd 100644 --- a/radicalpy/simulation.py +++ b/radicalpy/simulation.py @@ -595,13 +595,16 @@ 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 += D * (Sz @ Sz - (1 / 3) * Ssquared) + E * ((Sx @ Sx) - (Sy @ Sy)) + result += Dmod * (Sz @ Sz - (1 / 3) * Ssquared) + result += Emod * ((Sx @ Sx) - (Sy @ Sy)) return result def total_hamiltonian( From 232327fc0ca73662c0eb3b4b7d9937042009738d Mon Sep 17 00:00:00 2001 From: Emil VATAI Date: Wed, 24 Jan 2024 22:22:10 +0900 Subject: [PATCH 29/30] Cleanup --- radicalpy/experiments.py | 8 -------- 1 file changed, 8 deletions(-) diff --git a/radicalpy/experiments.py b/radicalpy/experiments.py index 4f28cad..ea82e6c 100644 --- a/radicalpy/experiments.py +++ b/radicalpy/experiments.py @@ -30,13 +30,5 @@ def steady_state_mary( sim.apply_liouville_hamiltonian_modifiers(H, kinetics) # + relaxations) rhos[i] = np.linalg.solve(H, Q.flatten()) - # Phi_s = sim.product_probability(obs, rhos) - print(f"{Q.shape=}") - print(f"{rhos.shape=}") - # Phi_s = np.sum(Q * rhos, axis=(-1, -2)) - print(f"{Q.shape=} {rhos.shape=}") Phi_s = rhos @ Q.flatten() - # print(f"{(Q * rhos).shape=}") - print(f"{Phi_s.shape=}") - return rhos, Phi_s From b52b33d60ecb2a63d7bd528b79c5579c832ad183 Mon Sep 17 00:00:00 2001 From: Emil VATAI Date: Wed, 24 Jan 2024 22:59:45 +0900 Subject: [PATCH 30/30] Move minuses inside ZFS Hamiltonian function --- radicalpy/experiments.py | 5 ++--- radicalpy/simulation.py | 4 ++-- 2 files changed, 4 insertions(+), 5 deletions(-) diff --git a/radicalpy/experiments.py b/radicalpy/experiments.py index ea82e6c..a3a56cc 100644 --- a/radicalpy/experiments.py +++ b/radicalpy/experiments.py @@ -3,8 +3,7 @@ import numpy as np from tqdm import tqdm -from .simulation import (HilbertIncoherentProcessBase, HilbertSimulation, - LiouvilleSimulation, State) +from .simulation import HilbertIncoherentProcessBase, LiouvilleSimulation, State def steady_state_mary( @@ -19,7 +18,7 @@ def steady_state_mary( kinetics: list[HilbertIncoherentProcessBase] = [], # relaxations: list[HilbertIncoherentProcessBase] = [], ) -> np.ndarray: - HZFS = sim.zero_field_splitting_hamiltonian(-D, -E) + 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) diff --git a/radicalpy/simulation.py b/radicalpy/simulation.py index 85adadd..e432ffc 100644 --- a/radicalpy/simulation.py +++ b/radicalpy/simulation.py @@ -595,8 +595,8 @@ 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 + 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")