Large-scale magnetic resonance simulations: A tutorial

Computational modeling is becoming an essential tool in magnetic resonance to design and optimize experiments, test the performance of theoretical models, and interpret experimental data. Recent theoretical research and software development made possible simulations of large spin systems, for example, proteins with thousands of spins, in reasonable time. In the last few years, the Fokker – Planck formalism also re-emerged due to its ability to handle spatial dynamics. The purpose of this tutorial is to describe advantages and disadvantages of the most common formalisms, the latest developments and strategies to improve the computational efficiency, and to guide users in the setting up of a simulation using the Spinach software.


| INTRODUCTION
A simulation software may be seen as a tool complementary to the spectrometer that grants a better understanding of the dynamical properties of the system, as well as for the optimization of experiments. Although the theories describing magnetic resonance experiments have been available since the 1960s, computations were inhibited by the lack of computer power and the exponential complexity scaling. Simulations involving more than 10 spins were not feasible. Calculations of this scale now take just a few seconds due to improved algorithms, an increase in computing power, and the simplicity of modern programming languages. The choice of the theoretical approach depends on the field of application. For example, the Bloch-Torrey equations 1,2 can be used for simulations of relaxation, diffusion, and flow, on isolated spin-1/2 nuclei. The density operator formalism, reliant on the Liouville-von Neumann equation 3 , can used in simulations of experiments on both single systems and ensembles 4,5 ; and allows the treatment of kinetics and sophisticated relaxation models, for example, the Redfield relaxation. 6 Models based on the Fokker-Planck formalism 7 are used in problems where spin dynamics coexists with classical dynamics (e.g., diffusion, flow, and sample rotation). Another time domain approach is the product operator formalism, which is used for calculating the evolution under free precession, chemical shift, and radio frequency (RF) pulses of the density operator (described in Section 3) in systems of weakly coupled spins. 8,9 Other methods use analytical descriptions, [10][11][12] but it is always convenient to use numerical simulations in the case of large spin systems and powders with orientational effects or spatial dynamics. The Spinach software 13 supports different formalisms in a single package, and due to the implementation of tensor structured objects 14,15 and polynomially scaling algorithms, 16,17 it enables quantum mechanical simulations previously considered unfeasible.

| CLASSICAL VECTOR MODEL BASED ON BLOCH EQUATIONS
The dynamics of isolated spins can be described in terms of the motion of a classical magnetization vector using the Bloch equations. 1,[18][19][20][21] Considering an ensemble of spins with quantum number equal to ½, in the presence of an external magnetic field, B 0 , along the Z direction; each spin occupies one of two distinct energy levels, aligned either parallel (the low energy, α state) or antiparallel (the high energy, β state) to the magnetic field. At the thermal equilibrium, there is slight excess of spins in the low-energy state that leads to net magnetization. Bloch equations connect the applied magnetic field vector, B = (B x , B y , B z ), to the bulk nuclear magnetization, M = (M x , M y , M z ). In the laboratory frame of reference, the equation of motion describing the time dynamics of the magnetization vector is: where γ is the gyromagnetic ratio, R is the relaxation matrix, T 1 and T 2 are the longitudinal and transverse relaxation times, M 0 is the magnetization vector at the thermal equilibrium, and B(t) = B 0 +B r. f. is the external magnetic field vector composed by a static part B 0 and an RF field B r. f. (t) that is time dependent. The B r. f. (t) is applied in the form of a linearly oscillating field with amplitude B 1 that can be decomposed into two counterrotating components with frequency ±ω r. f . Only the component rotating in the same sense as the spins is retained: B r:f: t ð Þ = 2B 1 cos ω r:f t + φ ð Þ e x + sin ω r: where φ is the phase acquired by the RF field and e x and e y are unit vectors in the x and y directions. The time dependence of Equation 2 can be eliminated through a transformation from the laboratory frame to a new frame, called rotating frame, where most of the simulations are performed (described in Section 3.4). In the rotating frame, the three components of the effective magnetic field vector correspond to: B r x = B 1 cosφ,B r y = B 1 sinφ, B r z = B 0 + ω r: where the superscript r indicates rotating frame, Ω = ω 0 − ω r. f is the effective precession frequency in the rotating frame, ω 0 = − γB 0 is the Larmor frequency in the laboratory frame. Notice that the time dependence has disappeared from Equation 3, but the x and y components might reacquire time dependence for pulses with variable phase and amplitude. In the rotating frame, the equation of motion can be written as: The free induction decay (FID) can then be written in terms of complex transverse magnetization according to: If pulse field gradients are present, the Larmor frequency of the spins becomes position dependent, and the field vector can be written as: where r is the displacement vector and G(t) = (G x (t), G y (t), G z (t)) is the time-dependent gradient. The effect of molecular self-diffusion can be included in Equation 1 in a generalized model proposed by Torrey 2 in 1956. The complex transverse magnetization in Equation 5 becomes a function of both of r and t; the resulting Bloch-Torrey equation can be written in terms of transport of transverse magnetization (assuming G x and G y equal to zero), in the rotating frame: where D is the diffusion tensor and the symbol r represents the gradient. Further generalizations were made by Stejskal and Tanner 22 in 1965 to incorporate anisotropic diffusion, restricted diffusion and flow; and by Henkelman et al. 23 in 1993 to describe chemical exchange. Bloch-Torrey equations can be used to describe magnetic resonance experiments involving isolated spins, for example, Fourier spectroscopy, T 1 measurements, and composite pulses. They are the pillar of several packages for magnetic resonance imaging (MRI) simulations 8,24 and can often be used to describe intravoxel dephasing effects due to RF pulses, field inhomogeneities, spin echoes, and stimulated echo sequences. Bloch-Torrey equations cannot be used to simulate pulse sequences that involve high-spin orders including multiple quantum coherence and coherence transfer processes. In Equations 1 and 7, there are no terms to describe the coupling between spins, these are included in the density operator formalism described in Section 3. In the context of MRI, the Bloch-Torrey equations fail to describe quantum mechanical effects (e.g., J-coupling evolutions) coexisting with spatial dynamics (e.g., diffusion, flow, and sample rotation). These can be simulated using the Fokker-Plank formalism described in Section 4.

| Density operator
The state of a spin system is represented by a state function or wave function. In the Dirac notation, the kth state function can be written as ket, |ki, and its transpose as a bra, hk|. A state function can be expanded in terms of orthonormal basis sets, {|φ k i,k = 1,2,. …,n}. The definition is given in Section A.2 in the supporting information: where c k are potentially time-dependent complex numbers, called expansion coefficients, and H is the overall dimension of vector space. A vector space where a norm and a scalar product (Equations A.1 and A.2 in the supporting information) are defined is called the Hilbert space. For a system of N spins, the dimensionality of the Hilbert space is: where S is the quantum number of the spin n. The state of an ensemble of spin systems is defined by the density operator,ρ , operators are specified with a single hat,Á The density operator is given by the direct product (described in Section A.3 in the supporting information) of a wave function and its conjugate transpose: The physical meaning of the density operator can be understood from the expression of its matrix elementsthe diagonal elements represent the populations of the states, and the off-diagonal elements represent the correlation coefficients. 5 In Spinach, the density operator defining the state of a single spin or a spin system is built using irreducible spherical tensor operators 13 (described in Section 3.3.2) and can be specified by calling the state.m function. For example, an initial stateρ 0 corresponding to longitudinal magnetization, represented by the operatorL Z , and a detection stateρ coil , corresponding to transverse magnetization, represented by the operatorL + , can be specified as: where the flag all represents the sum of a specified state on all the spins in the system. In the case of a spin system in the thermal equilibrium, where the energy levels are populated according to the Boltzmann law, the density operator at a given temperature T in Kelvin is given by: where k B is the Boltzmann constant, the symbol Tr is the trace of the matrix (corresponding to the sum of the diagonal elements of the matrix), andĤ is the spin Hamiltonian (defined in Section 3.3). A simulations where the initial condition is the thermal equilibrium can be set by calling the equilibrium.m function: where the spin_system object contains the spin system information (described in Section 3.2.1), H 0 is the isotropic part of the spin Hamiltonian, Q is the irreducible component of the anisotropic part (described in Section 3.3), and α, β, and γ are Euler angles defining the orientation of the system in specified powder grid. In the case of liquid state simulations, where the anisotropic part of the Hamiltonian is averaged out, the input parameters Q and [α β γ] do not need to be provided.

