radicalpy package

Submodules

radicalpy.classical module

Classical kinetics, diffusion, and utility helpers.

This module provides small, focused tools for building and simulating classical (first-order) kinetic schemes, along with diffusion step estimates, LaTeX formatting helpers, and basic 3D random-walk sampling.

Contents:

  • Diffusion:

  • get_delta_r(D, dt): RMS relative displacement √(6 D Δt).

  • LaTeX utilities:

  • latexify(rate_equations): Build d[Xi]/dt = … strings from a rate map.

  • latex_eqlist_to_align(eqs): Wrap equations in an align* environment.

  • Kinetics:

  • 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.

  • EquationRateResult: Convenience accessor to sum populations of one or more

    states over time.

  • Visualization:

  • reaction_scheme(path, rate_equations): Generate a LaTeX (dot2tex) diagram

    of the reaction network.

  • Random sampling:

  • 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).

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.

  • Time evolution assumes uniform spacing in time and advances via the matrix

exponential of the rate matrix over Δt.

Dependencies:

Relies on NumPy and SciPy (sparse) for numerics, graphviz + dot2tex for LaTeX diagram generation.

class radicalpy.classical.EquationRateResult(result, indices)[source]

Bases: object

Container for rate-equation time evolution results.

class radicalpy.classical.Rate(value, label)[source]

Bases: object

Scalar rate value paired with a LaTeX-formatted symbolic label.

Stores a numeric rate (or expression value) and a human-readable/LaTeX label that tracks algebraic manipulations via operator overloading.

label: str

LaTeX representation of the rate constant or expression.

value: float
class radicalpy.classical.RateEquations(rate_equations)[source]

Bases: object

Representation and solver for first-order kinetic rate equations.

property all_keys: list

List of all state labels in the network.

Returns:

Names of all states included in the equations.

Return type:

list[str]

are_valid_keys(keys)[source]

Check whether a set of state labels is valid.

Parameters:

keys (iterable[str]) – Collection of state labels to check.

Returns:

True if all labels are recognised states.

Return type:

bool

check_initial_states(initial_states)[source]

Validate an initial state population dictionary.

Ensures all keys are valid states and populations sum to 1.

Parameters:

initial_states (dict) – Mapping from state label to initial population (fractions).

Raises:

ValueError – If unknown state labels are included.

Notes

Prints a warning if the populations do not sum to unity.

time_evolution(time, initial_states)[source]

Evolve populations in time under the rate equations.

Solves dP/dt = M P using matrix exponentials at fixed time steps.

Parameters:

  • time (np.ndarray) – Array of time points (uniform spacing).

  • initial_states (dict) – Mapping of state → initial population. Must sum to 1.

Returns:

Object containing the full population time series and state index mapping.

Return type:

EquationRateResult

radicalpy.classical.get_delta_r(mutual_diffusion, delta_T)[source]

Root-mean-square step length for relative diffusion.

Uses ⟨Δr²⟩ = 6 D Δt for three-dimensional relative motion of two particles (mutual diffusion).

Parameters:

  • mutual_diffusion (float) – Mutual diffusion coefficient, D (m²/s).

  • delta_T (float) – Time interval, Δt (s).

Returns:

RMS displacement, √(6 D Δt) in meters (m).

Return type:

float

radicalpy.classical.latex_eqlist_to_align(eqs)[source]

Format LaTeX equations for alignment.

Wraps a list of LaTeX equations in an align* environment and inserts alignment markers before equal signs.

Parameters:

eqs (list[str]) – List of LaTeX equation strings.

Returns:

A single LaTeX string with all equations aligned.

Return type:

str

radicalpy.classical.latexify(rate_equations)[source]

Convert rate equations into LaTeX strings.

Produces a list of differential equations in LaTeX format, with terms written as sums over outgoing edges.

Parameters:

rate_equations (dict) – Dictionary mapping state labels (LHS) to dictionaries of transitions (RHS). Each RHS dict maps target states to Rate objects.

Returns:

List of LaTeX-formatted rate equations.

Return type:

list[str]

radicalpy.classical.random_theta_phi(num_samples=1)[source]

Randomly sample polar (θ) and azimuthal (φ) angles.

Samples directions uniformly over the unit sphere using inverse transform sampling.

Parameters:

num_samples (int, optional) – Number of random angle pairs to generate. Defaults to 1.

Returns:

Array of shape (2, num_samples) containing sampled angles in radians: - [0]: θ ∈ [0, π] (polar angle) - [1]: φ ∈ [0, 2π) (azimuthal angle)

Return type:

np.ndarray

radicalpy.classical.randomwalk_3d(n_steps, x_0, y_0, z_0, delta_r, r_min, r_max=0)[source]

Simulate a 3D Monte Carlo random walk of a radical pair.

Models diffusion of a radical in solution (r_max = 0) or inside a spherical microreactor (r_max > 0), enforcing reflective boundaries and minimum approach distance.

Parameters:

  • n_steps (int) – Number of simulation steps.

  • x_0 (float) – Initial x-coordinate (m).

  • y_0 (float) – Initial y-coordinate (m).

  • z_0 (float) – Initial z-coordinate (m).

  • delta_r (float) – Step length, e.g., from get_delta_r (m).

  • r_min (float) – Minimum allowed radical–radical distance (m).

  • r_max (float, optional) – Maximum radius of microreactor (m). Set to 0 for solution-based, and to a positive value for microreactor-based simulations.

Returns:

  • pos (ndarray): Radical positions at each step, shape (n_steps, 3).

  • dist (ndarray): Radical–radical distances at each step, shape (n_steps,).

  • angle (ndarray): Polar angles (θ) of displacements, shape (n_steps,).

Return type:

Tuple[np.ndarray, np.ndarray, np.ndarray]

radicalpy.classical.reaction_scheme(path, rate_equations)[source]

Generate a LaTeX reaction scheme diagram using Graphviz + dot2tex.

Builds a directed graph representation of the rate equations and writes the LaTeX TikZ code to a file.

Parameters:

  • path (str) – Output file path. Appends tex if not present.

  • rate_equations (dict) – Dictionary mapping states to outgoing transitions, with Rate objects providing LaTeX labels.

Returns:

Writes a tex file containing the diagram source.

Return type:

None

radicalpy.data module

Data models and utilities for spin-system simulations.

This module provides lightweight classes and helpers that describe the quantum‐mechanical ingredients used throughout the package (isotopes, nuclei, molecules) and common conversions between spin quantum numbers and multiplicities. It also includes accessors for bundled reference data and convenience builders for frequently used composite objects.

Contents:

Conversions

multiplicity 2S+1 (validates half‐integer/integer input). - multiplicity_to_spin() – invert the mapping, returning S = (M-1)/2.

Data access

  • get_data() – locate packaged data files (e.g. JSON databases).

Core data structures

  • Isotope – look up an isotope’s spin multiplicity and

magnetogyric ratio from the bundled database (CODATA-style values), with helpers like Isotope.available(). - Hfc – hyperfine coupling container that supports either an isotropic scalar value or a full 3×3 anisotropic tensor, and exposes Hfc.isotropic() and Hfc.anisotropic(). - Nucleus – a nucleus within a molecule, defined by its magnetogyric ratio, spin multiplicity, and hyperfine couplings; also provides spin (Pauli) operator matrices via Nucleus.pauli(). - 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. - Molecule – collection of nuclei plus a radical (electron‐like spin), with constructors for database-backed molecules (Molecule.fromdb(), Molecule.all_nuclei()) and ad-hoc assemblies (Molecule.fromisotopes()). - 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:

Isotope.available() to enumerate valid options. - Hfc.anisotropic raises ValueError when no tensor is available.

class radicalpy.data.FuseNucleus(magnetogyric_ratio, multiplicity, hfc, name=None, spinop=None, weight=None, direct_sum_info=None)[source]

Bases: Nucleus

Fuse identical nuclei into a single effective nucleus.

This class represents multiple identical nuclei that have been combined into a single effective nucleus for computational efficiency in spin dynamics calculations.

See also: Lindoy et al. Phys. Rev. Lett. 120, 220604 (2019).

Warning

This class should only be instantiated via the from_nuclei class method. Direct instantiation using __init__ is not recommended for end users. In addtion, for preparing the initial density matrix, the initial_density_matrix property combined with numpy.kron should be used instead of the pre-defined initial states such as radicalpy.simulation.State.SINGLET because the weights of each direct sum are not properly normalized.

Examples

Create a fused nucleus from three identical protons:

>>> protons = [Nucleus.fromisotope("1H", 1.5) for _ in range(3)]
>>> protons_fuse = FuseNucleus.from_nuclei(protons) # |J=3/2> ⊕ |J=1/2>
>>> assert protons_fuse.magnetogyric_ratio == protons[0].magnetogyric_ratio
>>> assert protons_fuse.multiplicity == 4 + 2
>>> assert protons_fuse.hfc.isotropic == protons[0].hfc.isotropic
>>> protons_fuse.name
'Fused1H(3)'

>>> protons_fuse.initial_density_matrix
array([[0.125, 0.   , 0.   , 0.   , 0.   , 0.   ],
       [0.   , 0.125, 0.   , 0.   , 0.   , 0.   ],
       [0.   , 0.   , 0.125, 0.   , 0.   , 0.   ],
       [0.   , 0.   , 0.   , 0.125, 0.   , 0.   ],
       [0.   , 0.   , 0.   , 0.   , 0.25 , 0.   ],
       [0.   , 0.   , 0.   , 0.   , 0.   , 0.25 ]])

>>> protons_fuse.pauli["z"]  # z-axis Pauli matrix
array([[ 1.5,  0. ,  0. ,  0. ,  0. ,  0. ],
       [ 0. ,  0.5,  0. ,  0. ,  0. ,  0. ],
       [ 0. ,  0. , -0.5,  0. ,  0. ,  0. ],
       [ 0. ,  0. ,  0. , -1.5,  0. ,  0. ],
       [ 0. ,  0. ,  0. ,  0. ,  0.5,  0. ],
       [ 0. ,  0. ,  0. ,  0. ,  0. , -0.5]])
classmethod from_nuclei(nuclei)[source]

Create a FuseNucleus from a list of identical nuclei.

This is the recommended way to create a FuseNucleus instance. The method validates that all nuclei are identical and computes the appropriate spin operators for the fused system.

Parameters:

nuclei (list[Nucleus]) – List of identical nuclei to fuse (minimum 2 nuclei)

Returns:

The fused nucleus with computed spin operators

Return type:

FuseNucleus

Raises:

ValueError – If fewer than 2 nuclei provided or nuclei are not identical

Examples

>>> h1 = Nucleus.fromisotope("1H", 1.5)
>>> h2 = Nucleus.fromisotope("1H", 1.5)
>>> h3 = Nucleus.fromisotope("1H", 1.5)
>>> h123 = FuseNucleus.from_nuclei([h1, h2, h3])
property initial_density_matrix: ndarray[tuple[int, ...], dtype[_ScalarType_co]]

Initial density matrix for the fused nucleus.

Returns:

Diagonal matrix with weights for initial state

Return type:

NDArray

static multiplicities_spin(N, I=Fraction(1, 2))[source]

Calculate the multiplicities for fusion of N identical spins.

Parameters:

  • N (int) – Number of spins to fuse

  • I (Fraction) – Spin quantum number of individual spins

Returns:

  • J is the total angular momentum quantum number

  • multiplicity is the number of times this J appears

Return type:

List of (J, multiplicity) tuples where

property pauli: dict

Return the spin operators dictionary.

class radicalpy.data.Hfc(hfc)[source]

Bases: object

The Hfc class represents isotropic and anisotropic HFC values.

Parameters:

hfc (float | list[list[float]]) – The HFC value. In case of a single float, only the isotropic value is set. In case of a 3x3 matrix both the isotropic and anisotropic values are stored.

Examples

Initialising the HFC with a 3-by-3 matrix (list of lists):

>>> with open(get_data("molecules/flavin_anion.json"), encoding="utf-8") as f:
...      flavin_dict = json.load(f)
>>> hfc_3d_data = flavin_dict["data"]["N5"]["hfc"]
>>> hfc_3d_obj = Hfc(hfc_3d_data)
>>> hfc_3d_obj
0.5141 <anisotropic available>

we can obtain both the isotropic value:

>>> hfc_3d_obj.isotropic
0.5141406139911681

and the anisotropic tensor:

>>> hfc_3d_obj.anisotropic
array([[-0.06819637,  0.01570029,  0.08701531],
       [ 0.01570029, -0.03652102,  0.27142597],
       [ 0.08701531,  0.27142597,  1.64713923]])

Initialising the HFC with a single float:

>>> with open(get_data("molecules/adenine_cation.json"), encoding="utf-8") as f:
...      adenine_dict = json.load(f)
>>> hfc_1d_data = adenine_dict["data"]["N6-H1"]["hfc"]
>>> hfc_1d_obj = Hfc(hfc_1d_data)
>>> hfc_1d_obj
-0.63 <anisotropic not available>

we can obtain both the isotropic value:

>>> hfc_1d_obj.isotropic
-0.63

but not the anisotropic tensor:

>>> hfc_1d_obj.anisotropic
Traceback (most recent call last):
...
ValueError: No anisotropic HFC data available.
property anisotropic: ndarray[tuple[int, ...], dtype[_ScalarType_co]]

Anisotropic value if available.

Returns:

The anisotropic HFC values.

Return type:

NDarray

property isotropic: float

Isotropic value.

Returns:

The isotropic HFC value.

Return type:

float

class radicalpy.data.Isotope(symbol)[source]

Bases: object

Class representing an isotope.

Parameters:

symbol (str) – The symbol/identifier of the isotope in the database.

Examples

Create an isotope using the database.

>>> E = Isotope("E")
>>> E
Symbol: E
Multiplicity: 2
Magnetogyric ratio: -176085963023.0
Details: {'name': 'Electron', 'source': 'CODATA 2018'}

Query the multiplicity:

>>> E.multiplicity
2

Query other details:

>>> E.details
{'name': 'Electron', 'source': 'CODATA 2018'}
classmethod available()[source]

List isotopes available in the database.

Returns:

List of available isotopes (symbols).

Return type:

list[str]

Examples

>>> available = Isotope.available()
>>> available[:5]
['G', 'E', 'N', 'M', 'P']

>>> Isotope(available[1])
Symbol: E
Multiplicity: 2
Magnetogyric ratio: -176085963023.0
Details: {'name': 'Electron', 'source': 'CODATA 2018'}

>>> Isotope(available[2])
Symbol: N
Multiplicity: 2
Magnetogyric ratio: -183247171.0
Details: {'name': 'Neutron', 'source': 'CODATA 2018'}
property gamma_mT

Return gamma value in rad/s/mT.

property spin_quantum_number: float

Spin quantum numer of Isotope.

class radicalpy.data.Molecule(name='', nuclei=[], radical=E(-176085963023.0, 2, 0.0 <anisotropic not available>), info={})[source]

Bases: object

Representation of a molecule for the simulation.

A molecule is described by a name and a list of nuclei (see Nucleus). Using the default constructor is cumbersome and not recommended. The preferred way is to use the following convenience methods (click on the method name to see its documentation):

Parameters:

  • name (str) – The name of the Molecule.

  • nuclei (list[Nucleus]) – List of nuclei/atoms which should be simulated (see Nucleus).

  • radical (Nucleus) – The radical of the molecule. (Default Nucleus.fromisotope("E", 0.0)).

  • info (dict[str, str]) – Dictionary of miscellaneous information about the molecule.

Examples

The default constructor takes an arbitrary name and a list of molecules to construct a molecule.

>>> gamma_1H = 267522.187
>>> gamma_14N = 19337.792
>>> Molecule("kryptonite", [Nucleus(gamma_1H, 2, Hfc(1.0), "Hydrogen"),
...                            Nucleus(gamma_14N, 3, Hfc(-0.5), "Nitrogen")])
Molecule: kryptonite
Nuclei:
  Hydrogen(267522186.99999997, 2, 1.0 <anisotropic not available>)
  Nitrogen(19337792.0, 3, -0.5 <anisotropic not available>)
Radical: E(-176085963023.0, 2, 0.0 <anisotropic not available>)

Or alternatively:

>>> gammas = [267522.187, 19337.792]
>>> multis = [2, 3]
>>> hfcs = [1.0, -0.5]
>>> names = ["Hydrogen", "Nitrogen"]
>>> params = zip(gammas, multis, map(Hfc, hfcs), names)
>>> Molecule("kryptonite", [Nucleus(*param) for param in params])
Molecule: kryptonite
Nuclei:
  Hydrogen(267522186.99999997, 2, 1.0 <anisotropic not available>)
  Nitrogen(19337792.0, 3, -0.5 <anisotropic not available>)
Radical: E(-176085963023.0, 2, 0.0 <anisotropic not available>)
classmethod all_nuclei(name)[source]

Construct a molecule from the database with all nuclei.

Parameters:

name (str) – A name of the molecule available in the database (see Molecule.available()).

Examples

>>> Molecule.all_nuclei("flavin_anion")
Molecule: flavin_anion
Nuclei:
  14N(19337792.0, 3, 0.5141 <anisotropic available>)
  14N(19337792.0, 3, -0.001275 <anisotropic available>)
  14N(19337792.0, 3, -0.03654 <anisotropic available>)
  1H(267522187.44, 2, 0.05075 <anisotropic available>)
  1H(267522187.44, 2, -0.1371 <anisotropic available>)
  1H(267522187.44, 2, -0.1371 <anisotropic available>)
  1H(267522187.44, 2, -0.1371 <anisotropic available>)
  1H(267522187.44, 2, -0.4403 <anisotropic available>)
  1H(267522187.44, 2, 0.4546 <anisotropic available>)
  1H(267522187.44, 2, 0.4546 <anisotropic available>)
  1H(267522187.44, 2, 0.4546 <anisotropic available>)
  1H(267522187.44, 2, 0.009597 <anisotropic available>)
  1H(267522187.44, 2, 0.4263 <anisotropic available>)
  1H(267522187.44, 2, 0.4233 <anisotropic available>)
  1H(267522187.44, 2, -0.02004 <anisotropic available>)
  14N(19337792.0, 3, 0.1784 <anisotropic available>)
Radical: E(-176085963023.0, 2, 0.0 <anisotropic not available>)
Info: {'units': 'mT', 'name': 'Flavin radical anion'}
classmethod available()[source]

List molecules available in the database.

Returns:

List of available molecules (names).

Return type:

list[str]

Examples

>>> available = Molecule.available()
>>> available[:4]
['2_6_aqds', 'adenine_cation', 'fad', 'flavin_anion']
custom: bool
property effective_hyperfine: float

Effective hyperfine for the entire molecule.

classmethod fromdb(name, nuclei=[])[source]

Construct a molecule from the database.

Parameters:

  • name (str) – A name of the molecule available in the database (see Molecule.available()).

  • nuclei (list[str]) – A list of nuclei names found in the molecule specified by name.

Examples

>>> Molecule.fromdb("flavin_anion", nuclei=["N14"])
Molecule: flavin_anion
Nuclei:
  14N(19337792.0, 3, -0.001275 <anisotropic available>)
Radical: E(-176085963023.0, 2, 0.0 <anisotropic not available>)
Info: {'units': 'mT', 'name': 'Flavin radical anion'}
classmethod fromisotopes(isotopes, hfcs, name='')[source]

Construct molecule from isotopes.

Parameters:

  • isotopes (list[str]) – A list of Isotope names found in the isotope database.

  • hfcs (list) – A list of HFC values (see Hfc constructor).

  • name (str) – An optional name for the molecule.

Examples

>>> Molecule.fromisotopes(isotopes=["1H", "14N"],
...                          hfcs=[1.5, 0.9],
...                          name="kryptonite")
Molecule: kryptonite
Nuclei:
  1H(267522187.44, 2, 1.5 <anisotropic not available>)
  14N(19337792.0, 3, 0.9 <anisotropic not available>)
Radical: E(-176085963023.0, 2, 0.0 <anisotropic not available>)
info: dict[str, str]
classmethod load_molecule_json(molecule)[source]

Load a molecule definition from the bundled JSON database.

Parameters:

molecule (str) – Molecule name (filename without json).

Returns:

Parsed JSON content with info and data sections.

Return type:

dict

Raises:
name: str
nuclei: list[Nucleus]
radical: Nucleus
property semiclassical_std: float

Standard deviation σ of the Schulten–Wolynes (SW) effective hyperfine field.

Definitions (SW isotropic case):

τ^{-2} = (1/6) Σ_k a_k^2 I_k (I_k + 1) σ = √2 / τ ⇒ σ^2 = (1/3) Σ_k a_k^2 I_k (I_k + 1)

Args:

units: {“omega”, “mT”}, default “omega”

  • “omega”: return σ_ω in angular-frequency units (same units as a_k).

  • “mT”return σ_B in millitesla, i.e., σ_B = σ_ω / γ, where

    γ = self.radical.gamma_mT (rad s^-1 mT^-1).

ak_units : {“angular”, “Hz”}, default “angular”

  • “angular”: each n.hfc.isotropic is already in angular frequency units (e.g., rad s^-1).

  • “Hz” : each n.hfc.isotropic is in cycles per second; converted internally via 2π.

Returns:

σ in the requested units.

Notes:

  • This computes the isotropic SW width (scalar). For anisotropic A-tensors, the isotropic equivalent variance is (1/3) Tr(Σ_ω), where

    Σ_ω = (1/3) Σ_k I_k(I_k+1) A_k A_k^T.

    If you need full anisotropic sampling, construct Σ_ω explicitly elsewhere.

class radicalpy.data.Nucleus(magnetogyric_ratio, multiplicity, hfc, name=None)[source]

Bases: object

A nucleus in a molecule.

Construct a nucleus from an Isotope and an Hfc.

>>> Nucleus.fromisotope("1H", 1.1)
1H(267522187.44, 2, 1.1 <anisotropic not available>)

The default constructor needs the magnetogyric ratio, the multiplicity and the HFC values.

>>> Nucleus(0.001, 2, Hfc(3.0))
Nucleus(1.0, 2, 3.0 <anisotropic not available>)

Additionally a name can also be added.

>>> Nucleus(0.001, 2, Hfc(3.0), "Adamantium")
Adamantium(1.0, 2, 3.0 <anisotropic not available>)
classmethod fromisotope(isotope, hfc)[source]

Construct a Nucleus from an Isotope.

Parameters:
Returns:

A nucleus with magnetogyric ratio, multiplicity

and name determined by the isotope and the hfc value.

Return type:

Nucleus

property gamma_mT

Return magnetogyric ratio, \(\gamma\) (rad/s/mT).

hfc: Hfc
magnetogyric_ratio: float
multiplicity: int
name: Optional[str]
property pauli

Generate Pauli matrices.

Generates the Pauli matrices corresponding to a given multiplicity.

Parameters:

mult (int) – The multiplicity of the element.

Returns:

A dictionary containing 6 np.array matrices of shape (mult, mult):

  • the unit operator result["u"],

  • raising operator result["p"],

  • lowering operator result["m"],

  • Pauli matrix for x axis result["x"],

  • Pauli matrix for y axis result["y"],

  • Pauli matrix for z axis result["z"].

Return type:

dict

property spin_quantum_number: float

Spin quantum numer of Isotope.

class radicalpy.data.Triplet[source]

Bases: Molecule

radicalpy.data.get_data(suffix='')[source]

Get the directory containing data files.

Return type:

Traversable

radicalpy.data.multiplicity_to_spin(multiplicity)[source]

Spin multiplicity to spin quantum number.

Parameters:

multiplicity (int) – Spin multiplicity.

Returns:

Spin quantum number.

Return type:

float

radicalpy.data.spin_to_multiplicity(spin)[source]

Spin quantum number to multiplicity.

Parameters:

spin (float) – Spin quantum number.

Returns:

Spin multiplicity.

Return type:

int

radicalpy.estimations module

Estimation utilities for fields, rates, interactions, and transport.

This module collects closed-form and semi-empirical estimators commonly used in spin chemistry and magnetic resonance modeling. Functions cover theoretical B₁/₂ values, T₁/T₂ relaxation rates (g-anisotropy, hyperfine tensor, and tumbling motion), diffusion/viscosity, dipolar and exchange interactions, kinetic rates (excitation, recombination, re-encounter, electron transfer), triplet relaxation, rotational correlation times, and helpers for fitting autocorrelation functions.

Major groups:
B½ / field scales
Relaxation (g-anisotropy / tumbling)
Transport & correlation
Spin interactions
Dephasing / mixing from trajectories & geometry
Kinetics
Units & conventions:

  • Fields: millitesla (mT) unless otherwise noted; some functions expect tesla (T).

  • Rates: s⁻¹; angular frequencies in rad·s⁻¹ where specified.

  • Distances: metres (m); protein electron-transfer separation is in Å as per literature fits.

  • Temperatures: kelvin (K) unless otherwise noted (e.g., viscosity uses °C input).

  • Time arrays (ts, time): seconds (s).

  • g-tensor inputs: principal components per radical, dimensionless.

  • Many expressions reuse physical constants from a shared C namespace and isotope data via Isotope("E") (electron); see your codebase for those definitions.

