Merge pull request #206 from Spin-Chemistry-Labs/205-add-soc-and-nqi-hamiltonians

205 add soc and nqi hamiltonians

Author: lmantill

examples/bloch-redfield.py

@@ -14,46 +14,37 @@
 def main():
     m1 = Molecule.fromdb("flavin_anion", ["N5"])
     m2 = Molecule.fromdb("tryptophan_cation", [])
-    simH = HilbertSimulation([m1, m2])
-    sim = LiouvilleSimulation([m1, m2])
-    H = simH.total_hamiltonian(B0=0, J=0, D=0)
+    sim = HilbertSimulation([m1, m2])
+    H = sim.total_hamiltonian(B0=0, J=0, D=0)
 
-    sigmax = simH.spin_operator(0, "x") + simH.spin_operator(1, "x")
+    rho0 = sim.projection_operator(State.SINGLET)
+    obs = sim.projection_operator(State.SINGLET)
 
-    S = 3e6
-    BRs = sigmax
-    blochredfield = sim.bloch_redfield_liouvillian(
-        H, channels=[(BRs, S)], secular=True, secular_cutoff=0.01
-    )
+    sigmax = sim.spin_operator(0, "x") + sim.spin_operator(1, "x")
+    sigmay = sim.spin_operator(0, "y") + sim.spin_operator(1, "y")
+    sigmaz = sim.spin_operator(0, "z") + sim.spin_operator(1, "z")
 
-    dB = 0.5
-    B0 = np.arange(0.0, 20.0 + 1e-9, dB)
-    time = np.arange(0, 3e-6, 10e-9)
-    kinetics = [Haberkorn(3e6, State.SINGLET), HaberkornFree(1e6)]
-    relaxations = [SingletTripletDephasing(1e7)]
+    Sxyz = sigmax + sigmay + sigmaz
 
-    init_state = State.SINGLET
-    obs_state = State.SINGLET
+    def S(omega: float, tau_c: float) -> float:
+        omega = float(omega)
+        return tau_c / (1.0 + (omega * tau_c) * (omega * tau_c))
 
-    sim.apply_liouville_hamiltonian_modifiers(blochredfield, kinetics + relaxations)
-    rhos = magnetic_field_loop(sim, init_state, time, blochredfield, B0, B_axis="z")
-    product_probabilities = sim.product_probability(obs_state, rhos)
+    tau_c = 1e7
+    NPS = lambda omega: S(omega, tau_c)
 
-    sim.apply_hilbert_kinetics(time, product_probabilities, kinetics)
-    k = kinetics[0].rate_constant if kinetics else 1.0
-    product_yields, product_yield_sums = sim.product_yield(
-        product_probabilities, time, k
-    )
+    time = np.arange(0, 5e-6, 1e-9)
 
-    x = time * 1e6
-    n = 0
+    L = sim.bloch_redfield_time_evolution(
+        H, rho0, time, bath=[Sxyz], noise=[NPS], obs=[obs]
+    )
+    L_result = L["expect"][0] / np.trace(rho0)
 
     fig = plt.figure(1)
-    plt.plot(x, product_probabilities[n, :], linewidth=2, label=r"$P_i(t)$")
-    plt.fill_between(x, product_yields[n, :], alpha=0.2, label=r"$\Phi_i$")
+    plt.plot(time * 1e6, L_result)
     plt.xlabel("Time / $\mu s$", size=14)
     plt.ylabel("Probability", size=14)
-    # plt.ylim([0, 1])
+    plt.ylim([0, 1])
     plt.legend(fontsize=14)
     fig.set_size_inches([7, 4])
     plt.show()

examples/energy_level.py

@@ -8,18 +8,14 @@
 
 
 def main():
-    flavin = rp.simulation.Molecule.fromdb("flavin_anion", ["H25"])
+    flavin = rp.simulation.Molecule.fromdb("flavin_anion", [])
     Z = rp.simulation.Molecule("Z")
     sim = rp.simulation.HilbertSimulation([flavin, Z])
-    H = sim.total_hamiltonian(B0=1, D=0, J=0)
 
-    eigval = np.linalg.eigh(H)
-    E = np.real(eigval[0])  # 0 = eigenvalues, 1 = eigenvectors
-
-    rp.plot.energy_levels(sim, B=np.arange(0.01, 1, 0.01), J=0, D=0)
-    # plt.show()
-    path = __file__[:-3] + f"_{0}.png"
-    plt.savefig(path)
+    rp.plot.energy_levels(sim, B=np.arange(0.01, 1, 0.01), J=-0.05, D=0)
+    plt.show()
+    # path = __file__[:-3] + f"_{0}.png"
+    # plt.savefig(path)
 
 
 if __name__ == "__main__":

radicalpy/classical.py

@@ -7,58 +7,55 @@
 sampling.
 
 