| Spin interactions
Most simulations require the couplings between spins to be considered, as their time dependence can provide valuable information about the system. Physical systems are composed of many interactions; but from an algebraic perspective, spin interactions can only be of three types: linear, bilinear, and quadratic. Linear interactions arise from the interaction between spins and the external magnetic field B 0 or a RF field (Equation 2); these are called Zeeman interactions. Bilinear interactions arise from homonuclear and heteronuclear couplings between spins, for example, J-coupling, hyperfine, dipolar, and exchange couplings. Finally, quadratic interactions, for example the nuclear quadrupolar interaction (NQI), that arise from the coupling between a spin and "itself". The total spin Hamiltonian can be built summing together individual interaction Hamiltonians calculated from the scalar product of spin vectors and an interaction matrix, for example, as: whereĤ z is the Zeeman Hamiltonian,Ĥ J is the Jcoupling Hamiltonian, andĤ Q is the quadrupolar Hamiltonian;Î ! andL ! denote spin operator vectors corresponding to the three Cartesian projections of the spin angular momentum operators in the space; B ! is the magnetic induction vector, and A represents an interaction tensor (a 3 × 3 matrix representing the orientation of the system in the space), that corresponds to the chemical shift tensor, hyperfine coupling tensor, and the quadrupolar coupling tensor.

| Setting up of spin interactions
Most simulation problems require the couplings between spins to be considered, as their time dependence can provide valuable information about the system. The structure array sys specifies the spin system and instrument configuration, whereas the structure inter specifies the interactions present in the system. For example, a spin system of two protons can be initialized as: Spin-field interactions require the external magnetic field strength to be supplied in units of Tesla: The full chemical shift tensors of each spin may be supplied as matrices with units of parts per million: where the indices {1,2} represents a pair of spins specified in sys.isotopes. The anisotropy of the J-coupling between two nuclear spins is usually ignored, and the corresponding orientational average is known as scalar coupling, which can be provided in Hertz: inter.coupling.scalar{1,2}=10; If Cartesian coordinates of every spin (in Angstroms) are provided, Spinach calculates interspin dipolar interactions internally: inter:coordinates = f½0:00:00:0½0:00:06:0g; If the spin system contains nuclei with spin quantum numbers greater than ½, the total spin Hamiltonian includes also quadratic terms arising from the interaction between a nuclear electric quadrupole moment and an electric field gradient (EFG). In Spinach, the quadrupole Hamiltonian in its principal axis frame (where the interaction tensor is diagonal) corresponds to: where C Q = e 2 qQ 4I 2I − 1 ð Þ h is the is the quadrupolar coupling constant in Hertz; and η is the quadrupolar tensor asymmetry parameter (calculated from the eigenvalues of the electric field gradient tensor), I spin quantum number of the quadrupolar nucleus and h is the Planck constant. The nuclear quadrupolar couplings matrix according to Equation 13 can be obtained by calling the eeqq2nqi.m function, for example as: inter.coupling.matrix{1,1}=eeqq2nqi(C Q ,η,I,eulers); where indices {1,1} indicates an interaction between a spin and itself; eulers is a vector containing three euler angles in radians, giving the orientation of the principal axis frame relative to the laboratory frame. The data structures sys and inter should be then supplied to the create.m function that builds the spin system object containing the spin system informationisotopes, spin coordinates, coupling, and so forth. The command is: This tutorial provides an example set in the supporting information, where interactions are set in accordance to the above description. Other options for the setting up of spin interactions can be found in Kuprov. 25

| Spin Hamiltonian
Molecules rotate in solution and spin in magic angle spinning (MAS) experiments (described in Section 4). Each interaction Hamiltonian contains a rotational invariant part ("isotropic") that is time independent and a part that changes with molecular orientation ("anisotropic") that is time dependent. The rotational modulation of a spin interaction tensor using Cartesian operators corresponds to: where R is a rotation matrix and R −1 is its inverse. The rotation matrix can be parameterized in various ways (using angle-axis parameterization, quaternions parameterization, and directional cosine matrix parameterization). The result of Equation 14 corresponds to a complicated trigonometric expression. A more practical representation for the description of three-dimensional rotations is based on irreducible spherical tensor operators (ISTs), indicated asT l,m , where l is the rank and m = − l, − l+1,…,l represents its 2l+1 components. Considering that the anisotropic part of all the spin interactions transforms virtually in according to a second-rank irreducible representation of the rotation group, the rotational modulation of these interactions may be treated in a general way using the IST operators. [26][27][28][29][30] A general rotationR α, β, γ ð Þwith Euler angles α, β, and γ can be expressed as: where D l ð Þ m 0 ,m α, β, γ ð Þ are elements of the Wigner rotation matrix of the order l . 31,32 The total Hamiltonian can then be expanded as a linear combination of IST operators corresponding to spin interactions according to:Ĥ whereĤ iso represents the isotropic part of the Hamiltonian, D l ð Þ km are the Wigner functions specifying the molecular orientation, andQ l ð Þ km contains rotational basis operators and represents the anisotropic part of the Hamiltonian. For a spin system of two spins, where one of them is quadrupolar, the rotational basis operators correspond to: where B indicates the static magnetic field, L and S indicate the spins, Φ m (B,L), Φ m (B,S), Φ m (L,S) and Φ m (S,S) are irreducible spherical tensor parameters for linear, bilinear, and quadratic interactions, respectively. The isotropic part of the spin Hamiltonian in Equation 16 corresponds to zero-rank terms that are invariant under rotation. The anisotropic part of the spin interactions corresponds to second-rank terms that are the only terms that have to be rotated in practical calculations. Firstrank terms in Equation 16 are very small and are often ignored. In Spinach, the spin Hamiltonian operator in the Hilbert space and superoperator in the Liouville space (described in Section 3.5), along with the rotational decomposition, can be obtained by calling the hamiltonian.m function: where H 0 represents the isotropic part of the spin Hamiltonian and Q is the irreducible basis of the anisotropic part. The spin Hamiltonian in Equation 16 can then be built with the following command: The orientation.m function builds the anisotropic part for a specified set of Euler angles [α β γ].

| The IST parameters
A total Hamiltonian in Equation 16 can be then written in terms of all the interaction tensor elements. In Equation 17, zero-rank IST parameters are given by: Second-rank IST parameters are given by: ÞÇi a yz + a zy À Á 2 , The irreducible spherical tensor parameters in Equations 18-19 can be obtained from a 3 × 3 interaction matrix by calling the mat2sphten.m function: [rank0,rank1,rank2]=mat2sphten(A); Spinach converts a 3 × 3 interaction matrix A into one zero-rank IST parameter, three first-rank IST parameters and five second-rank IST parameters according to Equations 18-19.