Notes

  • B½ estimators (Weller, Golesworthy) are empirical/approximate and assume model conditions stated in the cited works.

  • Tumbling-motion formulas (Bloembergen–Purcell–Pound) require a distance r and correlation time tau_c; results scale as r⁻⁶.

  • autocorrelation_fit uses a fixed set of τ values (log-spaced) with a non-negative multi-exponential fit; the effective tau_c is the amplitude-weighted sum.

  • Viscosity model (Volk) is accurate within stated temperature ranges; ensure frac_glyc ∈ [0, 1].

  • Some helpers accept scalars or NumPy arrays and broadcast accordingly.

References

Requirements:

  • numpy, scipy (for curve fitting in autocorrelation_fit), and project constants/utilities (C, Isotope, utils) available in scope.

See also

  • relaxation.py, kinetics.py for superoperators using several of these rates.

  • experiments.py for simulations (MARY, EPR/ODMR/OMFE) that can consume these estimates.

radicalpy.estimations.Bhalf_theoretical_hyperfine(sim)[source]

Theoretical B1/2 for radical pairs. Estimated with hyperfine interactions.

Source: Weller et al. Chem. Phys. Lett. 96, 1, 24-27 (1983).

Parameters:

sim (HilbertSimulation) – The sim object containing the hyperfine coupling constants. It should contain exactly two molecules.

Returns:

The B1/2 value (mT).

Return type:

float

Todo

Change sim to a list of molecules.

radicalpy.estimations.Bhalf_theoretical_relaxation(kstd, krec)[source]

Theoretical B1/2 for radical pairs. Estimated with spin dephasing rate.

Source: Golesworthy et al. J. Chem. Phys. 159, 105102 (2023).

Parameters:

  • kstd (float) – Singlet-triplet dephasing rate (1/s).

  • krec (float) – Recombination rate (1/s).

Returns:

The B1/2 value (mT).

Return type:

float

radicalpy.estimations.Bhalf_theoretical_relaxation_delay(kstd, krec, td)[source]

Theoretical B1/2 for radical pairs. Estimated with spin dephasing rate and pump-probe delay time.

Source: Golesworthy et al. J. Chem. Phys. 159, 105102 (2023).

Parameters:

  • kstd (float) – Singlet-triplet dephasing rate (1/s).

  • krec (float) – Recombination rate (1/s).

  • td (float or np.ndarray) – Pump-probe delay (s).

Returns:

The B1/2 value (mT).

Return type:

float

radicalpy.estimations.T1_relaxation_rate_g_tensor(g_tensors, B, tau_c)[source]

T1 relaxation rate.

Estimate T1 relaxation rate based on tau_c and g-tensor anisotropy.

Source: Hayashi, Introduction to Dynamic Spin Chemistry: Magnetic Field Effects on Chemical and Biochemical Reactions (2004).

Parameters:

  • g_tensors (list) – The principle components of g-tensor.

  • B (float or np.ndarray) – The external magnetic field strength (T).

  • tau_c (float or np.ndarray) – The rotational correlation time (s).

Returns:

The T1 relaxation rate (1/s)

Return type:

float or np.ndarray

radicalpy.estimations.T1_relaxation_rate_hyperfine_tensor(hyperfine_tensor, B, tau_c)[source]

T1 relaxation rate.

Estimate T1 relaxation rate based on tau_c and hyperfine tensor anisotropy.

Source: Carrington and McLachlan, Introduction to Magnetic Resonance with Applications to Chemistry and Chemical Physics (1967).

Parameters:

  • hyperfine_tensor (np.ndarray) – Hyperfine tensor (mT).

  • B (float) – The external magnetic field strength (mT).

  • tau_c (float) – The rotational correlation time (s).

Returns:

The T1 relaxation rate (1/s)

Return type:

float

radicalpy.estimations.T1_relaxation_rate_tumbling_motion(tau_c, B0, r)[source]

T1 relaxation rate.

Estimate T1 relaxation rate based on tau_c and r distance between radicals.

Source: Bloembergen, Purcell, and Pound, Phys. Rev. 73, 679 (1948).

Parameters:

  • tau_c (float or np.ndarray) – The rotational correlation time (s).

  • B0 (float or np.ndarray) – The external magnetic field strength (T).

  • r (float or np.ndarray) – The distance between radicals (m).

Returns:

The T1 relaxation rate (1/s).

Return type:

float or np.ndarray

radicalpy.estimations.T2_relaxation_rate_g_tensor(g_tensors, B, tau_c)[source]

T2 relaxation rate.

Estimate T2 relaxation rate based on tau_c and g-tensor anisotropy.

Source: Hayashi, Introduction to Dynamic Spin Chemistry: Magnetic Field Effects on Chemical and Biochemical Reactions (2004).

Parameters:

  • g_tensors (list) – The principle components of g-tensor.

  • B (float or np.ndarray) – The external magnetic field strength (T).

  • tau_c (float or np.ndarray) – The rotational correlation time (s).

Returns:

The T2 relaxation rate (1/s).

Return type:

float or np.ndarray

radicalpy.estimations.T2_relaxation_rate_hyperfine_tensor(hyperfine_tensor, B, tau_c)[source]

T2 relaxation rate.

Estimate T2 relaxation rate based on tau_c and hyperfine tensor anisotropy.

Source: Carrington and McLachlan, Introduction to Magnetic Resonance with Applications to Chemistry and Chemical Physics (1967).

Parameters:

  • hyperfine_tensor (np.ndarray) – Hyperfine tensor (mT).

  • B (float) – The external magnetic field strength (mT).

  • tau_c (float) – The rotational correlation time (s).

Returns:

The T2 relaxation rate (1/s)

Return type:

float

radicalpy.estimations.T2_relaxation_rate_tumbling_motion(tau_c, B0, r)[source]

T2 relaxation rate.

Estimate T2 relaxation rate based on tau_c and r distance between radicals.

Source: Bloembergen, Purcell, and Pound, Phys. Rev. 73, 679 (1948).

Parameters:

  • tau_c (float or np.ndarray) – The rotational correlation time (s).

  • B0 (float or np.ndarray) – The external magnetic field strength (T).

  • r (float or np.ndarray) – The distance between radicals (m).

Returns:

The T2 relaxation rate (1/s).

Return type:

float or np.ndarray

radicalpy.estimations.aqueous_glycerol_viscosity(frac_glyc, temp)[source]

Viscosity of aqueous glycerol solutions.

Gives a good approximation for temperatures in the range 0-100°C.

Source: Volk et al. Experiments in Fluids, 59, 76, (2018).

Parameters:

  • frac_glyc (float or np.ndarray) – The fraction of glycerol in solution (0.00-1.00).

  • temp (float) – The temperature in °C (0-100) (<0.07% accuracy between 15-30°C).

Returns:

The viscosity of the glycerol/water mixture in N s/m^2.

Return type:

float or np.ndarray

radicalpy.estimations.autocorrelation_fit(ts, trajectory, tau_begin, tau_end, num_exp=100, normalise=False)[source]

Fit multiexponential to autocorrelation plot and calculate the effective rotational correlation time.

Parameters:

  • ts (np.ndarray) – Time interval (x-axis of the trajectory) (s).

  • trajectory (np.ndarray) – The raw data which will be fit and used to calculate tau_c.

  • tau_begin (float) – Initial lag time (s).

  • tau_end (float) – Final lag time (s).

  • num_exp (int) – Number of exponential terms in the multiexponential fit (default=100).

  • normalise (bool) – When set to true, the autocorrelation is normalised (default=False).

Returns:

  • fit is the multiexponential fit to the autocorrelation.

  • tau_c is the effective rotational correlation time.

Return type:

dict

Thank you, Gesa Grüning!

radicalpy.estimations.diffusion_coefficient(radius, temperature, eta)[source]

Diffusion coefficient.

The Stokes-Einstein relation.

Source: Einstein, Ann. der Physik, 17, 549-560 (1905).

Parameters:

  • radius (float) – The radius of the molecule (m).

  • temperature (float) – The temperature of the solution (K).

  • eta (float) – The viscosity of the solution (kg/m/s).

Returns:

The diffusion coefficient (m^2/s).

Return type:

float

radicalpy.estimations.dipolar_interaction_MC(r, theta)[source]

Dipolar interaction for Monte Carlo trajectories.

Sources:

Parameters:

  • r (float | np.ndarray) – The interradical separation (m).

  • theta (float | np.ndarray) – The angle of molecular rotation (rad).

Returns:

The dipolar coupling constant in milli Tesla (mT).

Return type:

float | np.ndarray

radicalpy.estimations.dipolar_interaction_anisotropic(r)[source]

Anisotropic dipolar coupling.

Point dipole approximation is used.

Parameters:

r (float | ndarray) – The interradical separation (m).

Return type:

ndarray

Returns:

The dipolar coupling tensor in millitesla (mT).

Todo

np.ndarray not implemented. dipolar * diag fails.

radicalpy.estimations.dipolar_interaction_anisotropic_from_dipolar_vector(dipolar_vector, dipolar_prefactor=3.2698330438079847e-19, *, units='Å')[source]

Full electron–electron dipolar coupling tensor (angular frequency, rad/s).

Thank you, Luca Gerhards!

Computes

\[D = prefactor * [ (|r|^2 I) - 3 r rᵀ ] / |r|^5\]

where the default prefactor is

\[(μ0 / 4π) * ħ * γ_e^2\]

with γ_e the electron gyromagnetic ratio (rad s⁻¹ T⁻¹). The input displacement vector can be given in Å or m.

Parameters:

  • dipolar_vector (ndarray) – Displacement vector r from spin 2 to spin 1 (Cartesian).

  • dipolar_prefactor (float) – Physical prefactor. Default corresponds to two electrons and returns the tensor in angular frequency units (rad/s).

  • units (str) – Units of dipolar_vector. {“Å”,”m”}, optional. Default “Å” (Angstrom).

Return type:

ndarray

Returns:

Dipolar coupling tensor in angular frequency units (rad/s).

radicalpy.estimations.dipolar_interaction_anisotropic_from_dipolar_vector_without_prefactor(dipolar_vector, *, units='Å')[source]

Geometry-only dipolar tensor kernel (no physical prefactor).

Thank you, Luca Gerhards!

Given a displacement vector r between two spins, returns the 3×3 geometric tensor

\[T_{geom} = [ (|r|^2 I) - 3 r rᵀ ] / |r|^5 = (I - 3 r̂ r̂ᵀ) / |r|^3\]
Parameters:

  • dipolar_vector (ndarray) – Displacement vector r from spin 2 to spin 1 (Cartesian).

  • units (str) – Units of dipolar_vector “Å” or “m”. Default “Å” (Angstrom).

Return type:

ndarray

Returns:

Geometry-only tensor in SI length units (i.e. 1/m^3 factor embedded).

radicalpy.estimations.dipolar_interaction_isotropic(r)[source]

Isotropic dipolar coupling.

Point dipole approximation is used.

Source: Santabarbara et al. Biochemistry, 44, 6, 2119–2128 (2005).

Parameters:

r (float or np.ndarray) – The interradical separation (m).

Returns:

The dipolar coupling constant in millitesla (mT), which is negative.

Return type:

(float or np.ndarray)

radicalpy.estimations.dipolar_interaction_point_dipole(r12)[source]

Dipolar interaction for point dipole approximation.

Parameters:

r12 (tuple[float, float, float]) – The interradical separation (m).

Returns:

The 3x3 dipolar coupling tensor in millitesla (mT)

Return type:

np.ndarray

Example

>>> dipolar_interaction_point_dipole(r12=(1e-09, 0, 0))
array([[-3.71390565, -0.        , -0.        ],
        [-0.        ,  1.85695282, -0.        ],
        [-0.        , -0.        ,  1.85695282]])
radicalpy.estimations.exchange_interaction_in_protein(r, beta=14000000000.0, J0=9700000000.0)[source]

Exchange interaction for radical pairs in proteins.

Source: Moser et al. Nature 355, 796–802 (1992).

Parameters:

  • r (float or np.ndarray) – The interradical separation (m).

  • beta (float) – The range parameter (1/m).

  • J0 (float) – The strength of the interaction (mT).

Returns:

The exchange interaction (mT).

Return type:

(float or np.ndarray)

radicalpy.estimations.exchange_interaction_in_solution(r, beta=4.9e-11, J0rad=1.7e+17)[source]

Exchange interaction for radical pairs in solution.

Source: McLauchlan et al. Mol. Phys. 73:2, 241-263 (1991).

Parameters:

  • r (float or np.ndarray) – The interradical separation (m).

  • beta (float) – The range parameter (m).

  • J0rad (float) – The strength of the interaction (rad/s).

Returns:

The exchange interaction (mT).

Return type:

(float or np.ndarray)

radicalpy.estimations.exchange_interaction_in_solution_MC(r, beta=20000000000.0, J0=-570)[source]

Exchange interaction for Monte Carlo trajectories.

Sources:

Parameters:

  • r (np.ndarray) – The interradical separation (m).

  • beta (float) – The range parameter (1/m).

  • J0 (float) – The strength of the interaction (mT).

Returns:

The exchange coupling constant in milli Tesla (mT).

Return type:

np.ndarray

radicalpy.estimations.g_tensor_relaxation_rate(tau_c, g1, g2)[source]

g-tensor relaxation rate.

To be used with radicalpy.relaxation.RandomFields.

Source: Player et al. J. Chem. Phys. 153, 084303 (2020).

Parameters:

  • tau_c (float) – The rotational correlation time (s).

  • g1 (list) – The principle components of g-tensor of the first radical.

  • g2 (list) – The principle components of g-tensor of the second radical.

Returns:

The g-tensor relaxation rate (1/s).

Return type:

float

radicalpy.estimations.k_D(D, tau_c)[source]

D (dipolar)-dephasing rate for trajectories.

Source: Kattnig et al. New J. Phys., 18, 063007 (2016).

Parameters:

  • D (np.ndarray) – The dipolar interaction trajectory (mT).

  • tau_c (float) – The rotational correlation time (s).

Returns:

The D-dephasing rate (1/s).

Return type:

float

radicalpy.estimations.k_STD(J, tau_c)[source]

ST-dephasing rate for trajectories.

Source: Kattnig et al. New J. Phys., 18, 063007 (2016).

Parameters:

  • J (np.ndarray) – The exchange interaction trajectory (mT).

  • tau_c (float) – The rotational correlation time (s).

Returns:

The ST-dephasing rate (1/s).

Return type:

float

radicalpy.estimations.k_STD_microreactor(D, V, d=5e-10, J0=100000000000.0, alpha=20000000000.0)[source]

ST-dephasing rate for radical pairs in microreactors.

Source: Shushin, Chem. Phys. Lett., 181, 2–3, 274-278 (1991).

Parameters:

  • D (float) – The mutual diffusion coefficient (m^2/s).

  • V (float) – The volume of the microreactor (e.g. micelle) (m^3).

  • d (float) – The distance of closest approach (m).

  • J0 (float) – The maximum exchange interaction (1/s).

  • alpha (float) – The characteristic length factor (1/m).

Returns:

The ST-dephasing rate (1/s).

Return type:

float

radicalpy.estimations.k_ST_mixing(Bhalf)[source]

Singlet-triplet mixing rate.

Source: Steiner et al. Chem. Rev. 89, 1, 51–147 (1989).

Parameters:

Bhalf (float) – The theoretical B1/2 value (mT).

Returns:

The ST-mixing rate (1/s).

Return type:

float

radicalpy.estimations.k_constant(r, gamma)[source]

K constant used for T1 and T2 estimation.

K constant used to calculate T1 and T2 relaxation radicalpy.estimations.k_t1_relaxation_tumbling_motion radicalpy.estimations.k_t2_relaxation_tumbling_motion.

Source: Bloembergen, Purcell, and Pound, Phys. Rev. 73, 679 (1948).

Parameters:

  • tau_c (float or np.ndarray) – The rotational correlation time (s).

  • gamma (float) – The magnetogyric ratio (rad/s/T).

Returns:

K constant (1/s).

Return type:

float or np.ndarray

radicalpy.estimations.k_electron_transfer(separation, driving_force=-1, reorganisation_energy=1)[source]

Electron transfer rate.

The default values (when -driving_force == reorganisation_energy) return the maximum electron transfer rate.

Source: Moser et al. Biochim. Biophys. Acta Bioenerg. 1797, 1573‐1586 (2010).

Parameters:

  • separation (float) – The edge-to-edge separation, R (Å).

  • driving_force (float) – The driving force, \(\Delta G\) (eV).

  • reorganisation_energy (float) – The reorganisation energy, \(\lambda\) (eV).

Returns:

The electron transfer rate (1/s).

Return type:

float

radicalpy.estimations.k_excitation(wavelength, beam_radius, absorbance, concentration, laser_power, path_length)[source]

Groundstate excitation rate.

Parameters:

  • wavelength (float) – The excitation wavelength (m).

  • beam_radius (float) – Radius of the beam spot (m).

  • absorbance (float) – Absorbance of the sample (OD).

  • concentration (float) – Concentration of the sample (mol/m^3).

  • laser_power (float) – The excitation laser power (W).

  • pathlength (float) – The path length of the sample cell (m).

Returns:

The excitation rate (1/s).

Return type:

float

radicalpy.estimations.k_recombination(MFE, k_escape)[source]

Singlet recombination rate.

Source: Maeda et al. Mol. Phys., 117:19, 2709-2718 (2019).

Parameters:

  • MFE (float) – The magnetic field effect (0.00-1.00).

  • k_escape (float) – The free radical formation rate constant (1/s).

Returns:

The singlet recombination rate (1/s).

Return type:

float

radicalpy.estimations.k_reencounter(encounter_dist, diff_coeff)[source]

Re-encounter rate.

Source: Salikhov, J. Magn. Reson., 63, 271-279 (1985).

Parameters:

  • encounter_dist (float) – The effective re-encounter distance, R* (m).

  • diff_coeff (float) – The diffusion coefficient, D (m^2/s).

Returns:

The re-encounter rate (1/s).

Return type:

float

radicalpy.estimations.k_triplet_relaxation(B0, tau_c, D, E)[source]

Excited triplet state relaxation rate.

Source: Atkins et al. Mol. Phys., 27, 6 (1974).

Parameters:

  • B0 (float) – The external magnetic field (MHz).

  • tau_c (float) – The rotational correlation time (s).

  • D (float) – The zero field splitting (ZFS) parameter D (Hz).

  • E (float) – The zero field splitting (ZFS) parameter E (Hz).

Returns:

The excited triplet state relaxation rate (1/s).

Return type:

float

radicalpy.estimations.rotational_correlation_time_for_molecule(radius, temp, eta=0.00089)[source]

Rotational correlation time.

Rotational correlation time is the average time it takes for a molecule (smaller than a protein) to rotate one radian. For proteins see radicalpy.estimations.rotational_correlation_time_for_protein.

To calculate viscosity (eta) for glycerol-water mixtures see radicalpy.estimations.aqueous_glycerol_viscosity.

Parameters:

  • radius (float) – The radius of a spherical molecule (m).

  • temp (float) – The temperature of the solution (K).

  • eta (float) – The viscosity of the solution (N s/m^2) (default: 0.89e-3 corresponds to water).

Returns:

The rotational correlation time (s).

Return type:

float

radicalpy.estimations.rotational_correlation_time_for_protein(Mr, temp, eta=0.00089)[source]

Rotational correlation time.

Rotational correlation time is the average time it takes for a protein to rotate one radian. For small molecules see radicalpy.estimations.rotational_correlation_time_for_molecule.

To calculate viscosity (eta) for glycerol-water mixtures see radicalpy.estimations.aqueous_glycerol_viscosity.

Source: Cavanagh et al. Protein NMR Spectroscopy. Principles and Practice, Elsevier Academic Press (2007).

Parameters:

  • Mr (float) – The molecular weight of the protein (kDa).

  • temp (float) – The temperature of the solution (K).

  • eta (float) – The viscosity of the solution (N s/m^2) (default: 0.89e-3 corresponds to water).

Returns:

The rotational correlation time (s).

Return type:

float

radicalpy.experiments module

Experimental simulation routines for spin chemistry and magnetic resonance.

This module groups high-level experiment drivers and inner loops used to simulate anisotropy scans, magnetically affected reaction yields (MARY), EPR/ODMR/OMFE time-domain responses, RYDMR, steady-state signals, simple NMR lineshapes, OOP-ESEEM envelopes, and hybrid semiclassical/quantum MARY variants. Most functions assemble Hamiltonian/Liouvillian terms, apply kinetics/relaxation superoperators, propagate density matrices, and return structured results (dicts, arrays) ready for plotting and analysis.

Functions:

  • anisotropy: Full anisotropy experiment wrapper; returns time evolutions and yields.

  • anisotropy_loop: Inner loop over (θ, φ) orientations; propagates and returns product probabilities.

  • cidnp: CIDNP polarisation vs field for S–T0 mixing.

  • coherent_control: Coherent control of radical pair spin dynamics.

  • epr: CW/AC time-domain EPR vs B0 with B1 drive and frequency offset.

  • field_switching: Nanosecond field switching experiment (SEMF).

  • kine_quantum_mary: Hybrid kinetic+quantum MARY with stochastic hyperfine sampling.

  • magnetic_field_loop: Inner loop over swept field B; returns time-resolved density matrices.

  • magnetic_field_loop_semiclassical: Semiclassical inner loop over swept field B; returns time-resolved density matrices.

  • mary_lfe_hfe: Post-process product probabilities → MARY, low-/high-field effects (%).

  • mary: MARY vs B (time-domain propagation + yields + normalised response).

  • mary_semiclassical: Semiclassical MARY vs B (time-domain propagation + yields + normalised response).

  • modulated_mary_brute_force: Lock-in MARY via phase randomisation and numerical integration.

  • nmr: Simple 1D NMR synthesiser (FID → FFT) with multiplets and T2 decay.

  • odmr: ODMR vs RF frequency (B1_freq) at fixed B0.

  • omfe: Oscillating magnetic field effect vs RF frequency in transverse field.

  • oop_eseem: Out-of-phase ESEEM envelope via Gauss–Legendre quadrature.

  • rydmr: Reaction yield–detected magnetic resonance vs static field.

  • semiclassical_mary: Semiclassical MARY with explicit population channels.

  • steady_state_mary: Steady state via linear solve of Liouvillian (with ZFS/exchange/Zeeman).

Usage pattern:

  1. Build Hamiltonian components from your Simulation object (e.g., Zeeman, exchange J, dipolar D, hyperfine ± anisotropy).

  2. (Optional) Assemble kinetics/relaxation terms and apply to the Liouvillian via apply_liouville_hamiltonian_modifiers(...).

  3. Propagate using an inner loop (anisotropy_loop, magnetic_field_loop) or a high-level routine (mary, epr, odmr, omfe, rydmr).

  4. Post-process to yields, MARY/LFE/HFE, or spectra; the helpers return a dict containing inputs, intermediates, and outputs.

Args conventions (selected):

  • time (np.ndarray): Uniform time grid (s).

  • B, B0, B1, B1_freq (array/float): Fields in mT unless noted.

  • D, E, J (float): Dipolar/ZFS/Exchange parameters (mT) for Zeeman-scaled forms, or (rad/s) where indicated (e.g., oop_eseem).

  • theta, phi (float or np.ndarray): Polar/azimuthal angles (rad) for anisotropy.

  • kinetics, relaxations (list): Collections of (super)operators compatible with the simulation API; applied in Liouville space.

  • multiplets (NMR): [n_nuc, chem_shift_Hz, multiplicity, J_Hz] per line.

Returns (typical):

  • Dictionaries with keys such as: time, B/B0/B1_freq, angles, rhos, time_evolutions (product probabilities), product_yields, product_yield_sums, and derived metrics (MARY, LFE, HFE).

  • Arrays for specialised routines (e.g., (ppm, spectrum) for nmr, MARY tensors for lock-in models, steady-state vectors for linear solves).

Notes

  • Spaces and shapes: Propagation may occur in Liouville space; helpers such as _square_liouville_rhos are used to reshape back to Hilbert (N, N) for inspection. See each function’s docstring for shapes.

  • Units: Field parameters are mT unless otherwise specified; NMR chemical shifts are handled via Hz/MHz → ppm conversions; oop_eseem uses rad/s.

  • Performance: Inner loops precompute unit-field Zeeman operators and reuse sparse CSC matrices for repeated exponentials where possible.

  • CIDNP models: cidnp_model ∈ {“a”,”b”,”c”} selects exponential, truncated diffusion, or full diffusion variants and requires ks or alpha accordingly.

  • Lock-in MARY: modulated_mary_brute_force performs phase randomisation, envelope weighting, and numerical integration to estimate RMS harmonic responses.

References

Requirements:

  • numpy, scipy (sparse algebra, matrix exponentials), and a simulation object implementing APIs such as zeeman_hamiltonian, exchange_hamiltonian, dipolar_hamiltonian, hyperfine_hamiltonian, projection_operator, convert, time_evolution, product_probability, product_yield, and modifier application hooks.