-Contents
---------
+Contents:
 
-- Diffusion:
+    - Diffusion:
 
-  - get_delta_r(D, dt): RMS relative displacement √(6 D Δt).
+    - get_delta_r(D, dt): RMS relative displacement √(6 D Δt).
 
-- LaTeX utilities:
+    - LaTeX utilities:
 
-  - latexify(rate_equations): Build d[Xi]/dt = … strings from a rate map.
+    - latexify(rate_equations): Build d[Xi]/dt = … strings from a rate map.
 
-  - latex_eqlist_to_align(eqs): Wrap equations in an `align*` environment.
+    - latex_eqlist_to_align(eqs): Wrap equations in an `align*` environment.
 
-- Kinetics:
+    - Kinetics:
 
-  - Rate: Numeric value + LaTeX label with operator overloading that carries
-    symbolic expressions through +, −, ×, ÷.
+    - Rate: Numeric value + LaTeX label with operator overloading that carries
+        symbolic expressions through +, −, ×, ÷.
 
-  - RateEquations: Sparse representation and time evolution (matrix exponential)
-    for first-order networks; returns EquationRateResult for easy access.
+    - RateEquations: Sparse representation and time evolution (matrix exponential)
+        for first-order networks; returns EquationRateResult for easy access.
 
-  - EquationRateResult: Convenience accessor to sum populations of one or more
-    states over time.
+    - EquationRateResult: Convenience accessor to sum populations of one or more
+        states over time.
 
-- Visualization:
+    - Visualization:
 
-  - reaction_scheme(path, rate_equations): Generate a LaTeX (dot2tex) diagram
-    of the reaction network.
+    - reaction_scheme(path, rate_equations): Generate a LaTeX (dot2tex) diagram
+        of the reaction network.
 
-- Random sampling:
+    - Random sampling:
 
-  - random_theta_phi(n): Uniform sampling on the unit sphere.
+    - random_theta_phi(n): Uniform sampling on the unit sphere.
 
-  - randomwalk_3d(...): Monte-Carlo 3D random walk with min/ max-distance
-    constraints (solution or spherical microreactor).
+    - randomwalk_3d(...): Monte-Carlo 3D random walk with min/ max-distance
+        constraints (solution or spherical microreactor).
 
 
-Key conventions
----------------
+Key conventions:
 
-- Rate maps use:  sink_state -> { source_state: Rate|float, ... }
-  The resulting sparse matrix M stores rates at (sink, source), so dP/dt = M P.
+    - Rate maps use:  sink_state -> { source_state: Rate|float, ... }
+    The resulting sparse matrix M stores rates at (sink, source), so dP/dt = M P.
 
-- Time evolution assumes uniform spacing in `time` and advances via the matrix
-  exponential of the rate matrix over Δt.
+    - Time evolution assumes uniform spacing in `time` and advances via the matrix
+    exponential of the rate matrix over Δt.
 
 
-Dependencies
-------------
+Dependencies:
 
