The Spinach kernel in Python

Spynach derives from Spinach 2.13, the spin dynamics package of Ilya Kuprov and co-workers, released under the MIT licence. The translation keeps the Spinach attribution and notice in LICENSE, and keeps Spinach's function names, parameter structures and tolerances.

The deliberate differences are small: spin indices are 0-based, product and sum operators are told apart by the type of the first argument, spin system objects are immutable, and time series are 1-D arrays. The guide Coming from Spinach lists, feature by feature, what is the same, what differs and what is not yet ported.

The translation from MATLAB to Python was carried out largely by an AI coding agent working under the direction of Alexej Jerschow, who set the goals and reviewed the results; correctness rests on the test suite against Spinach reference data.

• An oracle harness. MATLAB dump scripts write the intermediate results of a Spinach calculation to HDF5 fixtures; parity tests load them and compare the Python results at each stage (oracle/README.md).
• Stage by stage. Spin system, basis, operators, Hamiltonian, propagation, relaxation, contexts, pulse sequences: each lands with its fixtures, tests and documentation before the next starts (roadmap).
• Two formalisms behind one basis abstraction. Dense Zeeman-Hilbert space for small systems, the spherical-tensor Liouville space with state-space restriction (IK-0, IK-1, IK-2) for large liquid-state systems.
• Nothing silently ignored. Any Spinach option that is not yet ported raises an error rather than being dropped.
• Fast where it matters. The sparse products of the propagator, the Taylor stepping and the rounding run as threaded Numba kernels (the fast extra) with the same numbers as the SciPy code, and the hot loops go through one backend adapter so that a GPU array library can be plugged in; that step is deferred until a target device is chosen.
• Executable documentation. Every guide's code blocks run in the test suite and are rendered into notebooks with their outputs.

After installing: two protons at 14.1 T, pulse-acquire, exponential apodisation, two doublets.

import numpy as np
import spynach

spin_system = spynach.create(
    {"isotopes": ["1H", "1H"], "magnet": 14.1},
    {"zeeman": {"scalar": [1.00, 1.10]}, "coupling": {"scalar": [[0, 12.0], [0, 0]]}},
)
spin_system = spynach.basis(spin_system, {"formalism": "zeeman-hilb", "approximation": "none"})
pars = {"spins": ["1H"], "offset": 630, "sweep": 200, "npoints": 1024,
        "rho0": spynach.state(spin_system, "L+", "1H"),
        "coil": spynach.state(spin_system, "L+", "1H")}
fid = spynach.liquid(spin_system, spynach.acquire, pars, "nmr")
fid = spynach.apodisation(fid, [("exp", 6)])
spectrum = np.fft.fftshift(np.fft.fft(fid, 4096))

The parity suites compare Spynach with Spinach 2.13 fixtures at every stage of a calculation. Of the 540 tests, 323 are fixture comparisons over 45 fixture cases: Liouville-space setup (132), Liouville-space propagation (54), Hilbert-space pulse-acquire (64), anisotropic interactions (25), relaxation (20), scalar relaxation (9), magic angle spinning (16), pulse sequences (3). Comparisons are relative with zero absolute tolerance; the largest measured difference is 8e-11 for quantities compared at 1e-9 (time evolution and FIDs) and 9e-15 for everything compared at 1e-10 or tighter (operators, states, Hamiltonians, propagators). Every case, tolerance and measured difference is listed in docs/parity.md (regenerated from the suite by oracle/parity_summary.py).

Full-size calculations are run on both sides under a profiler, on the same machine, and summarised in an HTML report with the stage timings, the normwise difference of every saved array, the figures and the profiler tops. The timings are wall clock for each implementation as shipped: Spinach with its parallel pool and multithreaded BLAS, Spynach as a single process with multithreaded BLAS inside dense calls. The table is generated from the reports' summary.json files.

Wall clock of the whole case on each side, Spynach in its fastest configuration; agreement is the largest relative difference over the arrays both sides saved (FIDs, spectra, trajectories).

Apple silicon Mac, 12 cores; MATLAB R2026b with an 11-worker pool; Python 3.12 with NumPy and SciPy on Accelerate. Spynach is timed with its compiled kernels (the fast extra), which thread the sparse products; the methylaziridine NOESY spends nine tenths of its time in them (26 squarings of a 38 628-dimensional sparse matrix), while the strychnine and ZULF cases are dominated by dense matrix products at BLAS speed. The protocol, including how to add a case, is in comparison/README.md.

The guides run reduced versions of Spinach examples in seconds. The scripts in examples/ are the same calculations at the Spinach size, for side-by-side comparison with MATLAB: each writes its results next to itself and can report the normwise error against a .mat file saved by the Spinach script. Python code contributed by users lives in user_py/ as notebooks saved with their outputs.

The four scripts are the cases of the timed comparisons above. Contributions to user_py/ keep their author's name in the file and in its README. The guides and notebooks are listed on the documentation page.