See also

  • plot.py for visualization helpers.

  • kinetics.py, relaxation.py for incoherent processes added to Liouvillians.

  • Individual function docstrings within this module for precise signatures, units, and return structures.

radicalpy.experiments.anisotropy(sim, init_state, obs_state, time, theta, phi, B0, D, J, kinetics=[], relaxations=[])[source]

Anisotropy experiment.

Parameters:

  • init_state (State) – Initial State of the density matrix.

  • obs_state (State) – Observable State of the density matrix.

  • time (np.ndarray) – An sequence of (uniform) time points, usually created using np.arange or np.linspace.

  • theta (np.ndarray) – rotation (polar) angle between the external magnetic field and the fixed molecule. See zeeman_hamiltonian_3d.

  • phi (np.ndarray) – rotation (azimuth) angle between the external magnetic field and the fixed molecule. See zeeman_hamiltonian_3d.

  • B0 (float) – External magnetic field intensity (milli Tesla) (see zeeman_hamiltonian).

  • D (np.ndarray) – Dipolar coupling constant (see dipolar_hamiltonian).

  • J (float) – Exchange coupling constant (see exchange_hamiltonian).

  • kinetics (list) – A list of kinetic (super)operators of type radicalpy.kinetics.HilbertKineticsBase or radicalpy.kinetics.LiouvilleKineticsBase.

  • relaxations (list) – A list of relaxation superoperators of type radicalpy.relaxation.LiouvilleRelaxationBase.

Returns:

  • time: the original time object

  • B0: B0 parameter

  • theta: theta parameter

  • phi: phi parameter

  • rhos: tensor of sequences of time evolution of density

    matrices

  • time_evolutions: product probabilities

  • product_yields: product yields

  • product_yield_sums: product yield sums

Return type:

dict

radicalpy.experiments.anisotropy_loop(sim, init_state, obs_state, time, B0, H_base, theta, phi)[source]

Inner loop of anisotropy experiment.

Parameters:

  • sim (LiouvilleSimulation) – Simulation object.

  • init_state (State) – Initial State of the density matrix.

  • obs_state (State) – Observable State of the density matrix.

  • time (np.ndarray) – An sequence of (uniform) time points, usually created using np.arange or np.linspace.

  • B0 (float) – External magnetic field intensity (milli Tesla) (see zeeman_hamiltonian).

  • H_base (np.ndarray) – A “base” Hamiltonian, i.e., the Zeeman Hamiltonian will be added to this base, usually obtained with total_hamiltonian and B0=0.

  • theta (np.ndarray) – rotation (polar) angle between the external magnetic field and the fixed molecule. See zeeman_hamiltonian_3d.

  • phi (np.ndarray) – rotation (azimuth) angle between the external magnetic field and the fixed molecule. See zeeman_hamiltonian_3d.

Returns:

A tensor which has a series of density matrices for each angle theta and phi obtained by running time_evolution for each of them (with time time-steps, B0 magnetic intensity).

Return type:

np.ndarray

radicalpy.experiments.cidnp(B0, deltag, cidnp_model, nucleus_of_interest, donor_hfc_spinhalf, acceptor_hfc_spinhalf, donor_hfc_spin1, acceptor_hfc_spin1, ks=None, alpha=None)[source]

CIDNP polarisation vs field for a radical pair with S-T0 mixing only.

Parameters:

  • B0 (ndarray) – External magnetic field (T).

  • deltag (float) – Difference in g-value between the acceptor and donor.

  • cidnp_model (str) – Choose between CIDNP kinetic models. a) Exponential model. b) Truncated diffusion model. c) Full diffusion model.

  • ks (Optional[float]) – Decay rate constant for the Exponential model (1/s).

  • alpha (Optional[float]) – Parameter for the full diffusion model.

  • nucleus_of_interest (int) – The nucleus chosen for the simulation.

  • donor_hfc_spinhalf (float) – spin 1/2 HFCs (1H) for the donor (mT).

  • acceptor_hfc_spinhalf (float) – spin 1/2 HFCs (1H) for the acceptor (mT).

  • donor_hfc_spin1 (float) – spin 1 HFCs (14N) for the donor (mT).

  • acceptor_hfc_spin1 (float) – spin 1 HFCs (14N) for the acceptor (mT).

Return type:

(ndarray, ndarray)

Returns:

B0 (T) and polarisation (polarisation at each field point)

radicalpy.experiments.coherent_control(sim, *, init_state, obs_state, time, sticks_A_freq, sticks_A_int, sticks_B_freq, sticks_B_int, B1_G=200.0, g_e=2.0023, k_s=2000000.0, J=0.0, dt_override=None, u_max_factor=1.0, u_smooth=0.2)[source]

Coherent microwave-feedback simulation.

It treats the radical pair as an ensemble of frequency-configurations. Each configuration (A-stick × B-stick) is propagated in Hilbert space under

\[H_{ij}(t) = \omega_{A,i} S_{zA} + \omega_{B,j} S_{zB} + H_J + u(t)(S_{xA}+S_{xB})\]

where

  • the frequencies \(\omega_{A,i}\), \(\omega_{B,j}\) come from the stick lists sticks_A_freq and sticks_B_freq (in rad/s),

  • the weights of each configuration are the products of the normalised stick intensities sticks_A_int and sticks_B_int,

  • u(t) is a global feedback field, common to all configurations, built from the ensemble signal.

The feedback field is computed at every time step as

\[u_\text{raw}(t) = B_{1c}\, \sum_{ij} w_{ij} \, \Im \operatorname{Tr} \left[ P_\text{target}( V \rho_{ij} - \rho_{ij} V ) \right],\]

then clipped to ±u_max_factor * B1c and optionally smoothed in time by an exponential smoother. Here

  • V = S_xA + S_xB is the microwave operator,

  • P_target = |αβ⟩⟨αβ| + |βα⟩⟨βα| is the “target” projector,

  • B1c = g_e μ_B B1_T / ħ is the on-resonance Rabi frequency, with B1_G given in gauss.

Time propagation is done configuration-by-configuration with an RK4 step on the user-supplied time grid.

Parameters:

  • sim – HilbertSimulation

  • init_state – State Spin state used both for the initial density.

  • obs_state – State Spin state to monitor as an observable population during the run. A projector onto this state is evaluated on every configuration and then averaged over the ensemble.

  • time (ndarray) – ndarray 1D array of time points (s). The integration advances along this grid. If dt_override is not given, the step size is inferred from time[1] - time[0].

  • sticks_A_freq (Sequence[float]) – sequence of float Stick frequencies for radical A and radical B, respectively, in angular frequency units (rad/s), e.g. 2*pi*1e6*[ ... MHz ... ].

  • sticks_B_freq (Sequence[float]) – sequence of float Stick frequencies for radical A and radical B, respectively, in angular frequency units (rad/s), e.g. 2*pi*1e6*[ ... MHz ... ].

  • sticks_A_int (Sequence[float]) – sequence of float Corresponding intensities for the two stick sets. Each list is normalised so that its entries sum to 1. The weight of configuration (i, j) is then the product of the two normalised intensities.

  • sticks_B_int (Sequence[float]) – sequence of float Corresponding intensities for the two stick sets. Each list is normalised so that its entries sum to 1. The weight of configuration (i, j) is then the product of the two normalised intensities.

  • B1_G (float) – float, optional Microwave amplitude in gauss. Default is 200.0.

  • g_e (float) – float, optional Electron g-value used in the Rabi-frequency conversion. Default 2.0023.

  • k_s (float) – float, optional Haberkorn recombination rate applied with respect to init_state (s⁻¹). Default 2e6.

  • J (float) – float, optional Exchange coupling (mT).

  • dt_override (Optional[float]) – float, optional If provided, use this as the integration time step instead of time[1] - time[0].

  • u_max_factor (float) – float, optional Hard clip on the feedback: |u(t)| ≤ u_max_factor * B1c. Useful to suppress large initial spikes. Default 1.0.

  • u_smooth (float) – float, optional Exponential smoothing factor for the feedback in (0, 1). Larger values give smoother control fields. Default 0.2.

Return type:

Dict[str, ndarray]

Returns:

dict of str -> ndarray

A dictionary with the following entries:

  • "time": Time axis (s)

  • "u": Feedback field (rad/s)

  • "target": nsemble-average target population

  • "each_target": Target population per configuration

  • "weights": Configuration weights

  • "obs": Ensemble-average population of obs_state

  • "each_obs": Per-configuration population of obs_state

  • "population": Ensemble-average trace of ρ

  • "each_population": Per-configuration trace

radicalpy.experiments.epr(sim, init_state, obs_state, time, D, J, B0, B1, B1_freq, B0_axis='z', B1_axis='x', kinetics=[], relaxations=[], hfc_anisotropy=False)[source]

Electron paramagnetic resonance (EPR) time‐domain simulation with CW/AC fields.

Parameters:

  • sim (HilbertSimulation) – Simulation object.

  • init_state (State) – Initial State of the density matrix.

  • obs_state (State) – Observable State of the density matrix.

  • time (np.ndarray) – Sequence of (uniform) time points.

  • D (float) – Dipolar coupling constant (mT).

  • J (float) – Exchange coupling constant (mT).

  • B0 (np.ndarray) – Static magnetic field sweep (mT) applied along B0_axis.

  • B1 (float) – RF/AC field amplitude (mT) applied along B1_axis.

  • B1_freq (float) – RF/AC angular frequency offset (mT). Implemented as an effective Zeeman term with negative sign on B0_axis.

  • B0_axis (str) – Axis for B0 Zeeman term ('x'|'y'|'z').

  • B1_axis (str) – Axis for B1 Zeeman term ('x'|'y'|'z').

  • kinetics (list) – List of kinetic superoperators.

  • relaxations (list) – List of relaxation superoperators.

  • hfc_anisotropy (bool) – Include anisotropic hyperfine Hamiltonian if True.

Returns:

  • time: original time

  • B0: sweep values

  • B0_axis: axis of B0

  • B1: RF amplitude

  • B1_axis: axis of B1

  • B1_freq: RF angular frequency offset

  • B1_freq_axis: axis used for the offset term

  • rhos: density matrices (squared to Hilbert shape)

  • time_evolutions: product probabilities vs time & field

  • product_yields: integrated product yields

  • product_yield_sums: scalar yield per field

  • MARY: normalised magnetoresponse

  • LFE: low field effect (%)

  • HFE: high field effect (%)

Return type:

dict

radicalpy.experiments.field_switching(sim, *, B_on=2.0, B_off=0.0, init_state=State.TRIPLET, dt=1e-09, n_offsets=100, offset_step=10, pulse_width_steps=800, k_rec=50000000.0, k_esc=2000000.0, J=0.0, D=0.0)[source]

Nanosecond field-switching experiment. Switched External Magnetic Field (SEMF) simulation.

Parameters:

  • sim (HilbertSimulation) – Sim object.

  • B_on (float) – Magnetic fields (mT) for the “before switch” part and the “after switch” part, respectively.

  • B_off (float) – Magnetic fields (mT) for the “before switch” part and the “after switch” part, respectively.

  • dt (float) – Time step (s).

  • n_offsets (int) – Number of different switch times to scan.

  • offset_step (int) – Number of time steps to increase the switch time by each iteration.

  • pulse_width_steps (int) – Length (in steps) of the second part of the evolution (the “pulse”).

  • k_rec (float) – Haberkorn recombination and escape rates (s⁻¹).

  • k_esc (float) – Haberkorn recombination and escape rates (s⁻¹).

  • J (float) – Exchange and dipolar couplings.

  • D (float) – Exchange and dipolar couplings.

Returns:

  • time: (T,),

  • switch_times: (n_offsets,),

  • TA_on: (T, n_offsets),

  • TA_off: (T, n_offsets),

  • TA_diff: (T, n_offsets),

Return type:

A dict with the following keys (dims)

radicalpy.experiments.kine_quantum_mary(sim, num_samples, init_state, radical_pair, ts, Bs, D, J, kinetics, relaxations)[source]

Kinetic + quantum MARY (hybrid) with sampling over stochastic hyperfine fields.

Source: Antill and Vatai, J. Chem. Theory Comput. 20, 21, 9488–9499 (2024).

Parameters:

  • sim (SemiclassicalSimulation) – Simulation object.

  • num_samples (int) – Number of stochastic hyperfine field realisations.

  • init_state (ArrayLike) – Initial population vector/state.

  • radical_pair (list) – Slice indices [i0, i1] defining the radical-pair subspace into which the Liouvillian is inserted.

  • ts (np.ndarray) – Time grid (s).

  • Bs (ArrayLike) – Magnetic field sweep values (mT).

  • D (float) – Dipolar coupling constant (mT).

  • J (float) – Exchange coupling constant (mT).

  • kinetics (ArrayLike) – Kinetic matrix to be added to the Liouvillian.

  • relaxations (list[ArrayLike]) – List of relaxation superoperators.

Returns:

  • ts: time grid

  • Bs: field sweep

  • yield: complex yield tensor with shape (len(ts), len(init_state), len(Bs))

Return type:

dict

radicalpy.experiments.magnetic_field_loop(sim, init_state, time, H_base, B, B_axis, theta=None, phi=None, hfc_anisotropy=False)[source]

Generate density matrices (rhos) for MARY.

Parameters:

  • sim (HilbertSimulation) – Simulation object.

  • init_state (State) – Initial State of the density matrix.

  • time (np.ndarray) – A sequence of (uniform) time points, usually created using np.arange or np.linspace (s).

  • H_base (np.ndarray) – A “base” Hamiltonian or Liouvillian that does not include the swept Zeeman term (e.g. the result of total_hamiltonian(...) with B0=0 and any static terms).

  • B (np.ndarray) – Magnetic-field sweep values (mT). Each value scales a precomputed unit-field Zeeman operator on B_axis.

  • B_axis (str) – Axis of the Zeeman term ('x'|'y'|'z').

  • theta (float, optional) – Rotation (polar) angle between the external magnetic field and the fixed molecular frame. See zeeman_hamiltonian_3d.

  • phi (float, optional) – Rotation (azimuth) angle between the external magnetic field and the fixed molecular frame. See zeeman_hamiltonian_3d.

  • hfc_anisotropy (bool) – Reserved for API symmetry; anisotropic hyperfine terms should be included in H_base if required (not used here).

Returns:

A tensor of density matrices with shape (len(B), len(time), *rho_shape), where rho_shape is (N, N) for Hilbert-space propagation and (N,) for Liouville-space propagation. For each field value B[i], the time evolution is obtained by propagating under H = H_base + B[i] * H_Z(1.0) with H_Z(1.0) the unit-field Zeeman operator on B_axis.

Return type:

np.ndarray

Notes

  • The Zeeman operator for a unit field is computed once and reused for all B[i], then scaled and added to H_base.

  • Propagation uses a CSC sparse representation for efficiency.

radicalpy.experiments.magnetic_field_loop_semiclassical(sim, init_state, time, B, H_base, kinetics=[], relaxations=[], theta=None, phi=None, num_samples=None)[source]

Generate density matrices (rhos) for MARY using a semiclassical ensemble of hyperfine realisations and incoherent processes.

Parameters:

  • sim (HilbertSimulation) – Simulation object.

  • init_state (State) – Initial State of the density matrix.

  • time (np.ndarray) – A sequence of (uniform) time points, usually created using np.arange or np.linspace (s).

  • B (np.ndarray) – Magnetic-field sweep values (mT). Each value scales a precomputed unit-field Zeeman operator summed over all radicals on the laboratory 'z' axis (after applying the optional theta/phi rotation).

  • H_base (np.ndarray) – A “base” Hamiltonian that does not include the swept Zeeman term, e.g., the result of total_hamiltonian(...) with B0=0 plus any static terms.

  • kinetics (list[HilbertIncoherentProcessBase]) – List of incoherent kinetic processes (e.g., recombination, intersystem crossing) to be added as Liouvillian modifiers. Defaults to [].

  • relaxations (list[HilbertIncoherentProcessBase]) – List of relaxation/dephasing processes to be added as Liouvillian modifiers. Defaults to [].

  • theta (float, optional) – Rotation (polar) angle between the external magnetic field and the fixed molecular frame. See zeeman_hamiltonian_3d.

  • phi (float, optional) – Rotation (azimuth) angle between the external magnetic field and the fixed molecular frame. See zeeman_hamiltonian_3d.

  • num_samples (int, optional) – Number of semiclassical hyperfine realizations used for ensemble averaging. If None, the simulator’s default is used by sim.semiclassical_HHs(...).

Returns:

A tensor of density matrices with shape (len(B), len(time), *rho_shape), where rho_shape is (N, N) for Hilbert-space propagation and (N,) for Liouville-space propagation. For each field value B[i], the result is the sample-average over num_samples semiclassical Hamiltonians H_HH[k] of the time evolution under

H_t = H_base + B[i] * H_Z(1.0) + H_HH[k],

with incoherent kinetics + relaxations applied as Liouvillian modifiers.

Return type:

np.ndarray

Notes

  • The unit-field Zeeman operator is constructed by summing the per-radical contributions on the 'z' axis (after rotation by theta/phi) and is then scaled by each B[i].

  • An ensemble HHs = sim.semiclassical_HHs(num_samples) provides the semiclassical hyperfine Hamiltonians used for averaging.

  • sim.convert(...) maps the total Hamiltonian to the working representation (Hilbert→Liouville if required); incoherent processes are then added via sim.apply_liouville_hamiltonian_modifiers(...).

  • Propagation uses a CSC sparse representation for efficiency, and the final array at index i is the mean across all samples for that field value.

radicalpy.experiments.mary(sim, init_state, obs_state, time, B, D, J, kinetics=[], relaxations=[], theta=None, phi=None, hfc_anisotropy=False)[source]

Magnetically affected reaction yield (MARY) simulation over a magnetic field sweep.

Parameters:

  • sim (HilbertSimulation) – Simulation object.

  • init_state (State) – Initial State of the density matrix.

  • obs_state (State) – Observable State of the density matrix.

  • time (np.ndarray) – Sequence of (uniform) time points (s).

  • B (np.ndarray) – Magnetic field sweep values (mT) along z-axis by default.

  • D (float) – Dipolar coupling constant (mT).

  • J (float) – Exchange coupling constant (mT).

  • kinetics (list) – List of kinetic superoperators.

  • relaxations (list) – List of relaxation superoperators.

  • theta (float, optional) – Polar angle for the Zeeman term (anisotropy).

  • phi (float, optional) – Azimuthal angle for the Zeeman term (anisotropy).

  • hfc_anisotropy (bool) – Include anisotropic hyperfine Hamiltonian if True.

Returns:

  • time: original time

  • B: sweep values

  • theta: theta parameter

  • phi: phi parameter

  • rhos: density matrices (squared to Hilbert shape)

  • time_evolutions: product probabilities

  • product_yields: product yields

  • product_yield_sums: product yield sums

  • MARY: normalised magnetoresponse

  • LFE: low field effect (%)

  • HFE: high field effect (%)

Return type:

dict

radicalpy.experiments.mary_lfe_hfe(obs_state, B, product_probability_seq, dt, k, c=0.0)[source]

Calculate MARY, LFE, HFE.

Return type:

(ndarray, ndarray, ndarray)

radicalpy.experiments.mary_semiclassical(sim, init_state, obs_state, time, B, D, J, kinetics=[], relaxations=[], theta=None, phi=None, num_samples=None, c=None)[source]

Compute a MARY (Magnetically Affected Reaction Yield) curve using the semiclassical simulation pipeline.

This routine constructs a field-independent base Hamiltonian \(H = H_J + H_D\) from the exchange and dipolar interactions, then performs a semiclassical magnetic-field sweep, evolving the state from init_state over the time grid for each field value in B. Product probabilities are converted to integrated yields and post-processed into the standard MARY, low-field effect (LFE), and high-field effect (HFE) metrics.

Parameters:

  • sim (SemiclassicalSimulation) – SemiclassicalSimulation Semiclassical simulator providing Hamiltonian terms and evolution back-ends (including stochastic sampling where applicable).

  • init_state (State) – State Initial state for time evolution (e.g., State.SINGLET).

  • obs_state (State) – State Observable (product) state whose probability/yield is reported.

  • time (ndarray) – ndarray, shape (T,) Monotonic time grid for propagation (units consistent with Hamiltonians).

  • B (ndarray) – ndarray, shape (NB,) Magnetic-field grid for the sweep (scalar magnitude in the chosen axis).

  • D (float) – float Dipolar interaction strength used to build H_D.

  • J (float) – float Exchange interaction strength used to build H_J.

  • kinetics (list[HilbertIncoherentProcessBase]) – list[HilbertIncoherentProcessBase], optional Kinetic processes to apply during evolution (e.g., Haberkorn).

  • relaxations (list[HilbertIncoherentProcessBase]) – list[HilbertIncoherentProcessBase], optional Relaxation processes to apply during evolution (semiclassical variants).

  • theta (Optional[float]) – float, optional Polar angle (radians) of the magnetic-field axis relative to the lab frame. If None, the simulator default is used.

  • phi (Optional[float]) – float, optional Azimuthal angle (radians) of the magnetic-field axis. If None, the simulator default is used.

  • num_samples (Optional[int]) – int, optional Number of stochastic samples/trajectories for the semiclassical averaging. If None, the simulator default is used.

  • c (Optional[float]) – float, optional Normalisation/contrast parameter forwarded to mary_lfe_hfe for MARY post-processing.

Return type:

dict

Returns:

dict

A results dictionary with keys: - time : ndarray, the input time grid. - B : ndarray, the field grid. - theta / phi : angles used for the sweep. - rhos : list/ndarray of density matrices (reshaped to square form). - time_evolutions : ndarray, product probabilities vs. time and B. - product_yields : ndarray, integrated yields vs. B. - product_yield_sums : float or ndarray, sum. - MARY : ndarray, normalised magnetoresponse vs. B. - LFE : float, low-field effect (%). - HFE : float, high-field effect (%).

radicalpy.experiments.modulated_mary_brute_force(Bs, modulation_depths, modulation_frequency, time_constant, harmonics, lfe_magnitude)[source]

Lock-in detected MARY via brute-force phase randomisation and numerical integration.

Source: Konowalczyk et al. Phys. Chem. Chem. Phys. 23, 1273-1284 (2021).

Parameters:

  • Bs (np.ndarray) – Array of bias fields at which to compute the MARY signal (G).

  • modulation_depths (list) – List of field modulation amplitudes (G).

  • modulation_frequency (float) – Field modulation angular frequency (Hz).

  • time_constant (float) – Lock-in amplifier time constant used for exponential weighting of the reference (s).

  • harmonics (list) – Harmonic indices (e.g., [1,2]) for which to compute responses.

  • lfe_magnitude (float) – Amplitude parameter controlling the Lorentzian LFE line shape.

Returns:

Tensor of shape (len(harmonics), len(modulation_depths), len(Bs)) containing RMS lock-in amplitudes at each harmonic, modulation depth, and bias field.

Return type:

np.ndarray

radicalpy.experiments.nmr(multiplets, spectral_width, number_of_points, fft_number, transmitter_frequency, carrier_position, linewidth, scale=1.0)[source]

Simple 1D NMR spectrum synthesiser (FID → FFT) with multiplets and T2 decay.

Parameters:

  • multiplets (list) – Each entry [number of nuclei, chemical shift, multiplicity, scalar coupling].

  • spectral_width (float) – Spectral width (Hz).

  • number_of_points (float) – Number of acquired points in the time domain.

  • fft_number (float) – Zero-filled length for the FFT (≥ number_of_points).

  • transmitter_frequency (float) – Transmitter frequency (MHz).

  • carrier_position (float) – Carrier position (ppm) used to define the reference.

  • linewidth (float) – Lorentzian linewidth (Hz) → T2 = 1/(π·linewidth).

  • scale (float) – Overall multiplicative scale on the final spectrum.

Returns:

  • ppm: Chemical shift axis (ppm).

  • spectrum: Complex spectrum after FFT (same length as fft_number).

Return type:

(np.ndarray, np.ndarray)

radicalpy.experiments.odmr(sim, init_state, obs_state, time, D, J, B0, B1, B1_freq, B0_axis='z', B1_axis='x', kinetics=[], relaxations=[], hfc_anisotropy=False)[source]

Optically detected magnetic resonance (ODMR) simulation vs RF frequency.

Parameters:

  • sim (HilbertSimulation) – Simulation object.

  • init_state (State) – Initial State of the density matrix.

  • obs_state (State) – Observable State of the density matrix.

  • time (np.ndarray) – Sequence of (uniform) time points (s).

  • D (float) – Dipolar coupling constant (mT).

  • J (float) – Exchange coupling constant (mT).

  • B0 (float) – Static magnetic field (mT) along B0_axis.

  • B1 (float) – RF/AC field amplitude (mT) along B1_axis.

  • B1_freq (np.ndarray) – RF angular frequency sweep (mT).

  • B0_axis (str) – Axis for B0 Zeeman term.

  • B1_axis (str) – Axis for B1 Zeeman term.

  • kinetics (list) – List of kinetic superoperators.

  • relaxations (list) – List of relaxation superoperators.

  • hfc_anisotropy (bool) – Include anisotropic hyperfine Hamiltonian if True.

