From a Hamiltonian to a calculation you can trust
Formal spin dynamics is only half the problem. The other half is numerical: how large is the Hilbert space, which representation should we propagate, how do we avoid constructing impossible matrices, and how do we know that a stochastic or time-discretized result is converged?
Estimate Hilbert-, operator- and Liouville-space scaling before choosing a numerical method.
Understand when sparse/Krylov or state-vector methods avoid impossible dense representations.
Design convergence tests for timestep, stochastic trace samples, orientations and molecular ensembles.
Representation
Count the Hilbert space before choosing an algorithm
The numerical representation is part of the physical modelling strategy
What it is: The vector space containing all pure spin states. Its dimension is the product of the dimensions of the individual spins.
What it changes: It sets the size of state vectors and the fundamental cost of exact propagation.
What you observe: Not a laboratory observable; it determines whether a proposed simulation is computationally feasible.
What it is: The density matrix stores populations and coherences for ensembles or mixed states; vectorizing it creates a Liouville-space object of dimension \(D^2\).
What it changes: It makes relaxation and ensemble dynamics natural to express but increases memory and propagation cost dramatically.
What you observe: It is the numerical object from which populations, coherences and expectation values are calculated.
For independent spins \(I_k\), the Hilbert-space dimension is
For \(N\) spin-\(\tfrac12\) particles this becomes \(D=2^N\). A state vector contains \(D\) complex amplitudes. A dense operator or density matrix contains \(D^2\) complex numbers. A dense Liouville-space superoperator would contain \(D^4\).
Interactive
See when dense representations become impossible
Spin-space scaling explorer
At this size a state-vector method is still modest, while a dense operator is already expensive and a dense Liouvillian is completely impractical.
Memory estimates assume 16 bytes per complex double. The Monte-Carlo number is \(\sigma/\sqrt M\): the standard error relative to the single-sample standard deviation, not a guaranteed relative error in the final observable.
Propagation strategies
Do not build a matrix just because the equation contains one
Numerical efficiency comes from applying the generator without representing every zero explicitly
What it is: Most spin Hamiltonians contain local or few-spin terms, so many matrix elements in a product basis are exactly zero.
What it changes: Sparse storage and matrix–vector products can reduce memory and work dramatically compared with dense linear algebra.
What you observe: Much larger spin systems become tractable without changing the underlying physical Hamiltonian.
What it is: A method that approximates \(e^{-iHt}|\psi\rangle\) inside a low-dimensional subspace generated by repeated actions of \(H\) on the current state.
What it changes: It avoids constructing the full matrix exponential and concentrates effort on the part of Hilbert space explored during that step.
What you observe: Controlled propagation error with far lower memory cost for large sparse systems.
What it is: A numerical accuracy parameter controlling how closely the discrete propagation follows the continuous equation of motion.
What it changes: Too coarse a step can distort phase, populations or decay even when the physical model is correct.
What you observe: A result that changes when \(\Delta t\) is halved is a numerical result that is not yet converged.
Simple and accurate for small systems. Becomes memory- and compute-limited quickly as \(D\) grows.
Exploit the action of \(H\) or a sparse generator on a vector without constructing a full matrix exponential.
Propagate \(D\)-component vectors rather than \(D^2\)-component density matrices when an appropriate stochastic or pure-state formulation exists.
Useful when the density matrix or superoperator structure is essential, but it requires more aggressive sparsity or symmetry exploitation.
Trace sampling
Stochastic trace estimation trades memory for variance
Random sampling here is a numerical approximation, not necessarily physical noise
What it is: A Monte-Carlo estimator that approximates a high-dimensional trace by averaging expectation values over random vectors.
What it changes: It replaces an exact sum over a huge basis by a controllable statistical error that decreases roughly as \(M^{-1/2}\).
What you observe: Run-to-run sampling fluctuations and convergence of the final observable as the number of samples increases.
What it is: Random trajectories used to represent environmental noise or a stochastic Schrödinger equation describe actual model dynamics rather than merely estimating a trace.
What it changes: They change the time evolution being modelled, not just how efficiently a deterministic quantity is evaluated.
What you observe: Noise-induced relaxation, dephasing or distribution of trajectories after ensemble averaging.
If random normalized states satisfy \(\mathbb E[\lvert r\rangle\langle r\rvert]=\mathbf 1/D\), then a trace can be estimated as
The standard error decreases asymptotically as \(M^{-1/2}\), but the prefactor is observable-dependent. Twelve samples do not mean “12-fold accuracy”, and they do not guarantee a specific percentage error.
Convergence
Converge the observable, not just the propagator
A reliable numerical result should be tested against the knobs that can change it:
Plotting a result as a function of \(M\), \(\Delta t\), orientation count or trajectory number is much more informative than reporting one calculation and assuming it is converged.
MolSpin
Use the software as a laboratory, not a black box
MolSpin is designed around general spin systems and multiple propagation strategies. The important skill is not memorizing an input file; it is understanding which Hamiltonian, state, relaxation model, numerical method and observable the input represents.
For released syntax and examples, use the Software page and the public MolSpin documentation rather than copying teaching pseudocode from this lecture.
Selected reading
Methods behind the calculations
Spin Dynamics of Radical Pairs Using the Stochastic Schrödinger Equation in MolSpin
State-vector propagation and stochastic methods for large radical-pair spin systems.
J. Chem. Theory Comput. (2024) →Modeling spin relaxation in complex radical systems using MolSpin
Relaxation theory and practical density-matrix spin dynamics.
J. Comput. Chem. (2023) →Key external literature
Where to read next
These are deliberately selected from outside my own work: foundational papers or reviews that are especially useful for this topic.
EXPOKIT: A Software Package for Computing Matrix Exponentials
R. B. Sidje · ACM Transactions on Mathematical Software (1998). A classic reference for matrix-free Krylov evaluation of exponential propagators.
Open DOI →A Stochastic Estimator of the Trace of the Influence Matrix for Laplacian Smoothing Splines
M. F. Hutchinson · Communications in Statistics—Simulation and Computation (1989). The historical source of the random-vector trace estimator now widely known as Hutchinson estimation.
Open DOI →Spinach – A software library for simulation of spin dynamics in large spin systems
H. J. Hogben et al. · Journal of Magnetic Resonance (2011). A useful comparison point for sparse and restricted-state-space strategies in large spin simulations.
Open DOI →