-Relies on NumPy and SciPy (sparse) for numerics, graphviz + dot2tex for
-LaTeX diagram generation.
+    Relies on NumPy and SciPy (sparse) for numerics, graphviz + dot2tex for
+    LaTeX diagram generation.
 
 """
 

radicalpy/data.py

@@ -8,52 +8,52 @@
 and multiplicities. It also includes accessors for bundled reference data
 and convenience builders for frequently used composite objects.
 
-Contents
---------
-Conversions
-    - :func:`spin_to_multiplicity` – convert spin quantum number ``S`` to
-      multiplicity ``2S+1`` (validates half‐integer/integer input).
-    - :func:`multiplicity_to_spin` – invert the mapping, returning
-      ``S = (M-1)/2``.
-
-Data access
-    - :func:`get_data` – locate packaged data files (e.g. JSON databases).
-
-Core data structures
-    - :class:`Isotope` – look up an isotope’s spin multiplicity and
-      magnetogyric ratio from the bundled database (CODATA-style values),
-      with helpers like :py:meth:`Isotope.available`.
-    - :class:`Hfc` – hyperfine coupling container that supports either an
-      isotropic scalar value or a full 3×3 anisotropic tensor, and exposes
-      :py:meth:`Hfc.isotropic` and :py:meth:`Hfc.anisotropic`.
-    - :class:`Nucleus` – a nucleus within a molecule, defined by its
-      magnetogyric ratio, spin multiplicity, and hyperfine couplings; also
-      provides spin (Pauli) operator matrices via :py:meth:`Nucleus.pauli`.
-    - :class:`FuseNucleus` – an effective nucleus formed by fusing several
-      identical nuclei into a direct-sum representation to reduce Hilbert-space
-      dimension; includes utilities for validation and operator construction.
-    - :class:`Molecule` – collection of nuclei plus a radical (electron‐like
-      spin), with constructors for database-backed molecules
-      (:py:meth:`Molecule.fromdb`, :py:meth:`Molecule.all_nuclei`) and
-      ad-hoc assemblies (:py:meth:`Molecule.fromisotopes`).
-    - :class:`Triplet` – convenience molecule containing a single S=1
-      (multiplicity 3) radical and no nuclei.
-
-Units & conventions
--------------------
-- Magnetogyric ratios are commonly expressed as ``rad/s/T`` in the database;
-  many class properties expose ``gamma_mT`` (``rad/s/mT``) for convenience.
-- Hyperfine couplings (HFCs) are in mT. When a 3×3 tensor is provided,
-  the isotropic value is computed as ``trace(D)/3``.
-- Spin operators follow the usual convention with raising (``p``),
-  lowering (``m``), and Cartesian components (``x``, ``y``, ``z``).
-
-Error handling
---------------
-- :func:`spin_to_multiplicity` validates that ``S`` is integer or half-integer.
-- :class:`Isotope` raises ``ValueError`` for unknown symbols and provides
-  :py:meth:`Isotope.available` to enumerate valid options.
-- :class:`Hfc.anisotropic` raises ``ValueError`` when no tensor is available.
+Contents:
+
+    Conversions
+        - :func:`spin_to_multiplicity` – convert spin quantum number ``S`` to
+        multiplicity ``2S+1`` (validates half‐integer/integer input).
+        - :func:`multiplicity_to_spin` – invert the mapping, returning
+        ``S = (M-1)/2``.
+
+    Data access
+        - :func:`get_data` – locate packaged data files (e.g. JSON databases).
+
+    Core data structures
+        - :class:`Isotope` – look up an isotope’s spin multiplicity and
+        magnetogyric ratio from the bundled database (CODATA-style values),
+        with helpers like :py:meth:`Isotope.available`.
+        - :class:`Hfc` – hyperfine coupling container that supports either an
+        isotropic scalar value or a full 3×3 anisotropic tensor, and exposes
+        :py:meth:`Hfc.isotropic` and :py:meth:`Hfc.anisotropic`.
+        - :class:`Nucleus` – a nucleus within a molecule, defined by its
+        magnetogyric ratio, spin multiplicity, and hyperfine couplings; also
+        provides spin (Pauli) operator matrices via :py:meth:`Nucleus.pauli`.
+        - :class:`FuseNucleus` – an effective nucleus formed by fusing several
+        identical nuclei into a direct-sum representation to reduce Hilbert-space
+        dimension; includes utilities for validation and operator construction.
+        - :class:`Molecule` – collection of nuclei plus a radical (electron‐like
+        spin), with constructors for database-backed molecules
+        (:py:meth:`Molecule.fromdb`, :py:meth:`Molecule.all_nuclei`) and
+        ad-hoc assemblies (:py:meth:`Molecule.fromisotopes`).
+        - :class:`Triplet` – convenience molecule containing a single S=1
+        (multiplicity 3) radical and no nuclei.
+
+Units & conventions:
+
+    - Magnetogyric ratios are commonly expressed as ``rad/s/T`` in the database;
+    many class properties expose ``gamma_mT`` (``rad/s/mT``) for convenience.
+    - Hyperfine couplings (HFCs) are in mT. When a 3×3 tensor is provided,
+    the isotropic value is computed as ``trace(D)/3``.
+    - Spin operators follow the usual convention with raising (``p``),
+    lowering (``m``), and Cartesian components (``x``, ``y``, ``z``).
+
+Error handling:
+
+    - :func:`spin_to_multiplicity` validates that ``S`` is integer or half-integer.
+    - :class:`Isotope` raises ``ValueError`` for unknown symbols and provides
+    :py:meth:`Isotope.available` to enumerate valid options.
+    - :class:`Hfc.anisotropic` raises ``ValueError`` when no tensor is available.
 """
 
 from __future__ import annotations

radicalpy/estimations.py

@@ -575,7 +575,6 @@ def dipolar_interaction_anisotropic_from_dipolar_vector_without_prefactor(
 
     I3 = np.eye(3, dtype=float)
     rrT = np.outer(r_m, r_m)
-    # [(|r|^2) I - 3 r r^T] / |r|^5
     T = ((r2 * I3) - 3.0 * rrT) / (rnorm**5)
     return T
 
@@ -626,7 +625,6 @@ def dipolar_interaction_anisotropic_from_dipolar_vector(
         dipolar_vector, units=units
     )
     D = float(dipolar_prefactor) * T_geom
-    # return complex dtype to match other Hamiltonian builders smoothly
     return np.asarray(D, dtype=np.complex128)
 
 
@@ -919,7 +917,6 @@ def k_excitation(
             pathlength (float): The path length of the sample cell
                 (m).
 
-
     Returns:
             float: The excitation rate (1/s).
     """