Returns:

  • time: original time

  • B0: static field value

  • B0_axis: axis of B0

  • B1: RF amplitude

  • B1_axis: axis of B1

  • B1_freq: RF sweep values

  • B1_freq_axis: axis used for the RF sweep

  • rhos: density matrices (squared to Hilbert shape)

  • time_evolutions: product probabilities vs time & field

  • product_yields: integrated product yields

  • product_yield_sums: scalar yield per field

  • MARY: normalised magnetoresponse

  • LFE: low field effect (%)

  • HFE: high field effect (%)

Return type:

dict

radicalpy.experiments.omfe(sim, init_state, obs_state, time, D, J, B1, B1_freq, B1_axis='x', B1_freq_axis='z', kinetics=[], relaxations=[], hfc_anisotropy=False)[source]

Oscillating magnetic field effect (OMFE) simulation vs RF frequency in transverse field.

Parameters:

  • sim (HilbertSimulation) – Simulation object.

  • init_state (State) – Initial State of the density matrix.

  • obs_state (State) – Observable State of the density matrix.

  • time (np.ndarray) – Sequence of (uniform) time points (s).

  • D (float) – Dipolar coupling constant (mT).

  • J (float) – Exchange coupling constant (mT).

  • B1 (float) – Oscillating field amplitude (mT) along B1_axis.

  • B1_freq (np.ndarray) – Angular frequency sweep (mT) applied along B1_freq_axis.

  • B1_axis (str) – Axis for the static B1 Zeeman term.

  • B1_freq_axis (str) – Axis used for the frequency-sweep effective term.

  • kinetics (list) – List of kinetic superoperators.

  • relaxations (list) – List of relaxation superoperators.

  • hfc_anisotropy (bool) – Include anisotropic hyperfine Hamiltonian if True.

Returns:

  • time: original time

  • B1: RF amplitude

  • B1_axis: axis of B1

  • B1_freq: RF sweep values

  • B1_freq_axis: axis used for the RF sweep

  • rhos: density matrices (squared to Hilbert shape)

  • time_evolutions: product probabilities vs time & field

  • product_yields: integrated product yields

  • product_yield_sums: scalar yield per field

  • MARY: normalised magnetoresponse

  • LFE: low field effect (%)

  • HFE: high field effect (%)

Return type:

dict

radicalpy.experiments.oop_eseem(tau, J, D, T1=inf, n_quad=200)[source]

Out-of-phase-electron-spin echo envelope modulation (OOP-ESEEM) simulation. Computes S(tau) ∝ exp(-tau/T1) * ∫_0^π sin( 2 [ J - D (cos^2 θ - 1/3) ] tau ) sinθ dθ using Gauss–Legendre quadrature on u = cosθ ∈ [-1, 1].

Parameters:

  • tau (float or np.ndarray) – Time (s).

  • J (float) – Exchange interaction (rad/s).

  • D (float) – Dipolar coupling constant (rad/s).

  • T1 (float) – Longitudinal relaxation time constant (s). Use np.inf to disable the exponential decay.

  • n_quad (int) – Number of Gauss–Legendre nodes (accuracy increases with n).

Returns:

OOP-ESEEM spectrum.

Return type:

S (np.ndarray)

radicalpy.experiments.rydmr(sim, init_state, obs_state, time, D, J, B0, B1, B1_freq, B0_axis='z', B1_axis='x', kinetics=[], relaxations=[], hfc_anisotropy=False)[source]

Reaction yield-detected magnetic resonance (RYDMR) vs static field.

Parameters:

  • sim (HilbertSimulation) – Simulation object.

  • init_state (State) – Initial State of the density matrix.

  • obs_state (State) – Observable State of the density matrix.

  • time (np.ndarray) – Sequence of (uniform) time points (s).

  • D (float) – Dipolar coupling constant (mT).

  • J (float) – Exchange coupling constant (mT).

  • B0 (np.ndarray) – Static field sweep (mT) along B0_axis.

  • B1 (float) – RF/AC amplitude (mT) along B1_axis.

  • B1_freq (float) – RF angular frequency offset (mT).

  • B0_axis (str) – Axis for B0 Zeeman term.

  • B1_axis (str) – Axis for B1 Zeeman term.

  • kinetics (list) – List of kinetic superoperators.

  • relaxations (list) – List of relaxation superoperators.

  • hfc_anisotropy (bool) – Include anisotropic hyperfine Hamiltonian if True.

Returns:

  • time: original time

  • B0: sweep values

  • B0_axis: axis of B0

  • B1: RF amplitude

  • B1_axis: axis of B1

  • B1_freq: RF angular frequency offset

  • B1_freq_axis: axis used for the offset term

  • rhos: density matrices (squared to Hilbert shape)

  • time_evolutions: product probabilities vs time & field

  • product_yields: integrated product yields

  • product_yield_sums: scalar yield per field

  • MARY: normalised magnetoresponse

  • LFE: low field effect (%)

  • HFE: high field effect (%)

Return type:

dict

radicalpy.experiments.semiclassical_mary(sim, num_samples, init_state, ts, Bs, D, J, triplet_excited_state_quenching_rate, free_radical_escape_rate, kinetics, relaxations, scale_factor)[source]

Semiclassical MARY with stochastic averaging and explicit population channels.

Source: Maeda et al. Mol. Phys. 104, 1779–1788 (2006).

Parameters:

  • sim (SemiclassicalSimulation) – Simulation object.

  • num_samples (int) – Number of stochastic hyperfine field realisations.

  • init_state (State) – Initial projection operator state.

  • ts (ArrayLike) – Time grid (s).

  • Bs (ArrayLike) – Magnetic field sweep (mT).

  • D (float) – Dipolar coupling constant (mT).

  • J (float) – Exchange coupling constant (mT).

  • triplet_excited_state_quenching_rate (float) – Quenching rate (1/s) feeding RP.

  • free_radical_escape_rate (float) – Escape rate from free radical (1/s).

  • kinetics (list[ArrayLike]) – Kinetic Liouvillian/superoperators.

  • relaxations (list[ArrayLike]) – Relaxation superoperators.

  • scale_factor (float) – Averaging weights for decay construction.

Returns:

  • ts: time grid

  • Bs: field sweep

  • MARY: response matrix with shape (len(ts), len(Bs))

Return type:

dict

radicalpy.experiments.steady_state_mary(sim, init, obs, Bs, D, E, J, theta, phi, kinetics=[])[source]

Steady-state MARY via linear solve of Liouvillian with ZFS, exchange, and Zeeman terms.

Parameters:

  • sim (LiouvilleSimulation) – Simulation object.

  • obs (State) – Observable State used to form the projection operator Q.

  • Bs (np.ndarray) – Magnetic field sweep (mT).

  • D (float) – Zero-field splitting parameter D (mT).

  • E (float) – Zero-field splitting parameter E (mT).

  • J (float) – Exchange coupling (mT).

  • theta (float) – Polar angle of the Zeeman field.

  • phi (float) – Azimuthal angle of the Zeeman field.

  • kinetics (list) – List of kinetic superoperators applied to the Liouvillian.

Returns:

  • rhos: Steady-state density vectors, one per field (shape (len(Bs), N)).

  • Phi_s: Projection rhos @ Q.flatten() giving steady-state signal per field.

Return type:

np.ndarray

radicalpy.kinetics module

Kinetics operators and superoperators for radical-pair dynamics.

This module implements common incoherent kinetic processes acting on spin systems in either Hilbert space (operators) or Liouville space (superoperators). These terms model recombination, product formation, and phenomenological decay channels that act alongside coherent spin evolution.

Classes:

  • HilbertKineticsBase: Base class for Hilbert-space kinetics operators.

  • LiouvilleKineticsBase: Base class for Liouville-space kinetics superoperators.

  • Exponential: Simple exponential decay of product probabilities in Hilbert space (Kaptein model).

  • Haberkorn: Haberkorn singlet/triplet selective recombination superoperator.

  • HaberkornFree: Haberkorn-style uniform decay (free radical / RP2 formation).

  • JonesHore: Jones–Hore two-site kinetics with separate S/T rates.

Usage pattern:

  1. Instantiate a kinetics object with the desired rate parameters (e.g., rate_constant, singlet_rate, triplet_rate, target).

  2. Call init(sim) with a simulation exposing the required API (projection operators, sizes, and Liouville conversion).

  3. Combine the returned term with your total generator: - Hilbert space: use the operator’s effect where appropriate. - Liouville space: add subH to the Liouvillian.

Args conventions (per class):

  • Exponential(rate_constant): Scales product probabilities as exp(-rate * t).

  • Haberkorn(rate_constant, target): Selective S/T recombination where target is one of State.SINGLET, State.TRIPLET, or specific triplet substates supported by the code.

  • HaberkornFree(rate_constant): Uniform (state-independent) decay proportional to identity in Liouville space.

  • JonesHore(singlet_rate, triplet_rate): Two-channel model with separate S and T rates and a coupling term constructed from S/T projectors.

Notes

  • Hilbert vs Liouville: Hilbert-space Exponential modifies product probabilities directly; Liouville-space classes construct superoperators (subH) via Kronecker products and projector conversions.

  • Projectors: QS, QT, and triplet sublevel projectors are obtained from the simulation (sim.projection_operator(...)); Haberkorn and Jones–Hore rely on these for state selectivity.

  • Rates and units: All rate parameters are in s⁻¹. Ensure time arrays passed to Exponential.adjust_product_probabilities are in seconds.

References

Requirements:

  • numpy for array operations.

  • A simulation object with: projection_operator(State.*), hamiltonian_size, and a Liouville conversion method (e.g., sim._convert(Q)).

raises ValueError:

Haberkorn validates that target is singlet or triplet (or supported triplet sublevels) and will raise if unsupported.

See also

radicalpy.simulation.HilbertIncoherentProcessBase, radicalpy.simulation.LiouvilleIncoherentProcessBase, and related relaxation modules for non-kinetic incoherent processes.

class radicalpy.kinetics.Exponential(rate_constant)[source]

Bases: HilbertKineticsBase

Exponential model kinetics operator.

Source: Kaptein et al. Chem. Phys. Lett. 4, 4, 195-197 (1969).

>>> Exponential(rate_constant=1e6)
Kinetics: Exponential
Rate constant: 1000000.0
adjust_product_probabilities(product_probabilities, time)[source]

See radicalpy.simulation.HilbertIncoherentProcessBase.adjust_product_probabilities.

class radicalpy.kinetics.Haberkorn(rate_constant, target)[source]

Bases: LiouvilleKineticsBase

Haberkorn kinetics superoperator for singlet/triplet recombination/product formation.

Source: Haberkorn, Mol. Phys. 32:5, 1491-1493 (1976).

Parameters:

  • rate_constant (float) – The kinetic rate constant (1/s).

  • target (State) – The target state of the reaction pathway (singlet or triplet states only).

>>> Haberkorn(rate_constant=1e6, target=State.SINGLET)
Kinetics: Haberkorn
Rate constant: 1000000.0
Target: S

>>> Haberkorn(rate_constant=1e6, target=State.TRIPLET)
Kinetics: Haberkorn
Rate constant: 1000000.0
Target: T
init(sim)[source]

See radicalpy.simulation.HilbertIncoherentProcessBase.init.

class radicalpy.kinetics.HaberkornFree(rate_constant)[source]

Bases: LiouvilleKineticsBase

Haberkorn kinetics superoperator for free radical/RP2 formation.

Source: Haberkorn, Mol. Phys. 32:5, 1491-1493 (1976).

>>> HaberkornFree(rate_constant=1e6)
Kinetics: HaberkornFree
Rate constant: 1000000.0
init(sim)[source]

See radicalpy.simulation.HilbertIncoherentProcessBase.init.

class radicalpy.kinetics.HilbertKineticsBase(rate_constant)[source]

Bases: HilbertIncoherentProcessBase

Base class for kinetics operators (Hilbert space).

class radicalpy.kinetics.JonesHore(singlet_rate, triplet_rate)[source]

Bases: LiouvilleKineticsBase

Jones-Hore kinetics superoperator for two-site models.

Source: Jones et al. Chem. Phys. Lett. 507, 269-273 (2011).

Parameters:

  • singlet_rate (float) – Singlet recombination rate constant (1/s).

  • triplet_rate (float) – Triplet product formation rate constant (1/s).

>>> JonesHore(1e7, 1e6)
Kinetics: JonesHore
Singlet rate: 10000000.0
Triplet rate: 1000000.0
init(sim)[source]

See radicalpy.simulation.HilbertIncoherentProcessBase.init.

property rate_constant: float

Average rate of the kinetic processes.

class radicalpy.kinetics.LiouvilleKineticsBase(rate_constant)[source]

Bases: LiouvilleIncoherentProcessBase

Base class for kinetics superoperators (Liouville space).

radicalpy.plot module

Plotting utilities for spin dynamics, tensors, and Monte Carlo trajectories.

This module collects convenience routines for visualising common objects in spin chemistry and molecular simulations, including: energy-level diagrams, anisotropy/tensor surfaces on the unit sphere, 3D random-walk trajectories (free and caged), density-matrix bar animations, and general helper plots for molecules and time–series.

Functions in this module build on NumPy and Matplotlib’s 3D toolkit and are intended for exploratory analysis and figure generation inside notebooks or scripts. Most functions draw directly to the active Matplotlib figure and return None, except where noted.

Contents:

Notes

  • Many plotting functions assume Matplotlib 3D axes; some create and manage figures/axes internally. To integrate into existing layouts, adapt the code to accept/use a provided Axes3D where appropriate.

  • Units are noted in axis labels (e.g., nm, mT, s). Some helpers apply scaling factors (e.g., 1e9) for visualization; inspect the source if absolute coordinates are required.

  • energy_levels assumes a HilbertSimulation API exposing total_hamiltonian and zeeman_hamiltonian. Replace or wrap as needed.

Dependencies:

  • numpy for array operations.

  • matplotlib (including mpl_toolkits.mplot3d, cm, and colors) for plotting.

  • matplotlib.animation.FuncAnimation for animated density-matrix plots.

raises ValueError:

spin_state_labels raises if sim does not contain exactly two radicals (labels are defined for radical pairs).

See also

matplotlib.pyplot, mpl_toolkits.mplot3d, and domain-specific simulation objects (e.g., HilbertSimulation) used by the energy-level routines.

radicalpy.plot.anisotropy_surface(theta, phi, Y)[source]

Anisotropy surface plot.

Renders a 3D surface on the unit sphere whose radius is proportional to the real part of a field Y(θ,φ), and colors the surface by Re(Y).

Parameters:

  • theta (np.ndarray) – Polar angles θ in radians (shape: (N,) or mesh-ready).

  • phi (np.ndarray) – Azimuthal angles φ in radians (shape: (M,) or mesh-ready).

  • Y (np.ndarray) – Complex field values sampled on (θ, φ); only the real component is used for coloring and scaling.

Returns:

Displays a Matplotlib 3D surface plot.

Return type:

None

radicalpy.plot.density_matrix_animation(rhos, frames, bar3d_kwargs, axes_kwargs)[source]

Density matrix bar-plot animation.

Builds a 3D bar plot over time from a sequence of density matrices, coloring bars by magnitude and returning a FuncAnimation.

Parameters:

  • rhos (Sequence[np.ndarray]) – Time-ordered density matrices ρ_t (each shape: (N, N)).

  • frames (int) – Number of frames to render in the animation.

  • bar3d_kwargs (dict) – Keyword arguments forwarded to Axes3D.bar3d.

  • axes_kwargs (dict) – Keyword arguments forwarded to Axes3D.set.

Returns:

The configured animation object.

Return type:

matplotlib.animation.FuncAnimation

radicalpy.plot.energy_levels(sim, B, J=0, D=0)[source]

Energy levels vs magnetic field.

Computes eigen-energies of the total Hamiltonian H = H_base + B0 * H_zee for each applied field in B and plots levels as functions of B0.

Parameters:

  • sim (HilbertSimulation) – Simulation in Hilbert space.

  • B (np.ndarray) – Magnetic field values (T).

  • J (float) – Isotropic exchange parameter (default: 0).

  • D (float) – Dipolar coupling parameter (default: 0).

Returns:

Displays a Matplotlib line plot of energy levels vs field.

Return type:

None

radicalpy.plot.linear_energy_levels(H, B, linecolour, title)[source]

Linear spectrum event plot.

Diagonalizes a Hamiltonian and shows its eigenvalues as vertical events, useful for compact level-spectrum visualization.

Parameters:

  • H (np.ndarray) – Hamiltonian matrix (Hermitian).

  • B (Any) – Unused placeholder (kept for API compatibility).

  • linecolour (str or tuple) – Color spec for event lines.

  • title (str) – Plot title.

Returns:

Displays a vertical event plot of eigenvalues.

Return type:

None

radicalpy.plot.monte_carlo_caged(pos, r_max)[source]

3D random-walk trajectory (caged).

Plots a Monte Carlo trajectory inside a spherical cage of radius r_max, rendering the wireframe boundary and the particle path.

Parameters:

  • pos (np.ndarray) – Cartesian positions along a trajectory (shape: (T, 3), in metres).

  • r_max (float) – Cage radius (metres) for the confining sphere.

Returns:

Displays a 3D wireframe and trajectory plot.

Return type:

None

radicalpy.plot.monte_carlo_free(pos)[source]

3D random-walk trajectory (free solution).

Plots a radical-pair Monte Carlo path in free solution, scaling coordinates to nanometres for visualization and marking start and origin.

Parameters:

pos (np.ndarray) – Cartesian positions along a trajectory (shape: (T, 3), in metres).

Returns:

Displays a 3D line plot of the walk.

Return type:

None

radicalpy.plot.plot_3d_results(xdata, ydata, zdata, xlabel, ylabel, zlabel, azim=-135, dist=10, elev=35, factor=1000000.0)[source]

3D surface plot of results.

Creates a colored surface z(x,y) with adjustable 3D view parameters and axis labels—useful for field maps or parameter sweeps.

Parameters:

  • xdata (np.ndarray) – X-axis coordinates.

  • ydata (np.ndarray) – Y-axis coordinates.

  • zdata (np.ndarray) – Z values on the meshgrid of (x, y).

  • xlabel (str) – Label for the X axis.

  • ylabel (str) – Label for the Y axis.

  • zlabel (str) – Label for the Z axis.

  • azim (float) – Azimuthal angle of the 3D view (deg).

  • dist (float) – Camera distance parameter.

  • elev (float) – Elevation angle of the 3D view (deg).

  • factor (float) – Multiplier applied to Y for unit scaling.

Returns:

Displays a Matplotlib 3D surface plot.

Return type:

None

radicalpy.plot.plot_autocorrelation_fit(t_j, acf_j, acf_j_fit, zero_point_crossing_j)[source]

Autocorrelation with fit (log-scale).

Plots an autocorrelation function up to its first zero crossing together with a fitted curve provided in acf_j_fit["fit"].

Parameters:

  • t_j (np.ndarray) – Lag times τ (s).

  • acf_j (np.ndarray) – Autocorrelation values g(τ).

  • acf_j_fit (Mapping[str, np.ndarray]) – Fit results containing key “fit”.

  • zero_point_crossing_j (int) – Index of first zero crossing for g(τ).

Returns:

Displays a Matplotlib line plot on a log-x axis.

Return type:

None

radicalpy.plot.plot_bhalf_time(ts, bhalf_time, fit_error_time, style='ro', factor=1000000.0)[source]

B_{1/2}-field width vs time.

Plots discrete estimates of the half-field width B_{1/2} against time with error bars from fitting uncertainties.

Parameters:

  • ts (np.ndarray) – Time points (s).

  • bhalf_time (np.ndarray) – Estimated B_{1/2} (mT) per time point.

  • fit_error_time (np.ndarray) – Fit errors (e.g., covariance-derived) where index [1, i] is used as y-error.

  • style (str) – Matplotlib line/marker style (default: “ro”).

  • factor (float) – Multiplier applied to x-values for unit scaling.

Returns:

Displays a scatter plot with error bars.

Return type:

None

radicalpy.plot.plot_exchange_interaction_in_solution(ts, trajectory_data, j)[source]

Exchange interaction and separation vs time.

Plots radical-pair separation and corresponding exchange interaction on twin y-axes against time.

Parameters:

  • ts (np.ndarray) – Time points (s or ns as labeled).

  • trajectory_data (np.ndarray) – Trajectory array where column 1 contains separations (m).

  • j (np.ndarray) – Exchange interaction values (mT), same length as ts.

Returns:

Displays a Matplotlib plot with twin y-axes.

Return type:

None

radicalpy.plot.plot_general(xdata, ydata, xlabel, ylabel, style='-', label=[], colors='r', factor=1)[source]

General line plot.

Convenience wrapper for a simple x–y line plot with label, style, color, and optional scaling of the x-values.

Parameters:

  • xdata (np.ndarray) – X-axis data.

  • ydata (np.ndarray) – Y-axis data.

  • xlabel (str) – Label for the X axis.

  • ylabel (str) – Label for the Y axis.

  • style (str) – Matplotlib style string (default: “-“).

  • label (list or str) – Legend label(s).

  • colors (Any) – Matplotlib color spec (default: “r”).

  • factor (float) – Multiplier applied to x-values.

Returns:

Displays a Matplotlib line plot.

Return type:

None

radicalpy.plot.plot_molecule(ax, labels, elements, coords, bonds, show_labels=True, show_atoms=False)[source]

Molecular stick model.

Draws a simple 3D molecular diagram from atom coordinates and bond pairs, optionally showing atom spheres and/or atom labels.

Parameters:

  • ax (mpl_toolkits.mplot3d.Axes3D) – Target 3D axes.

  • labels (Sequence[str]) – Per-atom labels (e.g., element symbols).

  • elements (Sequence[str]) – Per-atom element identifiers.

  • coords (np.ndarray) – Cartesian coordinates (Å or nm) of shape (N, 3).

  • bonds (Sequence[Tuple[int,int]]) – Index pairs of bonded atoms.

  • show_labels (bool) – Whether to render atom labels.

  • show_atoms (bool) – Whether to render atom scatter points.

Returns:

Draws onto the provided axes.

Return type:

None

radicalpy.plot.plot_sphere(ax, radius=1, color='black', alpha=0.05)[source]

Plot a wire-framed, lightly shaded sphere on a Matplotlib 3D axis.

The surface is rendered via Axes3D.plot_surface using a spherical parameterisation, and three faint coordinate axes (x, y, z) are overlaid through the sphere’s center for orientation. A thin great-circle outline is also drawn.

Parameters:

  • ax – matplotlib.axes._subplots.Axes3DSubplot or mpl_toolkits.mplot3d.Axes3D A 3D Matplotlib axes object on which to draw.

  • radius – float, optional Sphere radius. Default is 1.

  • color – str or tuple, optional Base color for the sphere surface and outlines. Default is 'black'.

  • alpha – float, optional Surface transparency in [0, 1]. Default is 0.05.

Returns:

None

The function draws on ax in place and returns None.

radicalpy.plot.set_equal_aspect(ax, X, Y, Z)[source]

Equalise 3D aspect by data range.

Sets axis limits so the plotted 3D data have equal scale along each axis.

Parameters:

  • ax (mpl_toolkits.mplot3d.Axes3D) – Target 3D axes.

  • X (np.ndarray) – X coordinates.

  • Y (np.ndarray) – Y coordinates.

  • Z (np.ndarray) – Z coordinates.

Returns:

Adjusts axis limits in place.

Return type:

None

radicalpy.plot.spin_state_labels(sim)[source]

Spin-state labels for a radical pair.

Generates human-readable kets for singlet/triplet sublevels combined with radical spin projections, suitable for plot tick labels.

Parameters:

sim (HilbertSimulation) – Simulation containing exactly two radicals.

Returns:

LaTeX-formatted labels like "$\vert T_+, +1/2 \rangle$".

Return type:

list[str]

radicalpy.plot.visualise_tensor(ax, tensor, rot_matrix, coords, colour)[source]

Rank-2 tensor visualisation on a sphere.

Projects a symmetric rank-2 tensor onto directions on the unit sphere, scales radius by the directional quadratic form, and plots the resulting surface in 3D after applying a rotation and translation.

Parameters:

  • ax (mpl_toolkits.mplot3d.Axes3D) – Target 3D axes.

  • tensor (np.ndarray) – 3×3 tensor (e.g., hyperfine or g-tensor).

  • rot_matrix (np.ndarray) – 3×3 rotation matrix to orient the tensor.

  • coords (np.ndarray) – 3-vector translation applied to the surface.

  • colour (Any) – Matplotlib color spec for the surface.

Returns:

The rendered surface coordinates with shape (Nθ, Nφ, 3).

Return type:

np.ndarray

radicalpy.relaxation module