| The IST operators
Zero-rank ITS operators correspond toT 0,0 =Ê , whereÊ is the unit matrix. Second-rank ITS operators in Equation 17, describing linear interactions, correspond to single operators and are given by: Second-rank ITS operators, describing bilinear and quadratic interactions between spins, correspond to combinations of two-spin product operator basis given by: Ladder operators appearing in Equations 20-21 are given byL AE =L X AE iL Y , where the spin operator basisL X , L Y , andL Z that can be built from the direct product of Pauli matrices. Single-spin IST operatorT k,m with rank k and multiplicity m can be obtained calling the irr_sph_ten.m function: T=irr_sph_ten(m,k); For example, considering a system with two ½-spins, the spin operator for each spin can be obtained from: whereσ x ,σ y , andσ z represent the Pauli matrices andÊ is the unit matrix, corresponding to: In Spinach, the spin operator basis in Equation 22 can be computed by calling the operator.m function that builds an operator in the Hilbert space and a commutation superoperator in the Liouville space (described in Section 3.5). For example, spin operators corresponding to the sum ofL X andL Y on all the protons of a spin system are computed as: The pulse operators can then be built: Lx = (Hp + Hp')/2; and Ly = (Hp-Hp')/i2; Pulse operators must be supplied in any Spinach calculations, so see the example set for details.

| Laboratory and rotating frame Hamiltonians
The spin Hamiltonian given in Equation 16 contains all the couplings present in the system. Considering a spin system of two spins L and S connected by J-coupling, the laboratory frame Hamiltonian used in liquid and solid state simulations corresponds to: whereĤ z,iso is the isotropic part of the Zeeman term corresponding to: where Z L is the Zeeman interaction matrix of the spin L and Z S is the Zeeman interaction matrix of the spin S. The termĤ J,iso is the isotropic coupling term corresponding tô where J LS is the J-coupling interaction matrix. The laboratory frame Hamiltonian contains also Zeeman carrier frequencies given by the product of the isotropic freeparticle gyromagnetic ratio and the magnetic field, corresponding to ω l,0 = − γ l B 0 and ω s,0 = − γ s B 0 for the spin L and S, respectively. The Zeeman carrier frequency does not carry any relevant information and dominates the cost of the calculation. For example, in a magnetic field of 14 Tesla, the Zeeman Larmor frequency is of the order of MHz and the time step is of the order of microseconds; processes occurring on time scales of seconds would require millions of points, leading to extremely long wall-clock times. Thus, it should be removed from the Hamiltonian. The procedure to remove Zeeman terms from the spin Hamiltonian is called interaction representation transformation or rotating frame transformations 5 that yield the secular Hamiltonian (the derivation is omitted):Ĥ In rotating frame simulations, Zeeman carrier frequencies are subtracted from the coefficient in front of theL Z andŜ Z operator in the isotropic Zeeman terms: In the case of couplings between spins of different types, the IST operator in the secular part of bilinear couplings, in Equation 27, simply corresponds tô The latter is called weak coupling term. The secular Hamiltonian is used in NMR simulations performed at high-field (≥1 Tesla). Spinach also allows more complicated higher-order corrections to the rotating frame transformation. 80 In rotating frame simulations, the offset with respect to the carrier Zeeman frequency can be specified as an input parameter: parameters.offset and can be added to the Hamiltonian by calling the frqoffset.m function. In Spinach, laboratory frame or rotating frame Hamiltonians can be set automatically using the assume.m function that builds the spin Hamiltonian in accordance with Equations 24 and 27: spin_system=assume (spin_system,'assumption'); or using a kernel context.m function: answer=context (spin_system, @pulse_sequence,parameters,'assumption'); where the flag assumption can be set equal to labframe to build a laboratory frame Hamiltonian, or equal to NMR to build a secular Hamiltonian. The MATLAB command window will display the following outputs: Summary of Zeeman tensors and their anisotropies (ppm for nuclei, g-tensor for electrons) +0.000000e+00 +0.000000e+00 +1.505100e+02 Summary of the tensors and their anisotropies in Hertz: where Tr/3 is the isotropic part of the Hamiltonian and norm[rank1] and norm[rank2] represent norms of the first-and second-rank terms, respectively. Context functions are responsible for running miscellaneous housekeeping (liquid state calculation, powder averages, spatial dynamics, etc.), creating spin operators or superoperators, kinetic and relaxation superoperators, and passing them to a pulse_sequence.m function. Available context functions are liquid.m (for simulation in liquid state), powder.m (for static simulations in solid state), imaging.m (for MRI and spatiotemporal NMR simulations), and singlerot.m (for single axis magic spin angle simulations). The pulse_sequence.m function must contain the following command: answer=pulse_sequence(spin_system, parameters, H,R,K); A spin Hamiltonian operator in the Hilbert space or a Liouvillian superoperator in the Liouville space (described in Section 3.5) can then be built with the following syntax: where H is the Hamiltonian, R is the relaxation superoperator, and K is the kinetic super operator that can be obtained by calling the functions relaxation.m and kinetics.m, respectively.

| Equations of motion
The density operator formalism is derived from the solution of the time-dependent Schrödinger equation. 5 The equation of the motion for the density operator is called the Liouville-von Neumann equation: whereĤ t ð Þ is a time-dependent Hamiltonian and |ψ(t)i is a time-dependent wave function. The Liouville-von Neumann equation is the pillar of the most magnetic resonance simulation software 9-12,33 and mostly used for large-scale liquid and solid-state simulations.

| Simulations in the Hilbert space
The Hilbert space can be used for frequency domain simulations (not described in this tutorial), for determining energy level diagrams and transition probabilities between energy levels, but they do not support time domain process. In the Hilbert space, Equation 29 has a solution corresponding to: whereρ 0 ð Þ is the initial density matrix,Û t ð Þ is the propagator from the time zero to the time t, andT is the Dyson time ordering operator that evaluates the exponential functions and reorders the Hamiltonians in a way that preserves the descending time sequence to ensure that the Hamiltonian commutes with itself at different time steps. 5 In the rotating frame, the Hamiltonian is made time independent within finite segments of time. The density matrix propagated in time intervals under a timeindependent Hamiltonian corresponds tô where Δt is the time step and e + iĤΔt is a matrix exponential called propagator. For example, considering a spin system with two ½-spins, the matrix dimension for operators, propagators, and states in the Hilbert space is equal to 4 by 4, Equation 31 has the following structure: involves a double-sided product that is inconvenient in the case of multiple nested products, so a single-sided simplification can be made by moving into Liouville space.

| Simulations in the Liouville space
Most of the simulations are performed in the Liouville space, because it supports sophisticated relaxation theories (e.g., the Redfield theory) and kinetics. The Liouville space is the adjoint representation of the Hilbert space. Adjoint operators (called superoperators) are denoted aŝ whereρ t ð Þ j i is the state vector andĤ is Hamiltonian commutation superoperator that corresponds to a commutator with the Hamiltonian:Ĥρ whereĤ is the spin Hamiltonian operator andĤ T is its transpose,Ê is the unit matrix of the same size of the spin Hamiltonian. In the Liouville space, the elements of the density matrix are arranged in a column vector, called state vector. In Spinach, Hilbert space operators can be converted into Liouville space superoperators by calling the hilb2liouv.m function: where the flag comm indicates a commutation superoperator. In the Liouville space, the evolution can be expressed as: The matrix dimension for propagators in the case of a spin system with two ½-spins is equal to 16 by 16, and the dimension of the state vector is 16 by 1, Equation 34 corresponds to: Although the matrix dimension is larger in the Liouville space, it allows polynomial scaling with the size of the system (described in Section 5). An example of setting up of simulation in Liouville space is shown below. Figure 1 depicts a simulation of a 1D NMR spectrum of a glycine system where the relaxation behavior is controlled by both the chemical shift anisotropy (CSA) and dipole-dipole (DD) coupling.
The DD-CSA interference causes the double members in the spectra to relax at different rates. Simulation parameters and the code (explained step by step) used to obtain Figure 1 are given in Example 1 in Section B.1 of the supporting information.