Relaxation superoperators in Liouville space for radical-pair dynamics.

This module defines a family of incoherent (relaxational) processes as Liouville-space superoperators that act on vectorised density matrices. Each process contributes an additive term subH (a superoperator) that can be combined with coherent Liouvillians to evolve spin systems subject to stochastic environmental interactions.

The implementations follow standard models in spin chemistry and magnetic resonance, including singlet–triplet dephasing, random local fields, triplet–triplet relaxation/dephasing, T₁/T₂ Bloch-type relaxation, and anisotropic g-tensor modulation under rotational diffusion.

Classes:
Helper functions:

  • _g_tensor_anisotropy_term: Per-radical superoperator contribution for g-anisotropy with Kivelson spectral densities.

Key concepts:

  • Liouville space / superoperators: Operators acting on vectorized density matrices. Additive relaxation terms are represented by subH.

  • Projection operators: Q_S, Q_T, Q_{T+}, Q_{T0}, Q_{T-} select singlet/triplet subspaces and enter several models.

  • Spin operators: Single-spin S_x, S_y, S_z per radical index, used to build isotropic/anisotropic relaxation channels via Kronecker products.

  • Spectral density: Frequency- and correlation-time–dependent scaling (e.g., Kivelson form) for stochastic modulation processes.

Usage pattern:

  1. Instantiate a relaxation class with its parameters (e.g., rate_constant, tau_c, g principal values, omega).

  2. Call init(sim) with a LiouvilleSimulation providing projection_operator(...) and spin_operator(radical_idx, axis).

  3. Retrieve the superoperator via the instance’s subH attribute and add it to your total Liouvillian.

Args conventions (per-class):

  • rate_constant (float, s⁻¹): Overall relaxation rate used in most models.

  • g1, g2 (list[float]): Principal components of the radicals’ g-tensors.

  • omega1, omega2 (float, rad·s⁻¹·mT⁻¹ × B): Larmor prefactors; the actual frequency is typically omega * B0. See note below.

  • tau_c1, tau_c2 (float, s): Rotational correlation times.

Notes

  • Magnetic-field dependence: For GTensorAnisotropy, the instantaneous Larmor frequency is proportional to the applied field B0. If B0 varies (e.g., in MARY scans), compute omega dynamically from B0 rather than treating it as a fixed constant.

  • Vectorisation size: Superoperators are built with Kronecker products sized to the underlying Hilbert spaces of the radicals in sim.

  • Normalisation and units: Rate-like parameters are in s⁻¹ unless noted. Axis labels/units are not handled here (these are dynamical, not plotting utilities).

References

Requirements:

  • numpy for array/Kronecker algebra.

  • A simulation object exposing: projection_operator(State.*) and spin_operator(idx, 'x'|'y'|'z'), as used in the init methods.

raises ValueError:

May be raised upstream if the simulation object does not support the required operators or radical indexing.

See also

radicalpy.simulation.LiouvilleSimulation, radicalpy.states.State, and the corresponding coherent Liouvillian construction in your codebase.

class radicalpy.relaxation.DipolarModulation(rate_constant)[source]

Bases: LiouvilleRelaxationBase

Dipolar modulation relaxation superoperator.

Source: Kattnig et al. New J. Phys., 18, 063007 (2016).

>>> DipolarModulation(rate_constant=1e6)
Relaxation: DipolarModulation
Rate constant: 1000000.0
init(sim)[source]

See radicalpy.simulation.HilbertIncoherentProcessBase.init.

class radicalpy.relaxation.GTensorAnisotropy(g1, g2, omega1, omega2, tau_c1, tau_c2)[source]

Bases: LiouvilleRelaxationBase

g-tensor anisotropy relaxation superoperator.

Source: Kivelson, J. Chem. Phys. 33, 1094 (1960).

Parameters:

  • g1 (list) – The principle components of g-tensor of the first radical.

  • g2 (list) – The principle components of g-tensor of the second radical.

  • omega1 (float) – The Larmor frequency of the first radical (rad/s/mT).

  • omega2 (float) – The Larmor frequency of the second radical (rad/s/mT).

  • tau_c1 (float) – The rotational correlation time of the first radical (s).

  • tau_c2 (float) – The rotational correlation time of the second radical (s).

>>> GTensorAnisotropy(g1=[2.0032, 1.9975, 2.0014],
...                   g2=[2.00429, 2.00389, 2.00216],
...                   omega1=-158477366720.7,
...                   omega2=-158477366720.7,
...                   tau_c1=5e-12,
...                   tau_c2=100e-12)
Relaxation: GTensorAnisotropy
g1: [2.0032, 1.9975, 2.0014]
g2: [2.00429, 2.00389, 2.00216]
omega1: -158477366720.7
omega2: -158477366720.7
tau_c1: 5e-12
tau_c2: 1e-10
init(sim)[source]

See radicalpy.simulation.HilbertIncoherentProcessBase.init.

class radicalpy.relaxation.LiouvilleRelaxationBase(rate_constant)[source]

Bases: LiouvilleIncoherentProcessBase

Base class for relaxation superoperators (Liouville space).

class radicalpy.relaxation.RandomFields(rate_constant)[source]

Bases: LiouvilleRelaxationBase

Random fields relaxation superoperator.

Source: Kattnig et al. New J. Phys., 18, 063007 (2016).

>>> RandomFields(rate_constant=1e6)
Relaxation: RandomFields
Rate constant: 1000000.0
init(sim)[source]

See radicalpy.simulation.HilbertIncoherentProcessBase.init.

class radicalpy.relaxation.SingletTripletDephasing(rate_constant)[source]

Bases: LiouvilleRelaxationBase

Singlet-triplet dephasing relaxation superoperator.

Source: Shushin, Chem. Phys. Lett. 181, 2,3, 274-278 (1991).

>>> SingletTripletDephasing(rate_constant=1e6)
Relaxation: SingletTripletDephasing
Rate constant: 1000000.0
init(sim)[source]

See radicalpy.simulation.HilbertIncoherentProcessBase.init.

class radicalpy.relaxation.T1Relaxation(rate_constant)[source]

Bases: LiouvilleRelaxationBase

T1 (spin-lattice, longitudinal, thermal) relaxation superoperator.

Source: Bloch, Phys. Rev. 70, 460-474 (1946).

>>> T1Relaxation(rate_constant=1e6)
Relaxation: T1Relaxation
Rate constant: 1000000.0
init(sim)[source]

See radicalpy.simulation.HilbertIncoherentProcessBase.init.

class radicalpy.relaxation.T2Relaxation(rate_constant)[source]

Bases: LiouvilleRelaxationBase

T2 (spin-spin, transverse) relaxation superoperator.

Source: Bloch, Phys. Rev. 70, 460-474 (1946).

>>> T2Relaxation(rate_constant=1e6)
Relaxation: T2Relaxation
Rate constant: 1000000.0
init(sim)[source]

See radicalpy.simulation.HilbertIncoherentProcessBase.init.

class radicalpy.relaxation.TripletTripletDephasing(rate_constant)[source]

Bases: LiouvilleRelaxationBase

Triplet-triplet dephasing relaxation superoperator.

Source: Gorelik et al. J. Phys. Chem. A 105, 8011-8017 (2001).

>>> TripletTripletDephasing(rate_constant=1e6)
Relaxation: TripletTripletDephasing
Rate constant: 1000000.0
init(sim)[source]

See radicalpy.simulation.HilbertIncoherentProcessBase.init.

class radicalpy.relaxation.TripletTripletRelaxation(rate_constant)[source]

Bases: LiouvilleRelaxationBase

Triplet-triplet relaxation superoperator.

Source: Miura et al. J. Phys. Chem. A 119, 5534-5544 (2015).

>>> TripletTripletRelaxation(rate_constant=1e6)
Relaxation: TripletTripletRelaxation
Rate constant: 1000000.0
init(sim)[source]

See radicalpy.simulation.HilbertIncoherentProcessBase.init.

radicalpy.shared module

Shared utilities for physical constants and data paths.

This module provides a lightweight mechanism to load numerically usable physical constants from a JSON file while preserving rich metadata (e.g., units, symbols, references):

  • Constant: a subclass of float that carries a details attribute (a types.SimpleNamespace) with auxiliary information about the constant. You can use a Constant anywhere a float is expected, and still access metadata via .details, e.g. constants.mu_B.details.units.

  • Constant.fromjson(path): load a JSON mapping of name → {value, …} into a SimpleNamespace of Constant objects, accessible by attribute.

  • DATA_DIR: path to bundled data files.

  • constants: the default namespace of constants loaded from DATA_DIR/constants.json.

class radicalpy.shared.Constant(details: dict)[source]

Bases: float

Constan class.

Extends float with the Constant.details member.

details: SimpleNamespace

Details (e.g. units) of the constant.

static fromjson(json_file)[source]

Read all constants from the JSON file.

Parameters:

json_file (str)

Returns:

A namespace containing all constants.

Return type:

SimpleNamespace

radicalpy.simulation module

Simulation utilities for spin dynamics.

This module provides classes and helpers to build quantum mechanical models in either Hilbert space (density matrices as square arrays) or Liouville space (vectorised densities / superoperators). It focuses on electron–nuclear spin systems typical of radical pairs and triplet pairs and supports common interactions and observables used in magnetic resonance and spin chemistry.

Main classes:

  • State :

    Enumerates common initial/observable spin states (singlet/triplet manifold, EPR observable, thermal equilibrium, etc.).

  • Basis :

    Choice of electron-pair basis: Zeeman or singlet–triplet (S/T).

  • HilbertSimulation :

    Core simulator that assembles Hamiltonians (Zeeman, hyperfine, exchange, dipolar, optional zero-field splitting), prepares initial density matrices, propagates them unitarily, and evaluates product probabilities/yields.

  • LiouvilleSimulation :

    Extends HilbertSimulation with Liouville-space evolution (superoperators, vectorised densities) for convenient inclusion of incoherent processes.

  • HilbertIncoherentProcessBase / LiouvilleIncoherentProcessBase :

    Base hooks to augment Hamiltonians or measured probabilities with phenomenological kinetics / relaxation.

  • SparseCholeskyHilbertSimulation :

    Hilbert-space variant optimised for large systems via sparse algebra and a Cholesky-factor time-stepping scheme.

  • SemiclassicalSimulation :

    Generates random semiclassical Hamiltonians for ensemble treatments.

Key features:

  • Spin operators for arbitrary particles and bases.

  • Zeeman (1D/3D), hyperfine (isotropic/tensor), exchange (J), dipolar (1D/3D),

and ZFS Hamiltonians. - Initial states: projection-based (S, T, etc.) and thermal equilibrium. - Time evolution: Hilbert (ρ → U ρ U†) and Liouville (ρ → e^{Lt} ρ). - Observables: projection operators and product probabilities/yields. - Shape conventions:

  • Hilbert density: (dim, dim)

  • Liouville density (vectorised): (dim**2, 1)

Units & conventions:

  • Magnetic fields in mT; gyromagnetic ratios provided as gamma_mT.

  • Tensors follow x/y/z Cartesian ordering.

  • S/T transform from the Zeeman basis.

See also

  • utils.spherical_to_cartesian for field orientation.

  • Module docstrings of related packages/classes for data structures (Molecule,

nuclei, radicals, hyperfine data).

class radicalpy.simulation.Basis(value)[source]

Bases: Enum

Spin basis choices for electron pair subspace (Zeeman vs. singlet–triplet).

ST = 'ST'
ZEEMAN = 'Zeeman'
class radicalpy.simulation.HilbertIncoherentProcessBase(rate_constant)[source]

Bases: object

adjust_hamiltonian(*args, **kwargs)[source]

Optionally modify the Hamiltonian in-place.

Subclasses can implement e.g. Redfield/Lindblad additions. The default implementation is a no-op and returns None.

adjust_product_probabilities(*args, **kwargs)[source]

Optionally post-process simulated product probabilities.

Can be used to fold in kinetic schemes in Hilbert simulations. The default implementation is a no-op and returns None.

init(sim)[source]

Hook to (optionally) initialise process state with the simulation.

Subclasses may cache shapes, indices, or precomputed operators derived from sim. The default implementation does nothing.

property rate_constant: float

Rate of the incoherent process.

class radicalpy.simulation.HilbertSimulation(molecules, custom_gfactors=False, basis=Basis.ST)[source]

Bases: object

Quantum simulation base class.

Parameters:

  • molecules (list[Molecule]) – List of Molecule objects.

  • custom_gfactor (bool) – Flag to use g-factors instead of the default gyromagnetic ratio gamma.

>>> HilbertSimulation([Molecule.fromdb("flavin_anion", ["N5"]),
...                    Molecule.fromdb("tryptophan_cation", ["Hbeta1", "H1"])])
Number of electrons: 2
Number of nuclei: 3
Number of particles: 5
Multiplicities: [2, 2, 3, 2, 2]
Magnetogyric ratios (mT): [-176085963.023, -176085963.023, 19337.792, 267522.18744, 267522.18744]
Nuclei: [14N(19337792.0, 3, 0.5141 <anisotropic available>), 1H(267522187.44, 2, 1.605 <anisotropic available>), 1H(267522187.44, 2, -0.5983 <anisotropic available>)]
Couplings: [0, 1, 1]
HFCs (mT): [0.5141 <anisotropic available>, 1.605 <anisotropic available>, -0.5983 <anisotropic available>]
ST_basis(M)[source]

Transform an operator from the Zeeman basis to the S/T basis.

Parameters:

M – Operator in the Zeeman basis.

Returns:

The operator expressed in the S/T basis.

static apply_hilbert_kinetics(time, product_probabilities, kinetics)[source]

Apply kinetic post-processing to product probabilities in Hilbert space.

Each object in kinetics may rescale or filter the probabilities as a function of time (e.g., recombination channels).

apply_liouville_hamiltonian_modifiers(H, modifiers)[source]

Apply (Hilbert) incoherent process modifiers to the Liouville Hamiltonian.

Each modifier is initialised with self once, then given the chance to adjust H in place. This is a no-op for Hilbert simulations.

bloch_redfield_liouvillian(H, bath, noise, secular=True)[source]

Bloch-Redfield tensor builder.

Parameters:

  • H (np.ndarray) – Hamiltonian (time-independent).

  • bath (list[np.ndarray]) – System operators coupling to the bath.

  • noise (list[np.ndarray]) – spectra[k](w) gives noise power at frequency w, for bath[k].

  • secular (bool) – If True, use the dw_min/10 frequency matching.

Returns:

Bloch-Redfield Liouvillian in the energy basis. eigvecs (np.ndarray): Columns are eigenkets of H.

Return type:

L (np.ndarray)

bloch_redfield_solver(L_e, eigvecs, rho0, time, obs=None, rtol=1e-07, atol=1e-07)[source]

Evolve Bloch-Redfield master equation in the energy basis.

Parameters:

  • L_e (np.ndarray) – Liouvillian in energy basis.

  • eigvecs (np.ndarray) – Energy eigenvectors.

  • rho0 (np.ndarray) – Initial density matrix in lab basis.

  • time (np.ndarray) – Time in seconds.

  • obs (list[np.ndarray]) – Expectation operator in lab basis.

Returns:

{‘states’: […]} or {‘expect’: [arrays…]}

Return type:

results (dict)

bloch_redfield_time_evolution(H, rho0, time, bath, noise=None, obs=None, secular=True, **kwargs)[source]

Bloch-Redfield Liouvillian time evolution solver.

static convert(H)[source]

Convert a Hilbert-space Hamiltonian to the simulator’s working space.

In Hilbert simulations this is the identity; Liouville overrides it.

Return type:

ndarray

property coupling

List mapping each nucleus to the index of its parent radical.

dipolar_hamiltonian(D)[source]

Construct the Dipolar Hamiltonian.

Construct the Dipolar Hamiltonian based on dipolar coupling constant or dipolar interaction tensor D between two electrons. Depending on the type of D, the 1D or the 3D version is invoked.

See:

Parameters:

D (float | np.ndarray) – dipolar coupling constant or dipolar interaction tensor.

Returns:

The Dipolar Hamiltonian corresponding to the system described by the HilbertSimulation object and dipolar coupling constant or dipolar interaction tensor D.

Return type:

np.ndarray

dipolar_hamiltonian_1d(D)[source]

Construct the 1D Dipolar Hamiltonian.

Construct the Dipolar Hamiltonian based on dipolar coupling constant D between two electrons.

The dipolar coupling constant can be obtained using radicalpy.estimations.dipolar_interaction_isotropic.

Parameters:

D (float) – dipolar coupling constant in mT.

Returns:

The 1D Dipolar Hamiltonian corresponding to the system described by the HilbertSimulation object and dipolar coupling constant D.

Return type:

np.ndarray

dipolar_hamiltonian_3d(dipolar_tensor)[source]

Construct the 3D Dipolar Hamiltonian.

Construct the Dipolar Hamiltonian based on dipolar interaction tensor D between two electrons.

The dipolar coupling tensor can be obtained using radicalpy.estimations.dipolar_interaction_anisotropic.

Parameters:

dipolar_tensor (np.ndarray) – dipolar interaction tensor in mT.

Returns:

The 3D Dipolar Hamiltonian corresponding to the system described by the HilbertSimulation object and dipolar interaction tensor D.

Return type:

np.ndarray

exchange_hamiltonian(J, prod_coeff=2)[source]

Construct the exchange Hamiltonian.

Construct the exchange (J-coupling) Hamiltonian based on the coupling constant J between two electrons.

The J-coupling constant can be obtained using:

Parameters:

  • 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:

The exchange (J-coupling) Hamiltonian corresponding to the system described by the HilbertSimulation object and the coupling constant J.

Return type:

np.ndarray

get_eye(shape)[source]

Return an identity matrix of the requested dimension (dense).

Return type:

ndarray

property hamiltonian_size

Dimension of the Hilbert space (product of particle multiplicities).

hyperfine_hamiltonian(hfc_anisotropy=False)[source]

Construct the Hyperfine Hamiltonian.

Construct the Hyperfine Hamiltonian. If hfc_anisotropy is False, then the isotropic hyperfine coupling constants are used. If hfc_anisotropy is True then the full hyperfine tensors are used (assuming they are available for all the nuclei of the molecule in the database, otherwise an exception is raised).

Parameters:

hfc_anisotropy (bool) – Use isotropic hyperfine coupling constants if False, use full hyperfine tensors if False.

Returns:

The Hyperfine Hamiltonian corresponding to the system described by the HilbertSimulation object.

Return type:

np.ndarray

initial_density_matrix(state, H)[source]

Construct the initial density matrix for time evolution.

Builds a properly normalised density operator in Hilbert space corresponding to the chosen initial state.

  • For most states (e.g. singlet, triplet), this is simply the normalised projection operator.

  • For the special case State.EQUILIBRIUM, the density is taken as the thermal equilibrium state, approximated via \(\rho_0 \propto e^{-i H \beta}\) where \(\beta = \hbar / (k_B T)\) is provided via projection_operator.

Parameters:

  • state (State) – Target spin state (e.g., SINGLET, TRIPLET, EQUILIBRIUM).

  • H (np.ndarray) – Spin Hamiltonian in Hilbert space.

Returns:

Normalised density matrix in Hilbert space.

Return type:

np.ndarray

lindblad_liouvillian(H, Ls=[])[source]

Assemble the Liouville-space generator for coherent + Lindblad dynamics.

This builds the superoperator \(\mathcal{L} = -i\,[H,\cdot] + \sum_k \mathcal{D}[L_k]\) using the column–stacking (vec) convention. The coherent part is implemented as

\[-i[H,\rho] \;\mapsto\; -i\,(I\otimes H - H^{\mathsf{T}}\otimes I),\]

and each dissipator is

\[\mathcal{D}[L](\rho) = L\,\rho\,L^{\dagger} - \tfrac{1}{2}\{L^{\dagger}L,\rho\} \;\mapsto\; L^{\ast}\otimes L - \tfrac{1}{2}\Big( I\otimes L^{\dagger}L + (L^{\mathsf{T}}L^{\ast})\otimes I \Big).\]
Parameters:

  • H (ndarray) – Hilbert-space Hamiltonian. For Hermitian problems, a complex dtype (e.g. np.complex128) is recommended.

  • Ls – Iterable of collapse operators L_k (each N×N) defining the Lindblad dissipators. If you have physical rates γ_k, scale your operators as L_k ← √γ_k · C_k before passing them.

Returns:

The full Liouville-space generator \(\mathcal{L}\) suitable for acting on vec(ρ) with the column-stacking convention.

Notes:

  • The mapping uses vec(O ρ) = (I ⊗ O^{\mathsf{T}}) vec(ρ) and vec(ρ O) = (O^{\mathsf{T}} ⊗ I) vec(ρ).

  • H enters only through the coherent commutator; the dissipator depends solely on Ls.

  • All arrays are combined with NumPy kron; performance-critical paths may prefer sparse representations.

property nuclei

Flattened list of all nuclei across molecules.

observable_projection_operator(state)[source]

Projection operator used to compute observable expectation values.

Simply forwards to projection_operator() in Hilbert space.

Return type:

ndarray

property particles

Concatenated list of electron spins followed by nuclei.

product_operator(idx1, idx2, h=1.0)[source]

Construct the (1D) product operator.

Construct the 1D (isotropic) product operator of two particles in the spin system.

Parameters:

  • idx1 (int) – Index of the first particle.

  • idx2 (int) – Index of the second particle.

  • h (float) – Isotropic interaction constant.

Returns:

Product operator for particles corresponding to idx1 and idx2 with isotropic interaction constant h.

Return type:

np.ndarray

product_operator_3d(idx1, idx2, h)[source]

Construct the 3D product operator.

Construct the 3D (anisotropic) product operator of two particles in the spin system.

Parameters:

  • idx1 (int) – Index of the first particle.

  • idx2 (int) – Index of the second particle.

  • h (np.ndarray) – Anisotropic interaction tensor.

Returns:

Product operator for particles corresponding to idx1 and idx2 with anisotropic interaction tensor h.

Return type:

np.ndarray

product_probability(obs, rhos)[source]

Calculate the probability of the observable from the densities.

Return type:

ndarray

static product_yield(product_probability, time, k)[source]

Calculate the product yield and the product yield sum.

projection_operator(state, T=298, chi=1.5707963267948966)[source]

Construct the projection operator corresponding to a state.

Parameters:

  • state (State) – The target state which is projected out of the density matrix.

  • T (float) – Temperature for the EQUILIBRIUM projection operator (K).

  • chi (float) – Degree of polarisation for the CISS projection operator.

Returns:

Projection operator corresponding to the State state.

Return type:

np.ndarray

propagate(propagator, rho)[source]

Propagate the density matrix (Hilbert space).

Propagates the density matrix using the propagator obtained using the unitary_propagator method.

\[\rho (t) = \mathbf{U} \rho_0 \mathbf{U}^*\]

See also: unitary_propagator and time_evolution.

Parameters:

  • propagator (np.ndarray) – Unitary operator obtained via the unitary_propagator method.

  • rho (np.ndarray) – (Initial) density matrix.

Returns:

The new density matrix after the unitary operator was applied to it.

Return type:

np.ndarray

property radicals

List of electron spins (one per molecule).

spin_operator(idx, axis)[source]

Construct the spin operator.

Construct the spin operator for the particle with index idx in the HilbertSimulation.

Parameters:

  • idx (int) – Index of the particle.

  • axis (str) – Axis, i.e. "x", "y" or "z".

Returns:

Spin operator for a particle in the HilbertSimulation system with indexing idx and axis axis.

Return type:

np.ndarray

spin_orbit_hamiltonian(soc_vectors, electron_spins=None, sparse=False)[source]

Spin–orbit interaction Hamiltonian.

H_SO = sum_i (Lx_i * S_i,x + Ly_i * S_i,y + Lz_i * S_i,z)

Parameters:

  • soc_vectors – list[int], optional Which electrons (0-based) to include. Default: all.

  • sparse (bool) – Return sparse CSR if True, else dense ndarray.

Returns:

ndarray or csr_matrix

Return type:

H_SO

time_evolution(init_state, time, H)[source]

Evolve the system through time.

See Also: - HilbertSimulation.unitary_propagator - HilbertSimulation.propagate - LiouvilleSimulation.unitary_propagator - LiouvilleSimulation.propagate

Parameters:

  • init_state (State) – Initial State of the density matrix (see projection_operator).

  • time (np.ndarray) – An sequence of (uniform) time points, usually created using np.arange or np.linspace.

  • H (np.ndarray) – Hamiltonian operator.

Returns:

Return a sequence of density matrices evolved through time, starting from a given initial state using the Hamiltonian H.

Return type:

np.ndarray

Examples

>>> molecules = [Molecule.fromdb("flavin_anion", ["N5"]),
...              Molecule("Z")]
>>> sim = HilbertSimulation(molecules)
>>> H = sim.total_hamiltonian(B0=0, J=0, D=0)
>>> time = np.arange(0, 2e-6, 5e-9)
>>> time.shape
(400,)
>>> rhos = sim.time_evolution(State.SINGLET, time, H)
>>> rhos.shape
(400, 12, 12)
total_hamiltonian(B0, J, D, theta=None, phi=None, hfc_anisotropy=False)[source]

Construct the total Hamiltonian.

The total Hamiltonian is the sum of Zeeman, Hyperfine, Exchange and Dipolar Hamiltonian.

Parameters:
Returns:

The total Hamiltonian.

Return type:

np.ndarray

tp_quintet_minus_one_projop(SAx, SAy, SAz, SBx, SBy, SBz)[source]

Projection operator onto the triplet-pair quintet minus one (TP-Q-1) subspace.

Builds the projector from electron spin operators of radicals A and B.

tp_quintet_minus_two_projop(SAx, SAy, SAz, SBx, SBy, SBz)[source]

Projection operator onto the triplet-pair quintet minus two (TP-Q-2) subspace.

Builds the projector from electron spin operators of radicals A and B.

tp_quintet_plus_one_projop(SAx, SAy, SAz, SBx, SBy, SBz)[source]

Projection operator onto the triplet-pair quintet plus one (TP-Q+1) subspace.

Builds the projector from electron spin operators of radicals A and B.

tp_quintet_plus_two_projop(SAx, SAy, SAz, SBx, SBy, SBz)[source]

Projection operator onto the triplet-pair quintet plus two (TP-Q+2) subspace.

Builds the projector from electron spin operators of radicals A and B.

tp_quintet_projop(SAx, SAy, SAz, SBx, SBy, SBz)[source]

Projection operator onto the triplet-pair quintet (TP-Q) subspace.

Builds the projector from electron spin operators of radicals A and B.

tp_quintet_zero_projop(SAx, SAy, SAz, SBx, SBy, SBz)[source]

Projection operator onto the triplet-pair quintet zero (TP-Q0) subspace.

Builds the projector from electron spin operators of radicals A and B.

tp_singlet_projop(SAx, SAy, SAz, SBx, SBy, SBz)[source]

Projection operator onto the triplet-pair singlet (TP-S) subspace.

Builds the projector from electron spin operators of radicals A and B.

tp_triplet_minus_projop(SAx, SAy, SAz, SBx, SBy, SBz)[source]

Projection operator onto the triplet-pair triplet minus (TP-T-) subspace.

Builds the projector from electron spin operators of radicals A and B.

tp_triplet_plus_projop(SAx, SAy, SAz, SBx, SBy, SBz)[source]

Projection operator onto the triplet-pair triplet plus (TP-T+) subspace.

Builds the projector from electron spin operators of radicals A and B.

tp_triplet_projop(SAx, SAy, SAz, SBx, SBy, SBz)[source]

Projection operator onto the triplet-pair triplet (TP-T) subspace.

Builds the projector from electron spin operators of radicals A and B.

tp_triplet_zero_projop(SAx, SAy, SAz, SBx, SBy, SBz)[source]

Projection operator onto the triplet-pair triplet zero (TP-T0) subspace.

Builds the projector from electron spin operators of radicals A and B.

static unitary_propagator(H, dt)[source]

Create unitary propagator (Hilbert space).

Create unitary propagator matrices U and U* for time evolution of the density matrix in Hilbert space (for the spin Hamiltonian H).

\[\begin{split}\mathbf{U} =& \exp( -i \hat{H} t ) \\ \mathbf{U}^* =& \exp( +i \hat{H} t )\end{split}\]

See also: propagate and time_evolution.

Parameters:

  • H (np.ndarray) – Spin Hamiltonian in Hilbert space.

  • dt (float) – Time evolution timestep.

Returns:

Two matrices (as tensor) in Hilbert space which are used by the propagate method to perform a single step in the time_evolution method.

Return type:

np.ndarray

Examples

>>> molecules = [Molecule.fromdb("flavin_anion", ["N5"]),
...              Molecule("Z")]
>>> sim = HilbertSimulation(molecules)
>>> H = sim.total_hamiltonian(B0=0, J=0, D=0)
>>> Up, Um = sim.unitary_propagator(H, 3e-9)
>>> Up.shape, Um.shape
((12, 12), (12, 12))
zeeman_hamiltonian(B0, B_axis='z', theta=None, phi=None)[source]

Construct the Zeeman Hamiltonian (1D or 3D).

Construct the Zeeman Hamiltonian based on the external magnetic field B0. If the angles theta and phi are also provided, the 3D Zeeman Hamiltonian is constructed (by invoking the zeeman_hamiltonian_3d), otherwise the 1D Zeeman Hamiltonian is constructed (by invoking the zeeman_hamiltonian_1d).

Parameters:
Returns:

The Zeeman Hamiltonian corresponding to the system described by the HilbertSimulation object and the external magnetic field intensity B0 and angles theta and phi.

Return type:

np.ndarray

zeeman_hamiltonian_1d(B0, axis)[source]

Construct the 1D Zeeman Hamiltonian.

Construct the 1D Zeeman Hamiltonian based on the external magnetic field B0.

Parameters:

  • B0 (float) – External magnetic field intensity (milli Tesla).

  • axis (str) – Axis of the magnetic field.

Returns:

The Zeeman Hamiltonian corresponding to the system described by the HilbertSimulation object and the external magnetic field intensity B0.

Return type:

np.ndarray

zeeman_hamiltonian_3d(B0, theta=0, phi=0)[source]

Construct the 3D Zeeman Hamiltonian.

Construct the 3D Zeeman Hamiltonian based on the external magnetic field B0 and angles theta and phi.

Parameters:

  • B0 (float) – External magnetic field intensity (milli Tesla).

  • theta (Optional[float]) – rotation (polar) angle between the external magnetic field and the fixed molecule.

  • phi (Optional[float]) – rotation (azimuth) angle between the external magnetic field and the fixed molecule.

Returns:

The Zeeman Hamiltonian corresponding to the system described by the HilbertSimulation object and the external magnetic field intensity B0 and angles theta and phi.

Return type:

np.ndarray

zero_field_splitting_hamiltonian(D, E)[source]

Build the zero-field splitting (ZFS) Hamiltonian.

Constructs the second-rank ZFS contribution for (typically triplet, S=1) spins using the conventional axial/rhombic parameters D and E. For each particle index idx, the operator \(D (S_z^2 - \tfrac{1}{3} \mathbf{S}^2) + E (S_x^2 - S_y^2)\) is formed and summed. The parameters are scaled into frequency units via the first radical’s electron gyromagnetic ratio.

Parameters:

  • D – Axial ZFS parameter (in the same field/frequency units used elsewhere; internally scaled by -radicals[0].gamma_mT).

  • E – Rhombic ZFS parameter (scaled identically to D).

Return type:

ndarray

Returns:

The ZFS Hamiltonian matrix in the current simulation basis.

Notes

  • This implementation multiplies D and E by -gamma_mT of the first radical to convert to angular frequency units consistent with the rest of the Hamiltonian.

  • The loop applies the operator to every particle index; in practice, ZFS is only meaningful for S≥1 spins (e.g., an excited triplet). Ensure your particle list reflects that physical situation.

class radicalpy.simulation.LiouvilleIncoherentProcessBase(rate_constant)[source]

Bases: HilbertIncoherentProcessBase

adjust_hamiltonian(H)[source]

Subtract the prebuilt incoherent sub-Hamiltonian from H in Liouville space.

Expects subclasses to define self.subH with the correct shape.

class radicalpy.simulation.LiouvilleSimulation(molecules, custom_gfactors=False, basis=Basis.ST)[source]

Bases: HilbertSimulation

static convert(H)[source]

Convert the Hamiltonian from Hilbert to Liouville space.

Return type:

ndarray

property hamiltonian_size

Dimension of the Liouville space (Hilbert dimension squared).

initial_density_matrix(state, H)[source]

Create an initial density matrix for time evolution of the spin Hamiltonian density matrix.

Parameters:

  • state (State) – a string = spin state projection operator

  • spins – an integer = sum of the number of electrons and nuclei

  • H (ndarray) – a matrix = spin Hamiltonian in Hilbert space

Return type:

ndarray

Returns:

A matrix in Liouville space

liouville_projection_operator(state)[source]

Column-vectorised projector for Liouville space.

Flattens the Hilbert-space projector associated with state to shape (N², 1) for use with Liouville evolution.

Return type:

ndarray

observable_projection_operator(state)[source]

Row-vectorised observable (projectorᵀ) used for expectation values.

Return type:

ndarray

propagate(propagator, rho)[source]

Propagate the density matrix (Liouville space).

Propagates the density matrix using the propagator obtained using the unitary_propagator method.

\[\rho (t) = \mathbf{U} \rho_0\]

See also: unitary_propagator and time_evolution.

Parameters:

  • propagator (np.ndarray) – Unitary operator obtained via the unitary_propagator method.

  • rho (np.ndarray) – (Initial) density matrix.

Returns:

The new density matrix after the unitary operator was applied to it.

Return type:

np.ndarray

static unitary_propagator(H, dt)[source]

Create unitary propagator (Liouville space).

Create unitary propagator matrix U for the time evolution of the density matrix in Liouville space (for the spin Hamiltonian H).

\[\mathbf{U} = \exp( \hat{\hat{L}} t )\]

See also: propagate and time_evolution.

Parameters:

  • H (np.ndarray) – Spin Hamiltonian in Liouville space.

  • dt (float) – Time evolution timestep.

Returns:

The matrix in Liouville space which is used by the propagate method to perform a single step in the time_evolution method.

Return type:

np.ndarray

Examples

>>> molecules = [Molecule.fromdb("flavin_anion", ["N5"]),
...              Molecule("Z")]
>>> sim = LiouvilleSimulation(molecules)
>>> H = sim.total_hamiltonian(B0=0, J=0, D=0)
>>> sim.unitary_propagator(H, 3e-9).shape
(144, 144)
class radicalpy.simulation.SemiclassicalSimulation(molecules, custom_gfactors=False, basis=Basis.ST)[source]

Bases: LiouvilleSimulation

external_field_mT()[source]

Lab field B0 in mT as a 3-vector. Override if you have a field in the simulation.

Return type:

ndarray

property nuclei

Semiclassical model replaces nuclei by random static fields; return empty list.

semiclassical_HHs(num_samples, *, anisotropic=False)[source]

Generate semiclassical Hamiltonians.

For each radical r:

draw a static random hyperfine field B_hf,r (3D Gaussian), add the external field B0, then convert to ω_r = γ_r * (B0 + B_hf,r), and construct H = Σ_r ω_r · S_r. Repeat num_samples times.

Parameters:

  • num_samples (int) – Number of random Hamiltonian realisations.

  • anisotropic (bool, optional) – If True, draws from a 3D Gaussian

  • tensors. (with covariance Σ_r determined by anisotropic A)

  • False (If)

Returns:

Array (num_samples, D, D) of complex Hamiltonians.

Return type:

np.ndarray

class radicalpy.simulation.SparseCholeskyHilbertSimulation(molecules, custom_gfactors=False, basis=Basis.ST)[source]

Bases: HilbertSimulation

A simulation class that exploits

  • Sparsity of the Hamiltonian

  • Cholesky decomposition of the density matrix

Since density matrices are positive semi-definite, they can be decomposed into

rho = X X^T

where X is a lower triangular matrix by Cholesky decomposition.

In particular, when rho is a diagonal matrix, X = sqrt(rho).

This class accelerates the time evolution for the system with large Hilbert space (> 10^3).

ST_basis(M)[source]

Sparse S/T-basis transform of an operator.

Accepts dense or sparse input, converts to sparse as needed, and applies the electron-only change-of-basis while preserving sparsity.

Return type:

sparray

get_eye(shape)[source]

Return a sparse identity matrix of the requested dimension.

Return type:

sparray

product_probability(obs, rhos)[source]

Calculate the probability of the observable from the densities.

Return type:

ndarray

propagate(propagator, rho)[source]

Propagate a Cholesky factorisation (X, Xᵀ) one step.

Uses Um @ X to advance the factor; the updated pair (UmX, (UmX)ᴴ) represents the new density without forming ρ = X Xᵀ explicitly. Automatically switches to dense math when sparsity falls below a threshold (~30% zeros in factors).

Return type:

tuple[sparray | ndarray, sparray | ndarray]

spin_operator(idx, axis)[source]

Construct the spin operator.

Construct the spin operator for the particle with index idx in the HilbertSimulation.

Parameters:

  • idx (int) – Index of the particle.

  • axis (str) – Axis, i.e. "x", "y" or "z".

Returns:

Spin operator for a particle in the HilbertSimulation system with indexing idx and axis axis.

Return type:

np.ndarray

time_evolution(init_state, time, H)[source]

Evolve the system through time.

See Also: - HilbertSimulation.unitary_propagator - HilbertSimulation.propagate

Parameters:

  • init_state (State) – Initial State of the density matrix (see projection_operator).

  • time (np.ndarray) – An sequence of (uniform) time points, usually created using np.arange or np.linspace.

  • H (np.ndarray) – Hamiltonian operator.

Returns:

Return a sequence of density matrices (X, X^T) evolved through time, starting from a given initial state using the Hamiltonian H. Density matrices are obtained by X X^T.

Return type:

tuple[np.ndarray, np.ndarray]

Examples

>>> molecules = [Molecule.fromdb("flavin_anion", ["N5"]),
...              Molecule("Z")]
>>> sim = SparseCholeskyHilbertSimulation(molecules)
>>> H = sim.total_hamiltonian(B0=0, J=0, D=0)
>>> time = np.arange(0, 2e-6, 5e-9)
>>> time.shape
(400,)
>>> rhos = sim.time_evolution(State.SINGLET, time, H)
>>> len(rhos)
400
>>> rhos[0][0].shape
(12, 12)
>>> rhos[0][1].shape
(12, 12)
unitary_propagator(H, dt)[source]

Sparse/dense unitary for one time step in Hilbert space.

Builds U = exp(-i H dt) as a sparse matrix; if the result becomes too dense (fill > 50%), it is converted to a dense ndarray.

Parameters:

  • H (sparray) – Sparse Hamiltonian (CSC preferred).

  • dt (float) – Time step (s).

Return type:

sparray | ndarray

Returns:

Sparse or dense unitary matrix depending on fill ratio.

class radicalpy.simulation.State(value)[source]

Bases: Enum

Enumeration of common spin/population states used as observables or initial conditions.

CISS = 'CISS'
EPR = 'EPR'
EQUILIBRIUM = 'Eq'
SECOND_ORDER_POLARISATION = '2ndOrderPol'
SINGLET = 'S'
TP_QUINTET = 'TP_Q'
TP_QUINTET_MINUS_ONE = 'TP_Q-1'
TP_QUINTET_MINUS_TWO = 'TP_Q-2'
TP_QUINTET_PLUS_ONE = 'TP_Q+1'
TP_QUINTET_PLUS_TWO = 'TP_Q+2'
TP_QUINTET_ZERO = 'TP_Q0'
TP_SINGLET = 'TP_S'
TP_TRIPLET = 'TP_T'
TP_TRIPLET_MINUS = 'TP_T-'
TP_TRIPLET_PLUS = 'TP_T+'
TP_TRIPLET_ZERO = 'TP_T0'
TRIPLET = 'T'
TRIPLET_MINUS = 'T_-'
TRIPLET_PLUS = 'T_+'
TRIPLET_PLUS_MINUS = 'T_\\pm'
TRIPLET_X = 'X'
TRIPLET_Y = 'Y'
TRIPLET_Z = 'Z'
TRIPLET_ZERO = 'T_0'

radicalpy.tensornetwork module

Tensor network simulations for radical pair dynamics.

This module provides comprehensive tensor network-based simulations for radical pair dynamics, implementing three different approaches:

  1. Stochastic Matrix Product State (SMPS): Monte Carlo sampling approach in Hilbert space

  2. Locally Purified Matrix Product State (LPMPS): Purification-based deterministic approach in Hilbert space

  3. Vectorised Matrix Product Density Operator (VMPDO): Deterministic approach in Liouville space

The simulations support quantum mechanical treatments of nuclear spins, allowing for accurate modeling of: - Hyperfine interactions (isotropic and anisotropic) - Exchange coupling between radicals - Dipolar interactions - Zeeman effects - Haberkorn recombination terms

Key Features:

  • Support for both Hermitian and non-Hermitian Hamiltonians

  • Multiple computational backends (NumPy, JAX)

  • Parallel processing capabilities

  • GPU acceleration support (via JAX)

Example

Basic usage examples are available in:

  • examples/lpmps.py - Linear product MPS simulations

  • examples/smps.py - Stochastic MPS simulations

  • examples/vmpdo.py - Vectorised MPDO simulations

from radicalpy.data import Molecule
from radicalpy.simulation import State
import radicalpy.tensornetwork as tn

# Create radical pair molecules
molecules = [Molecule(...), Molecule(...)]

# Define system parameters
J, D, B0, kS, kT, anisotropy = ...

# Initialize simulation
sim = tn.LocallyPurifiedMPSSimulation(
    molecules=molecules,
    bond_dimension=16,
    jobname="lpmps",
)

# Construct tensor network Hamiltonian
H = sim.total_hamiltonian(
    B0=B0,
    J=J,
    D=D,
    kS=kS,
    kT=kT,
    hfc_anisotropy=anisotropy,
)

# Run simulation
partial_trace = sim.time_evolution(
    init_state=State.SINGLET,
    time=np.linspace(0, 1e-6, 1000),
    H=H
)

# Get population of singlet state
singlet_population = sim.product_probability(
    obs=State.SINGLET,
    rhos=partial_trace,
)
Dependencies:

Install the required tensor network simulation libraries:

pip install git+https://github.com/QCLovers/PyTDSCF

References

class radicalpy.tensornetwork.BaseMPSSimulation(molecules, custom_gfactors=False, basis=Basis.ST, bond_dimension=16, backend='numpy', integrator='arnoldi', jobname='tensornetwork')[source]

Bases: HilbertSimulation, ABC

Abstract base class for matrix product state simulations.

This class provides the foundation for tensor network-based simulations of radical pair dynamics. It implements common functionality for constructing Hamiltonians and managing tensor network operations.

The class supports three concrete implementations: - StochasticMPSSimulation: Monte Carlo sampling approach - LocallyPurifiedMPSSimulation: Purification-based method - MPDOSimulation: Direct density operator approach

Key Features:

  • Support for various Hamiltonian terms (Zeeman, hyperfine, exchange, dipolar)

  • Automatic matrix product operator (MPO) construction

  • Multiple computational backends (NumPy, JAX)

Parameters:

  • molecules (list[Molecule]) – List of molecules in the radical pair.

  • custom_gfactors (bool) – Whether to use custom g-factors. Default is False.

  • basis (Basis) – Basis set for the simulation. Only ST basis is supported.

  • bond_dimension (int) – Bond dimension for the matrix product state. Default is 16.

  • backend (Literal["numpy", "jax", "auto"]) – Computational backend. ‘auto’ selects ‘jax’ for bond_dimension > 32, ‘numpy’ otherwise.

  • integrator (Literal["arnoldi", "lanczos"]) – Integration method for time evolution. Use ‘lanczos’ for Hermitian Hamiltonians, ‘arnoldi’ for non-Hermitian.

  • jobname (str) – Job name for output files. Default is “tensornetwork”.

Raises:

Note

Important differences from standard HilbertSimulation:

  • Nuclear spins are automatically sorted by hyperfine coupling strength to reduce entanglement in the tensor network representation

  • Only ST basis is supported (Zeeman basis not available)

  • Each term of the Hamiltonian returns SumOfProducts objects instead of numpy arrays

  • The total Hamiltonian returns a MPO (list of numpy arrays) instead of a single numpy array

  • Time evolution returns diagonal elements of reduced density matrix (shape: nsteps, 4) or full reduced density matrices (shape: nsteps, 4, 4) instead of full density matrices (shape: nsteps, dim, dim)

  • Requires external dependencies (PyTDSCF, PyMPO) for tensor operations

ST_basis(M)[source]

Transform an operator from the Zeeman basis to the S/T basis.

Parameters:

M – Operator in the Zeeman basis.

Returns:

The operator expressed in the S/T basis.

apply_liouville_hamiltonian_modifiers(H, modifiers)[source]

apply_liouville_hamiltonian_modifiers is not used for BaseMPSSimulation. For Haberkorn term, include kS and kT for total_hamiltonian().

dipolar_hamiltonian(D)[source]

Construct the dipolar Hamiltonian as a SumOfProducts.

Parameters:

D (float | np.ndarray) – Dipolar coupling constant or tensor in mT.

Returns:

The dipolar Hamiltonian as a sum of products of operators.

Return type:

SumOfProducts

exchange_hamiltonian(J, prod_coeff=2)[source]

Construct the exchange Hamiltonian as a SumOfProducts.

Parameters:

  • J (float) – Exchange coupling constant in mT.

  • prod_coeff (float) – Coefficient for the S1·S2 term. Default is 2.

Returns:

The exchange Hamiltonian as a sum of products of operators.

Return type:

SumOfProducts

haberkorn_hamiltonian(kS, kT)[source]

Construct the Haberkorn recombination term as a SumOfProducts.

Parameters:

  • kS (float) – Singlet recombination rate in s^-1.

  • kT (float) – Triplet recombination rate in s^-1.

Returns:

The Haberkorn recombination Hamiltonian as a sum

of products of operators.

Return type:

SumOfProducts

hyperfine_hamiltonian(hfc_anisotropy=False)[source]

Construct the hyperfine Hamiltonian as a SumOfProducts.

Parameters:

hfc_anisotropy (bool) – Whether to include anisotropic hyperfine coupling. Default is False.

Returns:

The hyperfine Hamiltonian as a sum of products of operators.

Return type:

SumOfProducts

initial_density_matrix(state, H)[source]

Initial density matrix is not used for BaseMPSSimulation.

Return type:

ndarray

observable_projection_operator(state)[source]

Projection operator used to compute observable expectation values.

Simply forwards to projection_operator() in Hilbert space.

Return type:

ndarray

product_probability(obs, rhos)[source]

Calculate the probability of the observable from the densities.

Return type:

ndarray

propagate(propagator, rho)[source]

Propagate is not used for BaseMPSSimulation. Use time_evolution() instead.

Return type:

ndarray

spin_operator(idx, axis)[source]

Construct the spin operator.

Construct the spin operator for the particle with index idx in the HilbertSimulation.

Parameters:

  • idx (int) – Index of the particle.

  • axis (str) – Axis, i.e. "x", "y" or "z".

Returns:

Spin operator for a particle in the HilbertSimulation system with indexing idx and axis axis.

Return type:

np.ndarray

abstract time_evolution(init_state, time, H)[source]

Evolve the system through time.

See Also: - HilbertSimulation.unitary_propagator - HilbertSimulation.propagate - LiouvilleSimulation.unitary_propagator - LiouvilleSimulation.propagate

Parameters:

  • init_state (State) – Initial State of the density matrix (see projection_operator).

  • time (np.ndarray) – An sequence of (uniform) time points, usually created using np.arange or np.linspace.

  • H (np.ndarray) – Hamiltonian operator.

Returns:

Return a sequence of density matrices evolved through time, starting from a given initial state using the Hamiltonian H.

Return type:

np.ndarray

Examples

>>> molecules = [Molecule.fromdb("flavin_anion", ["N5"]),
...              Molecule("Z")]
>>> sim = HilbertSimulation(molecules)
>>> H = sim.total_hamiltonian(B0=0, J=0, D=0)
>>> time = np.arange(0, 2e-6, 5e-9)
>>> time.shape
(400,)
>>> rhos = sim.time_evolution(State.SINGLET, time, H)
>>> rhos.shape
(400, 12, 12)
total_hamiltonian(B0, J, D, theta=None, phi=None, hfc_anisotropy=False, kS=0.0, kT=0.0)[source]

Construct the total Hamiltonian as a matrix product operator.

Parameters:

  • B0 (float) – Magnetic field strength in mT.

  • J (float) – Exchange coupling constant in mT.

  • D (float | np.ndarray) – Dipolar coupling constant or tensor in mT.

  • theta (float, optional) – Polar angle of the magnetic field in radians.

  • phi (float, optional) – Azimuthal angle of the magnetic field in radians.

  • hfc_anisotropy (bool) – Whether to include anisotropic hyperfine coupling. Default is False.

  • kS (float) – Singlet recombination rate in s^-1. Default is 0.0.

  • kT (float) – Triplet recombination rate in s^-1. Default is 0.0.

Returns:

Total Hamiltonian as a matrix product operator.

Return type:

list[np.ndarray]

static unitary_propagator(H, dt)[source]

Unitary propagator is not used for BaseMPSSimulation. Use time_evolution() instead.

Return type:

tuple[ndarray, ndarray]

zeeman_hamiltonian(B0, B_axis='z', theta=None, phi=None)[source]

Construct the Zeeman Hamiltonian as a SumOfProducts.

Parameters:

  • B0 (float) – Magnetic field strength in mT.

  • B_axis (str) – Axis of the magnetic field (‘x’, ‘y’, or ‘z’). Default is ‘z’.

  • theta (float, optional) – Polar angle of the magnetic field in radians.

  • phi (float, optional) – Azimuthal angle of the magnetic field in radians.