| Setting up of the Liouville and Hilbert space
Spinach supports both the Hilbert and the Liouville space simulation formalisms. The formalism can be specified using bas.formalism structure. Simulations in the Hilbert space use Pauli matrices as basis sets and can be specified with the command: bas.formalism='zeeman-hilb'; Simulations in the Liouville space use IST operators as basis and can be specified with:

bas.formalism='sphten-liouv';
Information about the formalism must be supplied to the basis.m function using the following command: spin_system=basis (spin_system,bas); Spinach supports complete and incomplete basis sets, as described in Section 5.

| Setting up of the propagation method
In Spinach, long-range propagations under a static Hamiltonian and Liouvillian in accordance to Equations 31 and 34, respectively, can be performed by calling the evolution.m function: rho_final=evolution (spin_system,L,coil,rho_zero,… time_step,nsteps,output,destination); where the rho_final corresponds to the final state vector (or density matrix), depending on the options specified in the function call, L is a spin Hamiltonian or a Liouvillian, the coil (optional argument) represents the detection state, rho_zero is the initial state vector or density matrix to be propagated, nsteps represents the number of simulation steps characterized by duration given in time_step, and the output specifies the form of rho_final that can be a single state vector (or density matrix) or a simulation trajectory, the destination flag is an optional to be used for destination screening (described in Section 5.1.7). This function should be used for trajectory calculations, detection periods of time domain experiments, and computation of the time integrals of observables. In Spinach, the propagation can also be performed using the Taylor series 13 : Equation 35 involves only matrix-vector multiplications and calculates the action by a matrix exponential on a vector without computing the matrix exponential. The scaling of the Taylor series is O (N 2 ), where N is the matrix dimension. This approach is particularly well suited for propagation under single hard pulses and timedependent Liouvillians or Hamiltonians. The time propagation with the Taylor series can be performed by calling the step.m function with the following command: where L can be a pulse operator (described in Section 3.3.2), a spin Hamiltonian operator or a Liouvillian and time_step is the time step to take. Considering that the Taylor series is polynomial, its convergence to the machine precision will depend on the product of the norm of the spin Hamiltonian operator F I G U R E 1 (a) Simulation of a 1D nuclear magnetic resonance spectrum of the glycine system ( 1 H- 15 N bond), showing the different intensity of the two lines of the 15 N doublet due to dipole-dipole and chemical shift anisotropy interference. (b) Simulation of the 1D nuclear magnetic resonance spectrum excluding the dipole-dipole contribution (or the Liouvillian) and the time step. When the norm is smaller than one, the Taylor series converges in less than 20 steps, and when the norm is higher than one, the propagation is performed using scaling and squaring techniques 37 that require matrix exponentialization. The matrix exponential is computed in n squaring steps that involve n matrix-matrix multiplications according to: The scaling of the matrix exponentialization is O (N 3 ), where N is the matrix dimension. The propagation of a spin system with more than 200 spins becomes problematic. An alternative is the Krylov propagation that computes the product between the matrix exponential propagator and the density matrix directly, thus avoiding matrix exponentialization. 38 In Spinach, the size of the Liouvillian is evaluated internally and the switch to Krylov propagation is made automatically in the case of a very large Liouvillians, like those found in protein spin system. The Krylov propagation requires less memory than the scaling and squaring technique but is slower and hence should only be used when exponential of the Liouvillian does not fit into memory, but the Liouvillian itself does. If memory is available, the switch to Kyrlov propagation can be disabled with the command: The krylov propagation can be also used within a pulse sequence by calling the krylov.m function: answer=krylov (spin_system,L,coil,rho_zero,… time_step,nsteps,output); This function uses the same input parameters as evolution.m.