Returns:

The Zeeman Hamiltonian as a sum of products of operators.

Return type:

SumOfProducts

zero_field_splitting_hamiltonian(D, E)[source]

Zero-field splitting Hamiltonian is not implemented.

Return type:

ndarray

class radicalpy.tensornetwork.LocallyPurifiedMPSSimulation(molecules, custom_gfactors=False, basis=Basis.ST, bond_dimension=16, integrator='arnoldi', backend='auto', jobname='lpmps')[source]

Bases: BaseMPSSimulation

Locally purified matrix product state simulation for radical pair dynamics.

This class implements locally purified matrix product state (LPMPS) simulations for mixed quantum states. LPMPS uses purification to represent mixed states as pure states in an enlarged Hilbert space, enabling efficient simulation of decoherence and relaxation processes.

The purification approach works by: 1. Doubling the physical Hilbert space with auxiliary degrees of freedom 2. Representing mixed states as pure states in the enlarged space 3. Propagating the purified state using tensor network methods 4. Tracing out auxiliary degrees of freedom to obtain physical observables

Parameters:

  • molecules (list[Molecule]) – List of molecules in the radical pair.

  • custom_gfactors (bool) – Whether to use custom g-factors. Default is False.

  • basis (Basis) – Basis set for the simulation. Only ST basis is supported.

  • bond_dimension (int) – Bond dimension for the matrix product state. Default is 16.

  • integrator (Literal["arnoldi", "lanczos"]) – Integration method for time evolution. Use ‘lanczos’ for Hermitian Hamiltonians, ‘arnoldi’ for non-Hermitian.

  • backend (Literal["numpy", "jax", "auto"]) – Computational backend. ‘auto’ selects ‘jax’ for bond_dimension > 32, ‘numpy’ otherwise.

  • jobname (str) – Job name for output files. Default is “lpmps”.

Note

Important differences from standard simulations:

  • Uses purification approach with auxiliary degrees of freedom (approximately twice the computational resources compared to pure state simulations)

  • Only State.SINGLET initial state is supported (other states raise NotImplementedError)

  • Returns diagonal elements of reduced density matrix instead of full density matrices

  • Requires proper site ordering for auxiliary degrees of freedom

  • Uses Hilbert space formulation with doubled physical space

Example

from radicalpy.tensornetwork import LocallyPurifiedMPSSimulation

sim = LocallyPurifiedMPSSimulation(
    molecules=molecules,
    bond_dimension=64,
    backend="jax"  # For GPU acceleration
)

result = sim.time_evolution(init_state=State.SINGLET, time=time, H=H)
time_evolution(init_state, time, H)[source]

Evolve the system through time.

Parameters:

  • init_state (State) – Initial State of the density matrix (see projection_operator).

  • time (np.ndarray) – An sequence of (uniform) time points, usually created using np.arange or np.linspace.

  • H (list[np.ndarray]) – Hamiltonian matrix product operator.

Returns:

Diagonal elements of the reduced density matrix.

Return type:

np.ndarray

total_hamiltonian(B0, J, D, theta=None, phi=None, hfc_anisotropy=False, kS=0.0, kT=0.0)[source]

Construct the total Hamiltonian as a matrix product operator.

Parameters:

  • B0 (float) – Magnetic field strength in mT.

  • J (float) – Exchange coupling constant in mT.

  • D (float | np.ndarray) – Dipolar coupling constant or tensor in mT.

  • theta (float, optional) – Polar angle of the magnetic field in radians.

  • phi (float, optional) – Azimuthal angle of the magnetic field in radians.

  • hfc_anisotropy (bool) – Whether to include anisotropic hyperfine coupling. Default is False.

  • kS (float) – Singlet recombination rate in s^-1. Default is 0.0.

  • kT (float) – Triplet recombination rate in s^-1. Default is 0.0.

Returns:

Total Hamiltonian as a matrix product operator.

Return type:

list[np.ndarray]

class radicalpy.tensornetwork.MPDOSimulation(molecules, custom_gfactors=False, basis=Basis.ST, bond_dimension=16, integrator='arnoldi', backend='auto', jobname='vmpdo')[source]

Bases: BaseMPSSimulation

Vectorised matrix product density operator simulation for radical pair dynamics.

This class implements vectorised matrix product density operator (VMPDO) simulations for mixed quantum states. VMPDO directly represents density operators as matrix product operators, enabling efficient simulation of decoherence, relaxation, and other open quantum system dynamics.

The VMPDO approach works by: 1. Reshape matrix product density operator into matrix product state 2. Propagating the density operator using Liouville space dynamics 3. Computing observables as traces of operator products 4. Supporting both Hermitian and non-Hermitian evolution

Parameters:

  • molecules (list[Molecule]) – List of molecules in the radical pair.

  • custom_gfactors (bool) – Whether to use custom g-factors. Default is False.

  • basis (Basis) – Basis set for the simulation. Only ST basis is supported.

  • bond_dimension (int) – Bond dimension for the matrix product operator. Default is 16.

  • integrator (Literal["arnoldi", "lanczos"]) – Integration method for time evolution. Use ‘lanczos’ for Hermitian Hamiltonians, ‘arnoldi’ for non-Hermitian.

  • backend (Literal["numpy", "jax", "auto"]) – Computational backend. ‘auto’ selects ‘jax’ for bond_dimension > 32, ‘numpy’ otherwise.

  • jobname (str) – Job name for output files. Default is “vmpdo”.

Note

  • Returns reduced density matrix instead of full density matrice

  • Trace preserving, positivity preserving, and completely positive evolution is not guaranteed

Example

from radicalpy.tensornetwork import MPDOSimulation

sim = MPDOSimulation(
    molecules=molecules,
    bond_dimension=64,
    backend="jax"  # For GPU acceleration
)

result = sim.time_evolution(init_state=State.SINGLET, time=time, H=H)
dipolar_hamiltonian(D)[source]

Construct the dipolar Hamiltonian as a SumOfProducts.

Parameters:

D (float | np.ndarray) – Dipolar coupling constant or tensor in mT.

Returns:

The dipolar Hamiltonian as a sum of products of operators.

Return type:

SumOfProducts

exchange_hamiltonian(J, prod_coeff=2)[source]

Construct the exchange Hamiltonian as a SumOfProducts.

Parameters:

  • J (float) – Exchange coupling constant in mT.

  • prod_coeff (float) – Coefficient for the S1·S2 term. Default is 2.

Returns:

The exchange Hamiltonian as a sum of products of operators.

Return type:

SumOfProducts

haberkorn_hamiltonian(kS, kT)[source]

Construct the Haberkorn recombination term as a SumOfProducts.

Parameters:

  • kS (float) – Singlet recombination rate in s^-1.

  • kT (float) – Triplet recombination rate in s^-1.

Returns:

The Haberkorn recombination Hamiltonian as a sum

of products of operators.

Return type:

SumOfProducts

hyperfine_hamiltonian(hfc_anisotropy=False)[source]

Construct the hyperfine Hamiltonian as a SumOfProducts.

Parameters:

hfc_anisotropy (bool) – Whether to include anisotropic hyperfine coupling. Default is False.

Returns:

The hyperfine Hamiltonian as a sum of products of operators.

Return type:

SumOfProducts

time_evolution(init_state, time, H)[source]

Evolve the system through time.

Parameters:

  • init_state (State) – Initial State of the density matrix (see projection_operator).

  • time (np.ndarray) – An sequence of (uniform) time points, usually created using np.arange or np.linspace.

  • H (list[np.ndarray]) – Hamiltonian matrix product operator.

Returns:

Reduced density matrix with shape (nsteps, 4, 4).

Return type:

np.ndarray

zeeman_hamiltonian(B0, B_axis='z', theta=None, phi=None)[source]

Construct the Zeeman Hamiltonian as a SumOfProducts.

Parameters:

  • B0 (float) – Magnetic field strength in mT.

  • B_axis (str) – Axis of the magnetic field (‘x’, ‘y’, or ‘z’). Default is ‘z’.

  • theta (float, optional) – Polar angle of the magnetic field in radians.

  • phi (float, optional) – Azimuthal angle of the magnetic field in radians.

Returns:

The Zeeman Hamiltonian as a sum of products of operators.

Return type:

SumOfProducts

class radicalpy.tensornetwork.StochasticMPSSimulation(molecules, custom_gfactors=False, basis=Basis.ST, bond_dimension=16, integrator='arnoldi', nsamples=128, max_workers=8, jobname='smps')[source]

Bases: BaseMPSSimulation

Stochastic matrix product state simulation for radical pair dynamics.

This class implements stochastic matrix product state (SMPS) simulations using Monte Carlo sampling to generate trajectories of radical pair states. SMPS is particularly effective for high-dimensional quantum systems where exact simulation becomes computationally prohibitive.

The simulation works by: 1. Sampling from spin coherent states for nuclear spins 2. Propagating each sample using tensor network methods 3. Averaging over all samples to obtain ensemble properties

nsamples

Number of Monte Carlo samples to generate.

Type:

int

max_workers

Maximum number of parallel workers for sampling.

Type:

int

Parameters:

  • molecules (list[Molecule]) – List of molecules in the radical pair.

  • custom_gfactors (bool) – Whether to use custom g-factors. Default is False.

  • basis (Basis) – Basis set for the simulation. Only ST basis is supported.

  • bond_dimension (int) – Bond dimension for the matrix product state. Default is 16.

  • integrator (Literal["arnoldi", "lanczos"]) – Integration method for time evolution. Use ‘lanczos’ for Hermitian Hamiltonians, ‘arnoldi’ for non-Hermitian.

  • nsamples (int) – Number of Monte Carlo samples to generate. Default is 128.

  • max_workers (int) – Maximum number of parallel workers for sampling. Default is 8.

  • jobname (str) – Job name for output files. Default is “smps”.

Raises:

RuntimeError – If thread number environment variables are not set.

Note

Important differences from standard simulations:

  • Only State.SINGLET initial state is supported (other states raise NotImplementedError)

  • Returns diagonal elements of reduced density matrix instead of full density matrices

  • Uses Monte Carlo sampling with statistical averaging over nsamples trajectories

  • Requires proper thread configuration. Set one of the following environment variables before importing: OMP_NUM_THREADS, OPENBLAS_NUM_THREADS, MKL_NUM_THREADS, VECLIB_MAXIMUM_THREADS, or NUMEXPR_NUM_THREADS

  • Uses multiprocessing with fork method for parallel execution

Example

import os
os.environ["OMP_NUM_THREADS"] = "1"

from radicalpy.tensornetwork import StochasticMPSSimulation

sim = StochasticMPSSimulation(
    molecules=molecules,
    bond_dimension=32,
    nsamples=256,
    max_workers=4
)

diagonal_elements = sim.time_evolution(init_state=State.SINGLET, time=time, H=H)
time_evolution(init_state, time, H)[source]

Evolve the system through time using stochastic matrix product states.

This method performs time evolution using Monte Carlo sampling of spin coherent states. Each sample is propagated independently using tensor network methods, and the results are averaged to obtain ensemble properties.

The simulation process: 1. Generates random spin coherent states for nuclear spins 2. Propagates each sample using the provided Hamiltonian 3. Computes reduced density matrix for each sample 4. Averages over all samples to obtain ensemble properties

Parameters:

  • init_state (State) – Initial state of the system. Only State.SINGLET is currently supported for stochastic simulations.

  • time (np.ndarray) – Array of time points for the evolution, typically created using np.linspace or np.arange. Must be uniformly spaced.

  • H (list[np.ndarray]) – Hamiltonian as a matrix product operator. Each element represents a tensor in the MPO.

Returns:

Diagonal elements of the reduced density matrix with

shape (nsteps, 4) representing [T+, T0, S, T-] populations.

Return type:

np.ndarray

Raises:

NotImplementedError – If init_state is not State.SINGLET.

Note

Critical differences from standard time evolution:

  • Returns diagonal elements of reduced density matrix (shape: nsteps, 4) instead of full density matrices (shape: nsteps, dim, dim)

  • Only State.SINGLET initial state is supported

  • Uses stochastic sampling over nuclear spin coherent states

  • Requires statistical averaging over multiple Monte Carlo samples

  • Uses parallel processing with multiprocessing for efficiency

Example

import numpy as np
from radicalpy import State

# Create time array
time = np.linspace(0, 1e-6, 1000)

# Run simulation
result = sim.time_evolution(
    init_state=State.SINGLET,
    time=time,
    H=hamiltonian
)

# result.shape is (1000, 4)
# result[:, 2] contains singlet population

radicalpy.utils module

Utilities for spin dynamics, magnetic‐resonance and molecular simulations.

This module collects small, reusable helpers used across the codebase: unit conversions, geometry on spheres, spectrum building blocks, I/O parsers for molecular files and quantum-chemistry outputs, and a few domain-specific utilities for MARY/NMR workflows.

Main contents:

Constants

  • COVALENT_RADII: Covalent radii (Å) for common elements; used by

bonding heuristics and label/element inference.

CLI/testing

  • is_fast_run(): Check --fast CLI flag to run lighter examples/tests.

Fitting & line shapes

  • Bhalf_fit(B, MARY): Fit MARY data to a Lorentzian to obtain

\(B_{1/2}\) and goodness-of-fit. - Lorentzian(B, amplitude, Bhalf): MARY Lorentzian model.

Unit conversions

  • Gauss_to_mT/MHz/…, MHz_to_mT/Gauss/…,

mT_to_MHz/Gauss/…, angular_frequency_in_MHz, *_to_angular_frequency: Consistent conversions between G, mT, MHz and rad·s⁻¹·T⁻¹. Consistent conversions between GHz, J, meV, mK.

Angular grids & spherical geometry

  • anisotropy_check(theta, phi): Validate angle grids and parity (θ odd/φ even).

  • _check_full_sphere(theta, phi): Ensure full-sphere linspace grids.

  • spherical_average(...): Weighted average over a θ/φ grid.

  • cartesian_to_spherical(...), spherical_to_cartesian(...):

Coordinate transforms.

Autocorrelation

  • autocorrelation(data, factor=1): FFT-based trajectory autocorrelation.

CIDNP helper models

  • cidnp_polarisation_*: Exponential, truncated diffusion, and full

diffusion models; see docstrings for details. - s_t0_omega(deltag, B0, hfc_star, onuc_all): ω₊/ω₋ for S–T₀ mixing.

NMR building blocks

  • nmr_chemical_shift_real/im...: cos/sin(2π f t) modulators.

  • nmr_scalar_coupling_modulation: J-modulation term.

  • nmr_t2_relaxation: exp(−t/T₂) envelope.

Lock-in MARY helpers

  • modulated_signal(...), reference_signal(...):

Time-domain modulation and reference waves. - mary_lorentzian(mod_signal, lfe_magnitude): Simple MARY line shape.

Molecular geometry & visualization

  • define_xyz(...): Construct an orthonormal (x,y,z) basis from

atom/point triplets. - get_angle_between_plane(A, B): Angle (deg) between plane normals. - get_rotation_matrix_euler_angles(A, B): Rotation matrix and ZXZ Euler angles mapping basis A→B. - rodrigues_rotation(v, k, theta): Rotate vectors around axis k by θ. - rotate_axes(A, x, y, z): Apply a rotation matrix to a basis. - enumerate_spin_states_from_base(base): Mixed-radix enumeration of spin projections (e.g., 2 for spin-½, 3 for spin-1). - infer_bonds(elements, coords, ...): Heuristic covalent bonding.

File I/O & parsing

  • parse_xyz(path), parse_label_xyz_txt(path), parse_pdb(path, ...):

Read coordinates and optional bonds/labels for plotting. - pdb_label(atom, scheme=...): Human-readable PDB atom labels. - write_xyz/mol_to_plot_arrays/write_pdb/write_sdf: Minimal writers and plotting helpers. - read_trajectory_files(dir, scale): Concatenate text trajectories from a folder.

Quantum-chemistry integration (ORCA)

  • read_orca_hyperfine(path, version=6): Dispatch to ORCA v6 parsers.

  • Internal: _hyperfine_from_orca6_out / _hyperfine_from_orca6_property_txt.

Return zero-based nucleus indices, isotope labels, and 3×3 HFC tensors (mT). - read_lines_utf8_or_utf16le: Robust text decoding for ORCA outputs. - write_orca_from_pdb: Creates customisable ORCA input files.

Quantum quantification methods
  • negativity: Computes the (logarithmic) negativity of a multipartite density matrix.

    A measure of entanglement.

  • purity: Calculates the purity of a density matrix.

  • von_neumann_entropy: Computes the von Neumann entropy of a quantum state.

Chemistry utilities (RDKit)

  • smiles_to_3d(smiles, add_h=True, opt='mmff'): Generate a 3D conformer

(ETKDG) with optional MMFF/UFF optimisation. - mol_to_plot_arrays(mol): Extract elements, labels, coords and bonds.

Spectral density

  • spectral_density(omega, tau_c): \(J(ω) = τ_c / (1 + ω^2 τ_c^2)\).

Notes

  • Many routines are vectorised and expect NumPy arrays; see individual

docstrings for shapes and units. - Angle conventions use radians: theta ∈ [0, π], phi ∈ [0, 2π]. - Unit conversions rely on physical constants from radicalpy.constants (C).

radicalpy.utils.Bhalf_LFEhalf_fit(B, MARY, LFE_position=1.0)[source]

B_1/2 and LFE_1/2 fitting for MARY spectra.

Parameters:

  • B (np.ndarray) – Magnetic field values (x-axis).

  • MARY (np.ndarray) – Magnetic field effect data (y-axis). Use the MARY entry in the result of radicalpy.simulation.HilbertSimulation.MARY.

  • LFE_position (float) – Initial guess for the low field effect position.

Returns:

  • Bhalf (float): The magnetic field strength at half the saturation magnetic field.

  • LFEhalf (float): The magnetic field strength at half the low field effect magnetic field.

  • fit_result (np.ndarray): y-axis from fit.

  • fit_error (float): Standard error for the fit.

  • R2 (float): R-squared value for the fit.

Return type:

(float, np.ndarray, float, float)

radicalpy.utils.Bhalf_fit(B, MARY)[source]

B_1/2 fit for MARY spectra.

Parameters:

  • B (np.ndarray) – Magnetic field values (x-axis).

  • MARY (np.ndarray) – Magnetic field effect data (y-axis). Use the MARY entry in the result of radicalpy.simulation.HilbertSimulation.MARY.

Returns:

  • Bhalf (float): The magnetic field strength at half the saturation magnetic field.

  • fit_result (np.ndarray): y-axis from fit.

  • fit_error (float): Standard error for the fit.

  • R2 (float): R-squared value for the fit.

Return type:

(float, np.ndarray, float, float)

radicalpy.utils.GHz_to_mK(GHz)[source]

Convert a frequency from GHz to mK using h ν = k_B T.

Relation:

h · (1e9 · ν_GHz) = k_B · (1e-3 · T_mK) ⇒ T_mK = 1e12 · (h / k_B) · ν_GHz

Parameters:

GHz (float or ndarray) – Value(s) in GHz.

Returns:

Value(s) in mK.

Return type:

float or ndarray

radicalpy.utils.GHz_to_meV(GHz)[source]

Convert an energy/frequency from GHz to meV.

Uses 1 GHz = 4.1357×10⁻³ meV.

Parameters:

GHz (float) – Value(s) in GHz.

Returns:

Value(s) in meV.

Return type:

float

radicalpy.utils.Gauss_to_MHz(Gauss)[source]

Convert Gauss to MHz.

Parameters:

Gauss (float) – The magnetic flux density in Gauss (G).

Returns:

The magnetic flux density converted to MHz.

Return type:

float

radicalpy.utils.Gauss_to_angular_frequency(Gauss)[source]

Convert Gauss to angular frequency.

Parameters:

Gauss (float) – The magnetic flux density in Gauss (G).

Returns:

The magnetic flux density converted to angular frequency (rad/s/T).

Return type:

float

radicalpy.utils.Gauss_to_mT(Gauss)[source]

Convert Gauss to millitesla.

Parameters:

Gauss (float) – The magnetic flux density in Gauss (G).

Returns:

The magnetic flux density converted to millitesla (mT).

Return type:

float

radicalpy.utils.J_to_meV(J)[source]

Convert an energy from Joule (J) to meV.

Uses 1 eV = 1.602176565×10⁻¹⁹ J.

Parameters:

J (float or ndarray) – Value(s) in Joule.

Returns:

Value(s) in meV.

Return type:

float or ndarray

radicalpy.utils.Lorentzian(B, amplitude, Bhalf)[source]

Lorentzian function for MARY spectra.

More information in radicalpy.utils.Bhalf_fit (where this is used).

Parameters:

  • B (np.ndarray) – The x-axis magnetic field values.

  • amplitude (float) – The amplitude of the saturation field value.

  • Bhalf (float) – The magnetic field strength at half the saturation field value.

Returns:

Lorentzian function for MARY spectrum.

Return type:

np.ndarray

radicalpy.utils.MHz_in_angular_frequency(MHz)[source]

Convert MHz into angular frequency.

Parameters:

MHz (float) – The angular frequency in MHz/T

Returns:

The angular frequency converted to rad/s/T.

Return type:

float

radicalpy.utils.MHz_to_Gauss(MHz)[source]

Convert Megahertz to Gauss.

Parameters:

MHz (float) – The frequency in Megahertz (MHz).

Returns:

Megahertz (MHz) converted to Gauss (G).

Return type:

float

radicalpy.utils.MHz_to_mT(MHz)[source]

Convert Megahertz to milltesla.

Parameters:

MHz (float) – The frequency in Megahertz (MHz).

Returns:

Megahertz (MHz) converted to millitesla (mT).

Return type:

float

radicalpy.utils.angular_frequency_in_MHz(ang_freq)[source]

Convert angular frequency into MHz.

Parameters:

ang_freq (float) – The angular frequency in rad/s/T.

Returns:

The angular frequency converted to MHz/T.

Return type:

float

radicalpy.utils.angular_frequency_to_Gauss(ang_freq)[source]

Convert angular frequency to Gauss.

Parameters:

ang_freq (float) – The angular frequency in rad/s/T.

Returns:

The angular frequency converted to Gauss (G).

Return type:

float

radicalpy.utils.angular_frequency_to_mT(ang_freq)[source]

Convert angular frequency to millitesla.

Parameters:

ang_freq (float) – The angular frequency in rad/s/T.

Returns:

The angular frequency converted to millitesla (mT).

Return type:

float

radicalpy.utils.anisotropy_check(theta, phi)[source]

Validate and normalise anisotropy angle grids.

Accepts scalar or array inputs for the polar (theta) and azimuthal (phi) angles, promotes scalars to 1-element arrays, and enforces domain and parity constraints that downstream integrators rely on.

Parameters:

  • theta (float or np.ndarray) – Polar angle(s) in radians. Must lie in [0, π]. If an array is provided and len(theta) > 1 while len(phi) > 1, then len(theta) must be odd.

  • phi (float or np.ndarray) – Azimuthal angle(s) in radians. Must lie in [0, 2π]. If an array is provided and len(theta) > 1 while len(phi) > 1, then len(phi) must be even.

Returns:

The validated angle arrays (theta, phi) (possibly promoted from scalars).

Return type:

(np.ndarray, np.ndarray)

Raises:

ValueError – If angles are out of range or the grid parity constraints are violated.

radicalpy.utils.autocorrelation(data, factor=1)[source]

Calculate the autocorrelation of a trajectory.

An FFT-based implementation of the autocorrelation for Monte Carlo or molecular dynamics trajectories (or any other time dependent value).

Parameters:

  • data (np.ndarray) – The time dependent trajectory.

  • factor (int) – Data length reduction factor.

Returns:

The autocorrelation of the trajectory.

Return type:

np.ndarray

radicalpy.utils.cartesian_to_spherical(x, y, z)[source]

Convert Cartesian coordinates to spherical coordinates.

Parameters:

  • x (float or np.ndarray) – Coordinate(s) in the x plane.

  • y (float or np.ndarray) – Coordinate(s) in the y plane.

  • z (float or np.ndarray) – Coordinate(s) in the z plane.

Returns:

The radial distance(s). theta (float or np.ndarray): The polar angle(s). phi (float or np.ndarray): The azimuthal angle(s).

Return type:

r (float or np.ndarray)

radicalpy.utils.cidnp_polarisation_diffusion_model(omega_plus, omega_minus, alpha=1.5)[source]

Compute the CIDNP polarisation using the full diffusion model.

Parameters:

  • omega_plus (np.ndarray) – Array of omega+ frequencies (rad/s).

  • omega_minus (np.ndarray) – Array of omega- frequencies (rad/s).

  • alpha (float) – Dimensionless parameter for the Adrian diffusion model 2p/m.

Returns:

CIDNP polarisation.

Return type:

p (float)

radicalpy.utils.cidnp_polarisation_exponential_model(ks, omega_plus, omega_minus)[source]

Compute the CIDNP polarisation using the exponential model.

Parameters:

  • ks (float) – Singlet recombination rate (s^-1).

  • omega_plus (np.ndarray) – Array of omega+ frequencies (rad/s).

  • omega_minus (np.ndarray) – Array of omega- frequencies (rad/s).

Returns:

CIDNP polarisation.

Return type:

p (float)

radicalpy.utils.cidnp_polarisation_truncated_diffusion_model(omega_plus, omega_minus)[source]

Compute the CIDNP polarisation using the truncated t^{-3/2} diffusion model.

Parameters:

  • omega_plus (np.ndarray) – Array of omega+ frequencies (rad/s).

  • omega_minus (np.ndarray) – Array of omega- frequencies (rad/s).

Returns:

CIDNP polarisation.

Return type:

p (float)

radicalpy.utils.define_xyz(x1, x2, z1, z2, z3, z4)[source]

Construct a right-handed orthonormal basis from atom/point triplets.

The z-axis is defined as the normal to the plane spanned by the vectors (z1 − z2) and (z3 − z4) (right-hand rule). The x-axis is the normalised direction from x2 to x1. The y-axis completes the right-handed triad via y = z × x, and x is re-orthogonalised by x = y × z.

Parameters:

  • x1 (array-like) – Points defining the x-axis direction.

  • x2 (array-like) – Points defining the x-axis direction.

  • z1 (array-like) – Points defining two in-plane vectors whose cross product yields the z-axis.

  • z2 (array-like) – Points defining two in-plane vectors whose cross product yields the z-axis.

  • z3 (array-like) – Points defining two in-plane vectors whose cross product yields the z-axis.

  • z4 (array-like) – Points defining two in-plane vectors whose cross product yields the z-axis.

Returns:

Unit vectors (x, y, z) forming a right-handed orthonormal basis.

Return type:

(np.ndarray, np.ndarray, np.ndarray)

radicalpy.utils.double_Lorentzian(B, amplitude, Bhalf, LFE_amplitude, LFEhalf)[source]

Double Lorentzian function for MARY spectra.

More information in radicalpy.utils.Bhalf_fit (where this is used).

Parameters:

  • B (np.ndarray) – The x-axis magnetic field values.

  • amplitude (float) – The amplitude of the saturation field value.

  • Bhalf (float) – The magnetic field strength at half the saturation field value.

  • LFE_amplitude (float) – The amplitude of the LFE component.

  • LFEhalf (float) – The magnetic field strength at half the LFE field value.

Returns:

Double Lorentzian function for MARY spectrum.

Return type:

np.ndarray

radicalpy.utils.eigensorter(H)[source]

Diagonalise a square matrix and return eigenvalues in ascending order together with the corresponding (row-stacked) eigenvectors.

This helper wraps numpy.linalg.eig(), sorts the eigenpairs by ascending eigenvalue, and returns the eigenvectors transposed so that each eigenvector appears as a row in the returned array.

Parameters:

H (ndarray of shape (N, N), complex or float) – Matrix to diagonalise.

:param For Hermitian inputs consider using numpy.linalg.eigh(): :param elsewhere for improved numerical stability: :param but this routine: :param intentionally uses eig to support non-Hermitian cases.:

Returns:

Eigenvalues sorted in ascending order.

evecs (ndarray of shape (N, N)): Row-stacked eigenvectors corresponding to evals. Each row k is the eigenvector v_k satisfying (approximately) H @ v_k = evals[k] * v_k.

Return type:

evals (ndarray of shape (N,))

radicalpy.utils.enumerate_spin_states_from_base(base)[source]

Return all spin-state patterns for a mixed-radix ‘base’ (e.g. [2,2,3,…]). Each row corresponds to one configuration. For base[i]=b, the digit d in [0..b-1] maps to spin projection m = (b-1)/2 - d.

Return type:

ndarray

radicalpy.utils.get_angle_between_plane(A, B)[source]

Angle between two plane normals in degrees.

Computes the angle between vectors A and B (interpreted as plane normals). The result is within [0, 180] degrees.

Parameters:

  • A (array-like) – First normal vector.

  • B (array-like) – Second normal vector.

Returns:

Angle between A and B in degrees.

Return type:

float

radicalpy.utils.get_rotation_matrix_euler_angles(A, B)[source]

Rotation matrix and ZXZ Euler angles mapping basis A → basis B.

Forms the rotation matrix R = Aᵀ B that maps coordinates expressed in the (column) basis A to the basis B, and extracts the ZXZ Euler angles (α, β, γ) in radians.

Parameters:

  • A (np.ndarray) – 3×3 column-basis matrix of frame A.

  • B (np.ndarray) – 3×3 column-basis matrix of frame B.

Returns:

  • R: 3×3 rotation matrix.

  • alpha (rad): First ZXZ angle.

  • beta (rad): Second ZXZ angle.

  • gamma (rad): Third ZXZ angle.

Return type:

(np.ndarray, float, float, float)

Notes

Handles the gimbal-lock case |R[2,2]| = 1 explicitly.

radicalpy.utils.hilbert_to_liouville(H)[source]

Liouvillian for unitary part in the basis of H (H is already in that basis): L[rho] = -i [H, rho] -> vec form: -i (I ⊗ H - H.T ⊗ I)

radicalpy.utils.infer_bonds(elements, coords, scale=1.2, max_dist=2.0)[source]

Heuristic covalent-bond detection from interatomic distances.

Bonds are inferred when the interatomic distance is less than scale × (r_i + r_j) using covalent radii (fallback 0.77 Å) and also less than max_dist to avoid spurious long bonds.

Parameters:

  • elements (list[str]) – Atomic symbols per atom.

  • coords (np.ndarray) – Cartesian coordinates, shape (N, 3) (Å).

  • scale (float) – Multiplicative tolerance on summed covalent radii.

  • max_dist (float) – Hard cutoff distance (Å).

Returns:

Pairs of atom indices (i, j) indicating bonds.

Return type:

list[tuple[int, int]]

radicalpy.utils.is_fast_run()[source]

Is the --fast parameter set at execution.

This function helps examples to be used as tests. By running the example with the --fast option, a faster version of main can be called (e.g., by setting fewer number of time steps etc.).

radicalpy.utils.mK_to_GHz(mK)[source]

Convert a temperature-like energy from mK to GHz using h ν = k_B T.

Inverse of convert_GHz_to_mK().

Parameters:

mK (float or ndarray) – Value(s) in mK.

Returns:

Value(s) in GHz.

Return type:

float or ndarray

radicalpy.utils.mK_to_meV(mK)[source]

Convert an energy from mK to meV.

Uses 1 mK = 8.61740×10⁻⁵ meV.

Parameters:

mK (float or ndarray) – Value(s) in mK.

Returns:

Value(s) in meV.

Return type:

float or ndarray

radicalpy.utils.mT_to_Gauss(mT)[source]

Convert millitesla to Gauss.

Parameters:

mT (float) – The magnetic flux density in millitesla (mT).

Returns:

The magnetic flux density converted to Gauss (G).

Return type:

float

radicalpy.utils.mT_to_MHz(mT)[source]

Convert millitesla to Megahertz.

Parameters:

mT (float) – The magnetic flux density in millitesla (mT).

Returns:

The magnetic flux density converted to Megahertz (MHz).

Return type:

float

radicalpy.utils.mT_to_angular_frequency(mT)[source]

Convert millitesla to angular frequency.

Parameters:

mT (float) – The magnetic flux density in millitesla (mT).

Returns:

The magnetic flux density converted to angular frequency (rad/s/T).

Return type:

float

radicalpy.utils.make_resonance_sticks(*, bins=1000, freq=9373000000.0, gval=2.0023, hfcH=(), hfcN=())[source]

Stick-field/intensity generator.

Returns:

  • fields_gauss ((M,) ndarray) – Non-zero histogram bin centers, in GAUSS.

  • intensities ((M,) ndarray) – Corresponding normalized intensities (sum = 1).

radicalpy.utils.mary_lorentzian(mod_signal, lfe_magnitude)[source]

Lorentzian MARY spectral shape.

Used for brute force modulated MARY simulations.

Parameters:

  • mod_signal (np.ndarray) – The modulated signal.

  • lfe_magnitude (float) – The magnitude of your low field effect (G).

Returns:

The modulated MARY signal.

Return type:

np.ndarray

radicalpy.utils.matrix_to_vector(M)[source]

Convert a matrix into a column vector.

radicalpy.utils.meV_to_GHz(meV)[source]

Convert an energy from meV to GHz.

Uses 1 meV = 1 / (4.1357×10⁻³) GHz.

Parameters:

meV (float or ndarray) – Value(s) in meV.

Returns:

Value(s) in GHz.

Return type:

float or ndarray

radicalpy.utils.meV_to_J(meV)[source]

Convert an energy from meV to Joules (J).

Uses 1 eV = 1.602176565×10⁻¹⁹ J.

Parameters:

meV (float or ndarray) – Value(s) in meV.

Returns:

Value(s) in Joules.

Return type:

float or ndarray

radicalpy.utils.meV_to_mK(meV)[source]

Convert an energy from meV to mK.

Uses 1 mK = 8.61740×10⁻⁵ meV.

Parameters:

meV (float or ndarray) – Value(s) in meV.

Returns:

Value(s) in mK.

Return type:

float or ndarray

radicalpy.utils.modulated_signal(timeconstant, theta, frequency)[source]

Modulated MARY signal.

Used for brute force modulated MARY simulations.

Parameters:

  • timeconstant (np.ndarray) – The modulation time constant (s).

  • theta (float) – The modulation phase (rad).

  • frequency (float) – The modulation frequency (Hz).

Returns:

The modulated signal.

Return type:

np.ndarray

radicalpy.utils.mol_to_plot_arrays(mol)[source]

Return labels, elements, coords[N,3], bonds[(i,j)] from an RDKit Mol with a conformer.

radicalpy.utils.negativity(rho, dims, subsys, method='tracenorm', logarithmic=False)[source]

Compute the (logarithmic) negativity of a multipartite density matrix.

This is a RadicalPy-compatible reimplementation of QuTiP’s negativity, using plain NumPy arrays. The partial transpose is taken over the subsystem(s) specified by subsys.

Parameters:

  • rho (ndarray of shape (N, N)) – Density matrix in the computational (lab)

  • prod (basis. Must be square with N =)

  • dims (Sequence[int]) – Local Hilbert-space dimensions for each subsystem,

  • [2 (or)

  • qubits (2] for two)

  • [2

  • 3

  • system. (2] for a 2×3×2)

  • N. (The product must equal)

  • subsys (int or Sequence[int]) – Index or indices of the subsystem(s)

  • example (on which to perform the partial transpose. For)

  • subsys=0

  • [d0 (on a bipartite)

  • subsystem. (d1] system transposes the first)

  • method ({"tracenorm", "eigenvalues"}, optional) –

    • “tracenorm” (default): use ||ρ^{T_A}||_1 via SVD (stable).

    • ”eigenvalues”: sum of negative eigenvalues of ρ^{T_A}.

  • logarithmic (bool, optional) – If True, return the logarithmic

  • log2 (negativity)

Returns:

Negativity (or logarithmic negativity if logarithmic=True).

Return type:

float

radicalpy.utils.nmr_chemical_shift_imaginary_modulation(freq_hz, t)[source]

Chemical-shift modulation (imag): sin(2π f t).

Return type:

ndarray

radicalpy.utils.nmr_chemical_shift_real_modulation(freq_hz, t)[source]

Chemical-shift modulation (real): cos(2π f t).

Return type:

ndarray

radicalpy.utils.nmr_scalar_coupling_modulation(j_hz, t, mult_minus_one)[source]

J-modulation: cos(π J t) ** (mult - 1).

Return type:

ndarray

radicalpy.utils.nmr_t2_relaxation(t, t2)[source]

T2 relaxation decay: exp(-t / T2).

Return type:

ndarray

radicalpy.utils.parse_label_xyz_txt(path)[source]

Parse a simple labeled XYZ-like text file.

Expected format per non-comment line: <label> <x> <y> <z> with whitespace separation. Lines starting with # are ignored.

Parameters:

path (str | Path) – Path to the labeled coordinate file.

Returns:

  • labels: Raw labels from file.

  • elements: Guessed elements from label alphabetic prefix.

  • coords: Cartesian coordinates, shape (N, 3) (Å).

Return type:

(list[str], list[str], np.ndarray)

Raises:

ValueError – If a line cannot be parsed into 4 fields.

radicalpy.utils.parse_pdb(path, use_rdkit_bonds=False, label_scheme='atom')[source]

Load a PDB and return plotting arrays (labels, elements, coords, bonds).

Builds coordinates from an existing conformer or embeds one if missing. Optionally uses RDKit’s bond topology.

Parameters:

  • path (str | Path) – Path to a .pdb file.

  • use_rdkit_bonds (bool) – If True, bonds are taken from RDKit connectivity; otherwise an empty list is returned.

  • label_scheme (str) – One of {"atom","res_atom","atom_serial","chain_res_atom"} passed to pdb_label().

Returns:

  • labels: Per-atom labels.

  • elements: Atomic symbols.

  • coords: Cartesian coordinates, shape (N, 3) (Å).

  • bonds: Pairs of atom indices (optional if use_rdkit_bonds=False).

Return type:

(list[str], list[str], np.ndarray, list[tuple[int,int]])

radicalpy.utils.parse_xyz(path)[source]

Parse a standard XYZ file (single frame).

The file must contain at least the header line with atom count, a comment line, and N subsequent atom lines.

Parameters:

path (str | Path) – Path to the .xyz file.

Returns:

  • labels: Auto-generated labels <element><index>.

  • elements: Atomic symbols.

  • coords: Cartesian coordinates, shape (N, 3) (Å).

Return type:

(list[str], list[str], np.ndarray)

Raises:

ValueError – If the file is too short or the atom count does not match.

radicalpy.utils.pdb_label(atom, scheme='chain_res_atom')[source]

Create a human-readable atom label from PDB residue info.

Parameters:

  • atom (rdkit.Chem.Atom) – RDKit atom with PDB residue info.

  • scheme (str) – Labeling scheme: - "atom" → atom name or element symbol. - "res_atom" → RESN<resid>:ATOM. - "atom_serial" → ATOM(serial). - "chain_res_atom" (default) → CHAIN:RESN<resid>:ATOM.

Returns:

Formatted label; falls back to element symbol if fields are missing.

Return type:

str

radicalpy.utils.purity(rho)[source]

Calculate the purity of a density matrix.

The purity is defined as \(P(\rho) = \operatorname{Tr}(\rho^2)\). It quantifies how mixed a quantum state is: pure states have \(P = 1\), while a maximally mixed state in a Hilbert space of dimension \(d\) has \(P = 1/d\).

Parameters:

  • rho (ndarray of shape (N, N)) – Density matrix of the quantum state.

  • and (The matrix is expected to be square)

  • trace (with unit)

  • inside (although these conditions are not enforced)

  • function. (the)

Returns:

Purity of the state, computed as real(trace(rho @ rho)).

Return type:

float

radicalpy.utils.read_lines_utf8_or_utf16le(path)[source]

Read text file as UTF-8 (with BOM) or UTF-16LE, returning lines.

Tries UTF-8 (BOM-aware) first, then UTF-16LE, which is common in ORCA outputs on Windows.

Parameters:

path (Path) – File path.

Returns:

Lines of the decoded file (including newlines).

Return type:

list[str]

Raises:

UnicodeError – If the file is neither valid UTF-8 nor UTF-16LE.

radicalpy.utils.read_orca_hyperfine(path, version=6)[source]

Read ORCA output file.

Parameters:

path (Path) – The path to the ORCA output file. Both .out and .property.txt are supported.

Returns:

The list of nucleus indices, isotopes, and HFC matrices. Note that indices starts from 0 (same as ORCA).

Return type:

tuple[list[int], list[str], list[np.ndarray]]

Examples

>>> import radicalpy as rp
>>> indices, isotopes, hfc_matrices = rp.utils.read_orca("orca.out") # or rp.utils.read_orca("orca.property.txt")
>>> nuclei = [rp.data.Nucleus.fromisotope(isotope, hfc_matrix.tolist()) for isotope, hfc_matrix in zip(isotopes, hfc_matrices)]
>>> molecule = rp.simulation.Molecule("mol1", nuclei)
Molecule: mol1
Nuclei:
14N(19337792.0, 3, 1.017 <anisotropic available>)
1H(267522187.44, 2, -2.199 <anisotropic available>)
1H(267522187.44, 2, -2.198 <anisotropic available>)
Radical: E(-176085963023.0, 2, 0.0 <anisotropic not available>)
radicalpy.utils.read_trajectory_files(path, scale=1e-10)[source]

Load and concatenate text trajectory files in a directory.

Each file under path is read with numpy.genfromtxt and concatenated along the first axis; the result is multiplied by scale (default converts Å → m if files are in Å).

Parameters:

  • path (Path) – Directory containing trajectory fragments.

  • scale (float) – Multiplicative scaling factor applied post-concat.

Returns:

Concatenated (and scaled) trajectory array.

Return type:

np.ndarray

radicalpy.utils.reference_signal(timeconstant, harmonic, theta, frequency)[source]

Modulated MARY reference signal.

Used for brute force modulated MARY simulations.

Parameters:

  • timeconstant (np.ndarray) – The modulation time constant (s).

  • harmonic (float) – The harmonic of the modulation.

  • theta (float) – The modulation phase (rad).

  • frequency (float) – The modulation frequency (Hz).

Returns:

The modulated reference signal.

Return type:

np.ndarray

radicalpy.utils.rodrigues_rotation(v, k, theta)[source]

Rotate vector(s) v about axis k by angle theta (Rodrigues’ formula).

Supports arrays of stacked row vectors v with shape (N, 3) or column-major layout (3, N). The axis k is normalised internally.

Parameters:

  • v (np.ndarray) – Vectors to rotate, shape (N, 3) or (3, N).

  • k (array-like) – Rotation axis (length-3).

  • theta (float) – Rotation angle (rad).

Returns:

Rotated vectors with the same shape as v.

Return type:

np.ndarray

radicalpy.utils.rotate_axes(A, x, y, z)[source]

Apply rotation matrix A to a basis (x, y, z).

Parameters:

  • A (np.ndarray) – 3×3 rotation matrix.

  • x (np.ndarray) – Input basis unit vectors (length-3).

  • y (np.ndarray) – Input basis unit vectors (length-3).

  • z (np.ndarray) – Input basis unit vectors (length-3).

Returns:

Rotated basis vectors (x', y', z').

Return type:

(np.ndarray, np.ndarray, np.ndarray)

radicalpy.utils.s_t0_omega(deltag, B0, hfc_star, onuc_all)[source]

Compute the two radical pair frequencies ω+ and ω- (in rad/s) for S-T0 mixing.

Parameters:

  • deltag (float) – Difference in g-factors of the two radicals.

  • b0 (float) – External magnetic field strength (Tesla).

  • hfc_star (float) – Hyperfine coupling constant (rad/s) of the nucleus of interest.

  • onuc_all (np.ndarray) – Array of total hyperfine contributions from all other nuclei (rad/s).

Returns:

The ω+ frequency (rad/s). omega_minus (float): The ω- frequency (rad/s).

Return type:

omega_plus (float)

radicalpy.utils.smiles_to_3d(smiles, add_h=True, opt='mmff')[source]

Create a 3D-embedded RDKit Mol from SMILES (ETKDG + optional MMFF/UFF).

radicalpy.utils.spectral_density(omega, tau_c)[source]

Frequency at which the motion of the particle exists.

Parameters:

  • omega (float) – The Larmor frequency of the electron (1/s).

  • tau_c (float) – The rotational correlation time (s).

Returns:

Spectral density frequency (1/s).

Return type:

float

radicalpy.utils.spherical_average(product_yield, theta, phi)[source]

Spherical average of anisotropic product yields.

Parameters:

  • product_yield (np.ndarray) – The anisotropic product yields.

  • theta (np.ndarray) – The angles theta by which the anisotropic product yields were calculated.

  • phi (np.ndarray) – The angles phi by which the anisotropic product yields were calculated.

Returns:

The spherical average of the anisotropic product

yields.

Return type:

float

radicalpy.utils.spherical_to_cartesian(theta, phi)[source]

Spherical coordinates to Cartesian coordinates.

Parameters:

  • theta (float or np.ndarray) – The polar angle(s).

  • phi (float or np.ndarray) – The azimuthal angle(s).

Returns:

The Cartesian coordinates.

Return type:

np.ndarray

radicalpy.utils.spin_character(state, character)[source]

Compute the “character” (overlap/content) of a spin state with respect to a chosen subspace, e.g. singlet, triplet, quintet, …

The subspace is specified by a list of kets (column vectors). The function builds the projector onto the span of those kets and evaluates the expectation value of that projector in the given state.

The state may be given either as a pure state vector |ψ⟩ or as a density matrix ρ:

  • If state is a 1-D array, the function returns \(⟨ψ|P|ψ⟩\), i.e. the probability that |ψ⟩ lies in the specified subspace.

  • If state is a 2-D array, the function returns \(\mathrm{Tr}(P ρ)\), i.e. the population of the subspace in the mixed state ρ.

Parameters:

  • state (ndarray) – Either a state vector of shape (N,) or a density

  • `` (matrix of shape)

  • character. (the kets in)

  • character (list of ndarray) – List of basis kets (each of shape (N,))

  • measured. (that span the subspace whose character is to be)

Returns:

The subspace weight / character, in the range [0, 1] for

normalised inputs.

Return type:

float

radicalpy.utils.vector_to_matrix(v, N)[source]

Convert a column vector into a matrix.

radicalpy.utils.von_neumann_entropy(rho)[source]

Compute the von Neumann entropy of a quantum state.

The von Neumann entropy is defined as \(S(\rho) = -\mathrm{Tr}\,(\rho\,\log \rho)\), which, in the eigenbasis of \(\rho\), reduces to \(S = -\sum_i \lambda_i \log \lambda_i\), where \(\{\lambda_i\}\) are the eigenvalues of \(\rho\). This implementation uses the natural logarithm, so the entropy is returned in nats.

Parameters:

  • rho (ndarray of shape (N, N)) – Density matrix.

  • states (For physical)

  • Hermitian (rho should be)

:param : :param positive semidefinite: :type positive semidefinite: Tr(rho)=1 :param and trace-normalised: :type and trace-normalised: Tr(rho)=1 :param the function does not explicitly enforce these constraints.:

Returns:

The von Neumann entropy \(S(\rho)\) in nats.

Return type:

float

radicalpy.utils.write_orca_from_pdb(pdb_path, out_dir, *, title='Tryptophan radical cation', charge=1, multiplicity=2, opt_filename=None, opt_method='RHF', opt_basis='def2-TZVP', opt_flags=('NormalPrint', 'NormalSCF'), scf_maxiter=125, scf_cnvdiis=True, scf_cnvsoscf=True, output_print_basis=2, output_print_mos=1, xyz_sidecar=None, epr_filename=None, epr_method='B3LYP', epr_basis='EPR-II', epr_autoaux=True, epr_gtensor=True, epr_nuclei_elements=('H', 'N', 'C', 'O'), epr_props=('AISO', 'ADIP'), pal_nprocs=None)[source]

Create two ORCA inputs that match your attached templates while retaining customisation knobs (methods, basis, SCF/output options, EPRNMR content).

Returns a list with paths to (opt_input, epr_input, xyz_file).

Return type:

list[Path]

radicalpy.utils.write_pdb(mol, path)[source]

Write the current conformer to a PDB file.

Parameters:

  • mol (rdkit.Chem.Mol) – Molecule to serialise.

  • path (str | Path) – Destination .pdb filepath.

Returns:

Prints a confirmation line with the saved path.

radicalpy.utils.write_sdf(mol, path)[source]

Write the molecule’s current conformer to an SDF file (one record).

Saves the active conformer of an RDKit Mol as a single SDF record, preserving topology (bonds, charges, stereochemistry) and 3D coordinates.

Parameters:

  • mol – RDKit Chem.Mol instance containing the conformer to export.

  • path – Destination file path for the sdf file.

Notes

  • Creates the file or overwrites it if it already exists.

  • Prints a short confirmation message on success.

radicalpy.utils.write_xyz(mol, path)[source]

Write the current conformer to an XYZ file (single frame).

Parameters:

  • mol (rdkit.Chem.Mol) – Molecule to serialise.

  • path (str | Path) – Destination .xyz filepath.

Notes

Writes a minimal XYZ with a generated comment line.

radicalpy.utils.yield_anisotropy(product_yield, theta, phi)[source]

Calculate the yield anisotropy.

Parameters:

  • product_yield (np.ndarray) – The anisotropic product yields.

  • theta (np.ndarray) – The angles theta by which the anisotropic product yields were calculated.

  • phi (np.ndarray) – The angles phi by which the anisotropic product yields were calculated.

Returns:

  • delta_phi (float): Maximum yield - minimum yield.

  • gamma (float): delta_phi / spherical average.

Return type:

(float, float)

Module contents

Radicalpy package root.