| Expectation value and trajectory analysis
In our simulations, we are often interested in the evolution of on observable property that can be extracted from the density operator. The term observable corresponds to properties of an individual molecule rather than ensembles. Operators corresponding to observables are indicated with the symbolÔ. The average outcome of many repeated experiments of an observable property of a quantum system is called the expectation value and is indicated with the angular brackets hÁi: The expectation value can be found by evaluating the trace of the product of an observable operator and the density operator. In practice, the expectation value related to an observable operator can be computed as: Hx=real (Hx'*trajectory); where Hx represents an observable spin operator and trajectory is a simulation trajectory obtained with the evolution.m function by setting the option with the flag trajectory. Considering a spin system of two coupled spins, 1 H and 13 C, an example of evolution of observable single spin operators over a simulation trajectory is shown in Figure 2.
The simulation starts with transverse magnetization on the proton, then it is transferred to the carbon during the spin lock. Simulation parameters and the code (explained step by step) used to obtain Figure 2 are given in Example 2 in Section B.2 of the supporting information. In Spinach, a trajectory analysis can be also performed by calling the trajan.m function: trajan (spin_system,trajectory,property,time_axis); where property corresponds to an observable property of the system, such as a correlation order, coherence order, or total state space population.

| FOKKER-PLANK FORMALISM
In the Fokker-Planck picture, the state of the system is not represented by the density matrix ρ(x,t) that describes the state of a specific particle at the point x at the time t. It considers the probability p(x, ρ, t) of having a particle in the point x in the spin state ρ at the time t in a very high dimensional space (three dimensions with respect to the spatial coordinates and 4 n dimensions with respect to spin). This approach can be used in simulations where there is the simultaneous presence of quantum mechanical dynamics (i.e., spin-spin coupling, spin relaxation, and chemical reactions) and spatial dynamics (e.g., diffusion, flow, or sample rotation), for example, singlet and metabolomics imaging, 39,40 pure-shift NMR, 41 and ultrafast NMR. 42 The Fokker-Planck equation can be derived from the continuity equation that states that the probability must be conserved and the local change in probability as a function of space and time is equal to the divergence of its flux, j(x,t): where r is the gradient containing the spatial derivatives r = ∂=∂x∂=∂y∂=∂z ½ T . For a given probability profile, the flux is given by the product of the probability and the velocity field: Considering that the velocity field in space can be described by the classical equation of the motion and the velocity field in quantum mechanics can be described by the Liouville-von Neuman equation, it is possible to write separately the divergence with respect to the spatial r x and spin r ρ coordinates. The former is given by the velocity in the space v x (x, ρ, t) multiplied by the probability p(x, ρ, t), whereas the latter is given by the velocity in the spin space v ρ (x, ρ, t) multiplied by the probability: The divergence represents the sum of all the derivatives with respect to all the coordinates: The average of any property can be obtained from the integral over the entire configurational space of the product between that property and its probability. The ensemble average density matrix corresponds to: where dV ρ is the volume of the spin space. The equation of the motion can be obtained from the time derivative of the average density matrix and is obtained by differentiating Equation 42 with respect to the time and then using the expression for the derivative from Equation 40, to obtain: Equation 43 is the Fokker-Planck equation, in which L(x,t) is the Liouvillian and M(x,t) is the spatial dynamics generator obtainable from classical mechanics. Time evolution generators corresponding to infinitesimal translation on a spatial or a phase grid are represented by where v(t) is the flow velocity, ω(t) is the angular velocity, and D is the diffusion coefficient. In cases where M(x,t) corresponds to the diffusion operator, Equation 43 is known as the stochastic Liouville equation (SLE), which was used for the first time in magnetic resonance in 1970 for simulations of EPR spectra. [43][44][45][46][47] The SLE could describe the rotational diffusion at all magnetic field strengths and correlation times, going from extreme narrowing to solid limits, and has recently been applied to MAS simulations. 48,49 The probability density can be moved forward in the time by applying an exponential with respect to a generator according to where a corresponds to linear velocity in the case of flow or angular velocity in the case of phase increment, p(t) is the initial probability, and ∂/∂t corresponds to one of the derivatives shown in Equations 44-46. An operator of infinitesimal translation squared is a generator of infinitesimal diffusion, Equation 46. Its solution is equivalent to Green's function of the heat equation describing the dynamics of a particle undertaking Brownian motion starting from x 0 at the time zero t 0 to the position x; the time t can be described by the conditional probability: In Spinach, spatial dynamics generators have a matrix representation. Flow and diffusion generators (Equations 44 and 46) are represented by central finite difference matrices, and spatial coordinates are discretized on a finite grid, 50 for example, as where f ' (x) is the derivative of a function defined on a finite grid, Δx is the spacing between discretisation points, and O(Δx 2 ) is the residual error determining the accuracy of the calculation. For example, if Δx is twice finer, the grid will become twice finder and the error will be reduced by a factor equal to Δx 2 . In Spinach, finitedifference derivative matrices can be set as parameters.deriv={'period',n}; where the flag period indicates finite difference approximations with periodic boundary conditions and n is the stencil size made by a point on the spatial grid with two, four, or six neighbors and can be set equal to 3, 5, and 7. The stencil size and the number of spatial points determine the accuracy of the calculation. Accurate simulations require proper digitization of the spatial dimension for the correct treatment of diffusion and flow. All the frequencies along the spatial dimension have to be below the Nyquist condition of the grid: where L is the sample size, kH z k 2 is the 2-norm corresponding to the largest eigenvalue of the Zeeman Hamiltonian, g(t) is the gradient amplitude in the chosen dimension as a function of the time, and B 0 is the external magnetic field. The main source of high frequencies along the spatial dimension are gradient fields, so the highest frequency in the strongest gradient present in the sequence has to be correctly digitized. In Spinach, the ngridpts.m function can be used to determine the minimum number of points required to satisfy the Nyquist condition of the grid: number_of_points = ngridpts (gradient_amplitude,gradient_duration, … isotope,max_coherence_order,sample_size); Accurate simulations require a number of points six or seven times greater than the minimum number of points given by Nyquist condition, or as many as required for convergence. 14 In the case of RF and rotor phase increment generators (Equation 45) used in simulations of shaped, frequency swept pulses and MAS, Spinach uses Fourier differentiation matrices 51 on a grid with N regularly spaced points: Fourier derivatives are calculated in the Spinach's internal kernels and do not require specification as input parameter. In the Fokker-Planck formalism, the time propagation occurs in the direct product between spin space and spatial degrees of freedom: Fokker-Planck = Liouville space Spatial degree of freedom : Considering a 10-spin system, the dimension of the Liouville space is of the order of 4 10~1 0 6 , and in the case of three-dimensional spatial dynamics that requires at least 100 points along the three dimensions, the result is a matrix of the order of 10 12 . Matrices of this dimension can not be stored without the polyadic decomposition, 14  At the time of writing, referring to the Spinach release 2.5, the polyadic is supported only within the imaging.m context (described later in the text) and should be set for three-dimensional MRI simulations.

| Fokker-Planck equation in imaging simulations
When we have multiple motions with the corresponding generators, then the overall dynamics is given by the sum of all the generators of individual dynamics. For example, the unidimensional Fokker-Planck equation that describes frequency swept pulses and shaped pulses in amplitude-frequency coordinates in the presence of magnetic field gradients, diffusion, and flow is given by where L is the Liouvillian, ∂/∂φ is the generator of the RF phase increment that increments the phase with rate ω RF (t), ∂ 2 /∂z 2 and ∂/∂z are the gerarators of spatial dyanmics, D T is the translational diffusion coefficient in case of isotropic diffusion, and v(t) is the flow velocity. In Fokker-Planck simulations, state vectors ρ(φ (n) , z (k) , t) are stacked vertically in order of increasing index of the spatial coordinates grid points, followed by the increasing index of the RF phase grid points. The matrix representation of Equation 53 is given in Figure 3. Each element of the Liouvillian corresponds to where H 0 is the time-independent part of the Hamiltonian and contains the static interactions; a(t) is the time-dependent amplitude of the RF pulse, φ 0 is the initial phase of the RF pulse, φ n is the phase of the RF pulse at the point n, {S X , S Y } are Cartesian radiofrequency operators, g(t) is the time dependent gradient amplitude,z k is the spatial coordinate in the point k, γ m is the gyromagnetic ratio of the m th spin in the system, and S Z is the longitudinal spin commutation superoperator for each spin m, R is the relaxation superoperator and K is the kinetics super operator, D z and D 2 z are finite-difference derivative matrices on the spatial grid, and D φ is a Fourier derivative matrix on the phase grid. The digitalization of the RF phase and spatial dynamics on a grid makes the Liouvillian a constant matrix, and generators of time-dependent process The parameter L is the Liouvillian, Sx and Sy are the pulse operators. Simulation parameters can be provided for the irradiation frequency, a value for the amplitude, a value for the phase of the pulse, a value for the overall pulse duration, and a value for the number of steps representing the pulse discretization points: The number of points in the RF phase grid is equal to 2k + 1, where k is termed the maximum rank and usually takes a value of 2 or 3. It can be set as:

parameters.max_rank
The function shaped_pulse_af.m is used to simulate soft off-resonance and slice-selective pulses. Figure 4 shows an example of a pulse acquire experiment using a Gaussian soft pulse, applied on a system of 39 coupled protons, to select a signal at a specific irradiation frequency.
In this example, the simulation starts with longitudinal magnetization then a 90 hard or soft pulse is applied at 0 Hz before the signal acquisition. The selective pulse is simulated with the Fokker-Planck equation (Equation 53). Diffusion, flow, and pulse field gradients are not included in this calculation. Simulation parameters and the code (explained step by step) used to obtain Figure 4 are given in Examples 3a and 3b in Section B.3 of the supporting information. In the case of an MRI simulation involving slice selection pulses, the simulation must include pulse field gradients. The gradient operator is computed by calling the g2fplanck.m function, with the syntax G=g2fplanck(spin_system,parameters); where parameters provide the number of points in the spatial grid and the sample dimension. A threedimensional MRI simulation can be set for example: An initial condition corresponding to longitudinal magnetization and a detection state corresponding to transverse magnetization in each point of the grid can be set, for example, as parameters.rho0_ph={ones (prod (parameters.npts,1))}; parameters.rho0_st={state(spin_system,'Lz','1H','cheap')}; parameters.coil_ph={ones (prod (parameters.npts,1))}; where parameters correspond to flow velocities and diffusion that can be set on a spatial grid by providing a vector of values describing the evolution of the velocity in time. For example, constant flow velocity along each dimension can be set by: parameters.u=flow_volocity_z*ones (parameters.npts); parameters.v=flow_volocity_y*ones (parameters.npts); parameters.w=flow_volocity_x*ones (parameters.npts); Isotropic diffusion at each point of the spatial grid can be set with the parameter: The Liouvillian, the relaxation and the kinetics superoperators, the Fokker-Planck spatial dynamics generator (including diffusion and flow), and the gradient operators can be generated internally by calling the imaging.m context function: answer=imaging (spin_system,@pulse_sequence, … parameters); that passes the operators to a pulse_sequence.m function: answer=pulse_sequence(spin_system,parameters, H,R,K,G,F); The Liouvillian is built according to L=H+F+1i*R+1i*K; The gradient operator can then be added to the Liouvillian with the following command: L=L+parameters.gradient_amplitude*G{1}; where G{1} indicates a unidimensional gradient. Figure 5 shows an example of a three-dimensional MRI simulation where a slice selection pulse is applied to a brain phantom to select a single slice at a specified frequency.
In this example, the simulation starts with longitudinal magnetization on all the protons, and then, a 90 slice selection pulse is applied at given frequency. Finally, the observable operatorL + can be computed using the fpl2phan.m function: Simulation parameters and the code (explained step by step) used to obtain Figure 5 are given in Example 4 in Section B.4 of the supporting information. An alternative approach for computing shaped pulses is the time slicing method, 52,53 made available in Spinach by calling the shaped_pulse_xy.m function, which computes the shaped pulse in Cartesian coordinates. In this approach, shaped pulses are simulated explicitly using a timedependent Liouvillian, whereas in Fokker-Planck F I G U R E 5 Simulation of a magnetic resonance imaging slice selection pulse using the Fokker-Planck formalism. On the left: An initial phantom representing a human brain. On the right: A slice of the human brain after the slice selection pulse formalism, the time dependence is eliminated. In the time slice calculation, the total duration of the pulse is divided into n slices of duration Δt; the density matrix is then propagated in each time slice in accordance with Equations 30-31. The total Liouvillian that contains time independent and RF terms is diagonalized to yield an eigenvalue matrix D and an eigenvector matrix, V, which are used to create a propagator P n in each time instance n according to The overall propagator P is built from the product of all the propagators P n and then multiplied by the initial state vector. The time discretization is also performed in the Fokker-Plank formalism, but algorithms implemented in Spinach avoid the diagonalization in Equation 55 because it scales according to O (N 3 ) with the matrix dimension N. A frequency-swept pulse in the presence of a gradient can be simulated, for example as: where the parameter L is the Liouvillian, Sx and Sy are the pulse operators, G{1} indicates a unidimensional gradient (optional), gradient_amplitudes is a vector representing the variation of the gradient in the time, and Cx and Cy are the Cartesian amplitudes in X and Y coordinates at each time slice. Chirp pulses in the presence of a gradient are used in ultrafast NMR for spatial encoding; however, not all the full pulse is active for excitation and inversion. Simulations can show the active region in a chirp profile used to apply spatial encoding to a quantum coherence in a spin system of three spins, see Figure 6.
In this simulation, a spin echo block generates multiple quantum coherence, then the coherence +1 is selected, and the chirp pulse inverts the coherence order from +1 to −1. The profile is linear in a region of about 10 mm, but the diffusion adds a distortion. Simulation parameters and the code (explained step by step) used to obtain Figure 6 are given in Example 5 in Section B.5 of the supporting information.

| Comparison between propagation methods
In Spinach, the propagation of the Fokker-Planck equation can be performed either with the exponential propagation method or with the Taylor series 13 (Equation 35). The former approach requires the calculation of a matrix exponential propagator that involves matrix-matrix multiplications that scales by O (N 3 ) with the matrix dimension. A shaped pulse in amplitudefrequency or Cartesian coordinates can be simulated using this approach with the following setting: parameters.method='expm'; The Taylor series involves only matrix-vector multiplications that scale by O (N 2 ). A shaped pulse in amplitude-frequency or Cartesian coordinates can be simulated using this apporach with the following setting:  Table 1 shows a comparison between the exponential propagator method and Taylor series method for the simulation of a shaped pulse in frequency-amplitude using the Fokker-Plank formalism (with the shaped_pulse_af.m function) and Cartesian coordinates (with the shaped_pulse_xy.m function) using the time slicing approach. The memory requirement for the storage of the Liouvillian superoperator in computations using the Taylor series is smaller than the requirement for the storage of the matrix exponential propagator in computations using the exponential propagation method, details are given in Table S1. Table 1 shows an improvement in efficiency when calculations are performed using the Taylor method. Computations using the time slicing method are faster than calculations performed using the Fokker-Planck formalism. However, a pulse in Cartesian coordinates requires thousands of points to correctly digitize its rapidly oscillating x and y components, whereas fewer points (typically 100) are required to digitize a pulse in amplitude-frequency coordinates. 48 Therefore, calculations become faster using the Fokker-Planck formalism when the propagation is computed with the Taylor series.

| Fokker-Planck equation in MAS simulations
The Fokker-Planck formalism can be also used in the context of NMR MAS. Its advantage relies in the disappearance of the time dependence from the MAS Hamiltonian, leading to an improvement in efficiency. 49 The time dependence is also eliminated in the Floquet theory, and previous work has shown comparable performance between the two theories. [54][55][56] A lot of work has been invested in improving the efficiency of the MAS calculations through the optimization of two-and three-angle spherical grids. [57][58][59][60][61][62][63][64][65] Considering that the gamma angles in the powder grid and the sample rotation angle both correspond to a rotation about the rotor axis, the gamma angles can be averaged over the rotor phase, and only a twoangle grid is required. This is a feature of several algorithms including the γ-COMPUTE algorithm. 65,66 The Fokker-Planck equation in MAS simulations corresponds to: where L is the Liouvillian, n = 111 ½ = ffiffi ffi 3 p is the spinning axis direction, and ω MAS and φ MAS are the MAS spinning rate and phase, respectively. In the Fokker-Planck formalism, the generator of time-dependent processes becomes time independent. The derivative ∂/∂φ is the phase of the rotor incremented by a time-independent Fourier derivative (Equation 51). The rotor grid corresponds to a circular grid with 2k + 1 points, where k is termed maximum rank of the Fokker-Planck theory, usually set equal to 2 or 3. This parameter can be set as:

parameters.max_rank
The minimum number of points in the rotor grid for these simulations is given by the ratio between the largest interaction present in the system divided by the spinning rate. The MAS Liouvillian becomes time independent through the digitalization the rotor phase on a periodic two-angle spherical grid {φ j }: Wall-clock time (wct) in computing a single time slice using one single central processing unit core considering a spin system composed of 10 coupled protons (J-couplings were set equal to 10 Hz). The simulation included different terms of Equation 53, isotropic diffusion was set equal to 18 × 10 −10 m 2 s −1 , the flow velocity was set equal to 1 × 10 −2 ms −1 , the number points in the spatial grid was set equal to 100, the maximum rank of the Fokker-Planck theory was set equal to two.

Wct (in seconds) to compute the RF term in Equation 53
Wct (in seconds) to compute the RF term + gradient term in Equation 53

Wct (in seconds) to compute RF term + gradient term + hydrodynamics (flow and diffusion) in Equation 53
Shaped_pulse_af where Q km is the rotational basis operators and R and K are the relaxation and kinetics superoperators, respectively. The matrix form of Equation 57 is given below: The Wigner D-matrices in Equation 58 accomplish the following transformations: 1 D crystal α PC ,β PC , γ PC ð Þ rotates the molecule from the principal axis frame, in which the interactions are defined, into the crystallite-fixed frame, where α PC , β PC , and γ PC represent the Euler angles between the two frames. In Spinach, a single crystallite orientation can be specified as:

parameters.orientation
A set of crystallite orientations on a spherical averaging grid can be specified as: parameters.grid 2 D rotor φ j rotates the crystallite-fixed frame into the rotor-fixed frame defined by an orientation with a specified rotor phase. The rotor axis direction can be specified as:

parameters.axis=[sqrt(2/3) 0 sqrt(1/3)];
The spinning rate ω MAS is defined by: Þrepresents the transformation between the rotor-fixed frame and the laboratory frame corresponding to the rotor axis direction n, and α RL , β RL , and γ RL represent the Euler angles describing the transformation between the rotor-fixed frame and the laboratory frame where the simulation is performed.
The transformations described above are performed internally by the singlerot.m context function. Figure 7 shows an example of simulation of a crosspolarization experiment between a proton and a 15 N nucleus with and without MAS.
In this simulation, the polarization is transferred from the spin with higher γ to the spin with smaller γ through a dipolar coupling. The transfer occurs at the Hartmann-Hahn matching condition, which is fulfilled when the RF power at the nuclei is equalized, in according to γ1 H B 1, 1 H = γ15 N B 1, 15

| Fokker-Planck equation in overtone MAS simulations
Overtone simulations are usually performed using a laboratory frame Hamiltonian because there is no convenient rotating frame due to the large difference between frequencies involved in the spin evolution process (order~Hz for 1 H and order~MHz for 14 N overtone), and the spinning frequency is of the order of kHz; and makes the spin Hamiltonian time dependent. Laboratory frame simulations require two points per period of the fast oscillation present in the system. Therefore, millions of time steps have to be averaged over a threeangle spherical grid. The convergence in an overtone pattern requires a number of points in the spherical grid much higher than those required to digitize single quantum transitions with the same number of spinning sidebands. Consider a spin system composed of a nitrogen 14 N and a proton, the Fokker-Planck equation in overtone MAS simulations is: where L is the Liouvillian: where H 0 is the isotropic part of the spin Hamiltonian; Q km is the rotational basis operators; R and K are the relaxation and kinetics superoperators, respectively; α, β, and γ are the three Euler angles in the powder grid; θ is spinning axis at the magic spin angle; and ω MAS and ω N RF are the MAS spinning rate and overtone frequency, respectively, specified by: parameters.rate and parameters.overtone_frquency The terms φ MAS and φ N RF correspond to the rotor and nitrogen RF phases, respectively, specified by: The terms a H RF and a N RF are the RF power amplitudes of the proton and nitrogen specified by:

parameters.rf_pwr
This parameter requires a vector of two values where the first value corresponds to the RF power on the 14 where exp (O) is a time-ordered exponential that can be computed by slicing the RF period T into 16 slices and then multiplying the slice propagator using the propagator squaring approach. 69 The Hamiltonian in Equation 61 becomes time-independent, the time propagation for the duration of the cross-polarization period can be computed according to: whereρ 0 andρ t CP ð Þ are the initial density matrix and final density matrix after the cross-polarization period, respectively; and t CP is the duration of the cross-polarization period. This parameter can be set, for example, as:

parameters.rf_dur
Considering that only few hundred points occupied by the overtone spectrum have to be acquired, the acquisition is performed in the frequency domain detection that is much faster than the time domain detection 67,68 : The Fourier transform of the free induction decay is performed on 14 N and corresponds to: whereŜ N + is the detection state on 14 N,Ê is the unit matrix, and r is the damping rate parameter. The latter parameter is required to avoid singularities in the complex plane. The matrix-inverse-times-vector operation can be performed with an iterative solver called ILU preconditioned generalized minimal residual (GMRES). 70  The initial condition and the detection state corresponding to the proton and nitrogen magnetization at the MAS can be set, for example, as: parameters.rho0=cos(atan (sqrt (2)  The cross-polarization process has already been described above in this section, Figure 8 shows the 14 N overtone excitation at 48 kHz. Simulation parameters and the code (explained step by step) used to obtain Figure 8 are given in Example 7 in Section B.7 of the supporting information.

| EFFICIENCY OPPORTUNITIES
The computational efficiency can be improved reducing the dimension of the state vector/Liouvillian including only spin states that practically contribute to the spin dynamics. Alternatively, parallel computing architectures, such as graphics processing units (GPU) or computer clusters can be used.

| State space restriction
As compared with the Hilbert space, simulations in the Liouville space are more time-consuming since matrices have larger dimension. However, the dimension of the state vector can be reduced, using restricted state space and low correlation order techniques. 16,17,71,72 In Section 3, we saw that the spin space is given by the direct product of Pauli matrices, P 1 P 2 … P n , where P m is F I G U R E 8 Hartmann-Hahn condition profile as a function of 1 H radio frequency power. The matching condition occurs at 48 kHz under magic angle spinning the IST operator or a Pauli matrix of the spin m or the identity operator. In the case of 1/2-spins, the matrix dimension is equal to 2 N in the Hilbert space and equal to 4 N in the Liouville space, where N is the number of spins. The matrix dimension therefore grows exponentially with the number of spins in the system. However, the dimension of the space spanned by a realistic spin system trajectory is usually much smaller than the dimension of the full Liouville space. 16,17,71,73 The initial state vector is usually composed of the sum of single spin operators corresponding to theL Z orL AE states. For example, for a spin system composed of three spins, the state vector representing the longitudinal magnetization corresponds toρ =L After the action of pulses and evolutions, the single spin operators evolve into multispin operator (or high-spin order) corresponding to product states equal toL (three-spin order), and so on. In liquid state simulations, the presence of relaxation prevents high-order coherences and high levels of spin correlation from being populated to a significant extent. 73 Figure 9 shows an example of dynamics from one to seven-spin order during in a pulse-acquire experiment on a strychnine spin system. Figure 9 shows that spin order higher than four do not contribute to the dynamics and can be neglected. Simulation parameters and the code (explained step by step) used to obtain Figure 9 are given in Example 8 in Section B.8 of the supporting information. If we exclude correlations of more than k < n spins, the operators would get the form E P 1 … P 2 E … P k , corresponding to a number of IST operators or Pauli matrices scattered among the unit matrix operators, E. The dimension of the new state space is The cost of the system propagation in this new restricted space scales polynomially with respect to the total number of spins, and the simulation time is reduced. There are methods implemented in Spinach that can be applied a priori in order to decide which state in practice would contribute to the dynamics. This analysis is performed internally, these methods include:

| Interaction graph and connected subgraph analysis
The interaction graph is based on the local coupling density and on the fact that high-spin orders relax faster than low-spin orders and therefore do not accumulate in large quantity. 16 Even complicated magnetic resonance experiments do not often include spin order higher than four and correlations between remote spins so they can be dropped from the basis. For a given spin, the relevant states and coherences are those that involve its close neighbors on an interaction graph that can be expanded in a complete set of connected subgraphs (spin clusters) of a specified size composed by a F I G U R E 9 Evolution of spin correlations during a simulation trajectory spin and its nearest neighbors. The size of the interaction matrices is specified by the users in the spin system specification, these can include the product states between directly coupled spins or between spins closer than a given cutoff. For example, a protein system in liquid state has billions of J-couplings and dipolar interactions, but only couplings larger than 5 Hz and interactions within 5 Å between three or four closest spins are important. Hence, setting the interaction cutoff up to 5 Hz and the proximity cutoff up to 5 Å, the number of J-couplings and dipolar interactions is reduced to less than hundreds of thousands. In this way, the important states are included and the entanglements higher than the size of biggest cluster are excluded, and the simulation time gets drammatically reduced. In Spinach, the option to request complete basis set, which is only practical up to about 20 spins in Hilbert space and 10 spins in the Liouville space, can be specified as: bas.approximation='none'; The option to use incomplete basis sets can be specified as: bas.approximation='IK-0' or 'IK-1' or 'IK-2'; where the option 'IK − 0 0 includes all product statescorrelations between indirectly coupled spins up to a number of spins specified by a cutoff level:

bas.level
This basis set is used to drop correlations between very remote spins. The option 'IK − 1' includes all product states between directly coupled spins up to a number specified by the cutoff bas.level and correlation between indirectly coupled spins up to a number specified by a cutoff space level: bas.space_level that are closer together than a user-specified tolerance, representing the proximity cutoff radius in Angstroms:

sys.tols.prox_cutoff
The option 'IK − 2 0 includes, for every spin, all correlations with all directly coupled spins and correlations with up to a number of spins specified by a cutoff space level: bas.space_level that are closer together than the user-specified proximity cutoff. Fastest simulations use incomplete basis sets that are only available within the Liouville space.

| Accuracy of the state space restriction approximation
The applicability for the state space restriction approximation requires that the fraction of the magnetization leaking outside the restricted state space does not exceed the user-specified tolerances. 73 This can be fulfilled when where k is the cutoff level specified by the users-the maximum spin correlation order to be included in the basis, ξ is the user-specified tolerance, r is the average relaxation rate in the system, and h is the average spinspin coupling. The critical parameter is the h/r ratio of the largest spin-spin coupling to the smaller relaxation rate. Spin systems characterized by small h/r ratio, where the time scale of the relaxation process is smaller or comparable with the time scale of the dynamics, can be accurately described by low correlation order and require small basis sets (e.g., in liquid-state NMR systems). Spin systems characterized by large h/r ratio, where the time scale of the relaxation process is slower than the time scale of the dynamics, require larger basis sets (e.g., in solid-state NMR systems). Consequences of the state space restrictions on the simulation time were already shown in a previous work. 73 It is advisable to run trajan.m function (described in Section 3.6) to confirm that the basis set is sufficiently large, as shown in Figure 9.

| Trajectory level pruning and Krylov subspace analysis
After the state space restriction using the interaction graph analysis, if operators are still sparse, the state space can be further reduced by considering only the states that contribute to the system evolution under a pulse sequence. 17 The idea of the trajectory level pruning in restricted state spaces consists of finding the minimal basis sets that contribute to the simulation trajectory. All simulation trajectories have a finite duration and finite step length. The space spanned by trajectory vectors is called the Krylov subspace. Many spin states are not populated in a typical magnetic resonance simulation. Therefore, the dimension of the Krylov subspace is smaller than the dimension of the full state space, and the size of the simulation can be reduced by finding a way to map these subspaces. Once a trajectory under a given Liouvillian has been obtained, the minimal basis sets spanning this trajectory are obtained using orthogonalization techniques, for example, Gram-Schmidt, 74 single value decomposition (SVD), 75 and Lanczos. 46 In Spinach, the minimal basis set construction is considered complete when the next trajectory vector is entirely contained in the previous, or within a user-specified tolerance. This option is switched on by default. Systems involving fewer than five spins do not require sophisticated trajectory-level state space analysis, and the option should be disabled with the command: sys.disable={'trajlevel'}; An example in this setting in given in Example 7 in the supporting information.

| Zero track elimination
The problem described in Section 5.1.3 can be reformulated in a way that instead of considering states that appear in the Krylov subspace, and then projecting them into that subspace, we look at states that do not appear in the Krylov subspace and project them out. Liouvillians and the typical initial states (i.e.,Ŝ z andŜ AE ) are sparse matrices, in a typical pulse sequence, many elements of the state vector will remain zero, and dimensions corresponding to those states can be removed. The zero track elimination (ZTE) techniques 17 detect states that do not appear in the Krylov subspace and prunes them out. In Spinach, the ZTE is switched on by default and works by inspecting the first few steps in the system trajectory, end drops states that did not get populated. The ZTE is equivalent to the Krylov subspace analysis, but is more efficient. State space restrictions in addition to ZTE are useful in simulations of large systems of strongly coupled spins.

| Symmetry pruning
The symmetry pruning is based on symmetries and connectivities in the Liouville space and considers that only fully symmetric irreducible representations of a symmetry group contributing to the simulation. 76 Symmetric systems, like magnetically equivalent nuclei, are represented by symmetry groups that commute with the Hamiltonian. A simulation that starts in a state with a particular symmetry and evolves under a Hamiltonian that conserves that symmetry would never populate states of different symmetries. The Spinach kernels detect these states and generate a block-diagonal Liouvillian in symmetry-adapted basis sets. Each block can be simulated separately, leading to improvement in the simulation efficiency.

| Liouvillian path tracing
In simulations involving spins that do not interact with each other, the coherence order of each spin is conserved independently. The Liouvillian path tracing (PT) technique is used to identify disconnected subspaces in the state space, introduced by accidental symmetries and conservation laws, and propagates them individually. 13,76 This option is switched on by default in Spinach, and it can be disabled with following command:

sys.disable={'pt'};
It is used in Liouville space simulations of large densely coupled spin systems encountered in electron paramagnetic resonance spectroscopy and spin chemistry.

| Destination state screening
Further improvements can be made to the previous methods considering only spin system trajectories that evolve to or affect the detection state. The destination state screening 77 identifies system trajectories containing contribution from both the initial and the detection states and retain them. It can be set with the destination flag in the evolution.m function (described in Section 3.5.4).

| Parallelization in multicore systems
It is believed that parallel computing architectures, including GPUs, multicore central processing units (CPUs), and computer clusters, can improve the efficiency of simulations. This is not always true. The equation of the motion are integrated into software packages in a matrix-vector representation; parallel architectures can accelerate matrix-vector multiplications per second. A single NVidia Titan V card can perform 10 13 multiplications per second, and a modern CPU with 16 cores and 4-GHz clock frequency and 4 multiplications per cycle can perform 256 multiplication per second; these are indeed very good numbers. However, the key to the efficiency relies on latency and bandwidth. The former is the time taken for the data requested by the CPU to start arriving from the different memory banks, and the latter is the rate at which that data arrives. Every bite of memory is always 1 ns away from the CPU, and during that time, a CPU would wait rather than calculate. Parallel computing systems can improve the efficiency of calculations that involve linear memory access (e. g. fast Fourier transform), matrices larger than 200 × 200 (smaller matrices would spend more time in transit than multiplying), simulations involving billions of identical calculations. Few examples include: different orientations in powder average simulations 78 , different points of on a grid, 79 and different time slices in multidimensional NMR simulations. Several functions in Spinach that take advantage of parallel architectures are switched off by default. The GPU arithmetic on NVidia GPU cards can be switched on with the command: sys.enble={'gpu'}; Situations where the state spaces are dominated by a single large subgraph, for example, multidimensional liquid state NMR simulations, can be accelerated using the greedy parallelization that can be switched on with the command sys.enble={'greedy'}; Situations where there are repeated requests of expensive functions of the same argument, for example, matrix exponential in the propagator function, can be solved using hashing and caching algorithms. 78 These can be switched on with the command sys.enble={'chaching'}; In this procedure, a matrix is first hashed, meaning that a record is kept into the disk and then fetched when needed, rather than recomputed. This approach can be used when sufficient storage is available, in cases where the matrix exponentialization dominates the numerical cost of the simulation (e.g., optimal control, or in MAS simulations, where each rotor period contains identical Hamiltonian).

| CONCLUSION
In summary, it is possible to conclude that compared with the Bloch equations, the density operator and the Fokker-Planck formalisms represent a more general way to simulate the spin dynamics process like relaxation and kinetics. The Fokker-Planck theory allows simulations of complicated magnetic resonance experiments (SPEN-NMR, MRI, and MAS NMR,) where spin dynamics coexist with spatial dynamics. In the Liouville space, efficiency improvements can be obtained with state space restriction and low correlation orders methods, or through the utilization of parallel computer architectures. Current challenges in spin dynamics remain related to the matrix dimension, particularly in the solid state where the relaxation is slow, and in MRI involving metabolomics (composed by more than six spins) diffusing in a three-dimensional space. The problem in solid state may be solved considering that the powder averaging has the same effect as relaxation in liquid state. 71 In MRI simulations, there are some correlations between spin and spatial degrees of freedom that may be redundant in the presence of relaxation and diffusion, and may be removed from the Fokker-Planck space. This is a work in progress in our team. 14 Algorithms implemented in Spinach allow large-scale magnetic resonance simulations, including MRI, 14 powder average 67 and relaxation theories. 80 It also provides a software environment for pulse sequence development with the possibility of using the optimal control theory, 80,81 and supports parallel processing on GPU and computer clusters. Another advantage of Spinach is that it is an open source and allows direct access to an intuitive and well-documented code. Instructions on installation and practical information can be found in Kuprov 25 or on the website http://spindynamics.org/. Information on each function described in this tutorial can be found on the Spinach wiki page http:// spindynamics.org/wiki/index.php?title=Main_Page.