Linear-scaling ab initio simulations of spin diffusion in rotating solids

Spin dynamics simulations are used routinely in many areas of magnetic resonance, including nuclear magnetic resonance (NMR) and electron paramagnetic resonance (EPR), for the simulation of spectra, the prediction of the result of a particular pulse sequence, the design of new experiments, and the investigation of novel phenomena such as the mechanisms involved in dynamic nuclear polarization (DNP). Although many approximate spin dynamics methods have been developed, such as perturbation theory and average Hamiltonian theory,¹ in contrast to other areas of chemical physics, the exact solution to the spin-only Schrödinger equation is usually accessible. As a result of this, brute force exact calculations have become the norm.^2–11 For most applications, satisfactory results can be obtained from exact simulations using small spin systems (≤12 spin-1/2 nuclei; far less for quadrupoles), and these can be performed on most personal computers. Certain applications, however, require simulation of the concerted dynamics of far more spins. Examples include spin diffusion,^12,13 cross-relaxation,^14–16 DNP,^17–22 multiple-quantum spectroscopy,²³ and quantum computing,²⁴ for which hundreds, if not thousands, of spins are required to obtain physically meaningful results. Tackling such large spin systems using exact, or even perturbative, methods is problematic to say the least given that the dimensions of the density matrix scale as $∏spins(2I+1)$ in Hilbert space and $∏spins(2I+1)2$ in Liouville space. If we consider a spin system consisting of N spin-1/2 particles, this corresponds to matrix dimensions of 2^N and 4^N, respectively.

The exponential scaling of spin dynamics simulations drastically limits the size of the spin systems that can be studied; as a result, researchers have developed innovative approximate solutions to the many-body dynamics involved in large groups of spins. Notable examples include the concept of spin temperature,²⁵ the modeling of spin diffusion using classical diffusion equations^26–30 or as a multistep rate process,^31–34 and the simulation of DNP in large spin systems using a kinetic Monte Carlo algorithm.^35–37 

Recent efforts have been increasingly directed at reducing the size of the density matrix by removing unimportant, or unpopulated, states from the basis set such that fully ab initio spin dynamics simulations can be performed in spin systems containing hundreds of spin-1/2 nuclei.^38–45 These are known as restricted-state space methods and have been used to a great success for the simulation of liquid state NMR spectra of large spin systems, including whole proteins,¹⁶ and have also been applied in solids for the simulation of spin diffusion,^46–52 13C cross-polarization (CP) magic-angle-spinning (MAS) NMR spectra,⁵³ and ¹H MAS NMR spectra.⁵⁴ More recently, we have also applied one of these methods, in combination with a Monte Carlo optimization scheme, for the simulation of DNP in moderately sized spin systems.²² 

Of particular interest for this work is the specific application of restricted-state space methods that Dumez and co-workers used for the simulation of spin diffusion, under MAS conditions, in spin systems containing up to 144 ¹H spins.^47–52 Their method, coined LCL (low-order correlations in Liouville space), makes two key reductions to the Liouville space density matrix. First, since spin diffusion is a purely zero-quantum process, all >0-quantum coherences can be safely dropped from the basis set. Second, since the initial and final states of the density matrix are composed of single-spin operators, the numerous product operators with large spin orders are assumed to contribute less to the overall result and are thus dropped from the basis set. This same approximation is used in liquid-state spin dynamics simulations for which its success has been attributed to the slower buildup and faster decay of higher-spin order product operators.⁴² In the solid state, this is no longer the case with the high-spin order operators being very rapidly populated. When a powder average is performed under MAS conditions, however, the importance of these higher-spin order operators is drastically reduced. This was originally assigned to an unknown cancellation in the presence of MAS,⁵¹ but Edwards et al. were later able to show that this was in fact caused by the reduced propensity of the higher-spin order operators to being reconverted into observable single-spin terms.⁵³ 

In this work, we explore the use of further state space reduction for the accurate calculation of spin diffusion. Our approach is guided by the previously developed phenomenological spin diffusion models, which suggest that far smaller basis sets can be used for the ab initio simulation of spin diffusion. We then compare our results to those obtained from full LCL calculations. The approximations we propose are very similar to the restricted-state space methods of Kuprov, in particular, the IK-1 basis set,⁵⁵ and allow for several orders of magnitude faster simulations than LCL. Notably, these calculations scale linearly with the spin system size and should thus enable the ab initio simulation of spin diffusion in spin systems containing thousands of nuclei.

The density matrix is used to represent the instantaneous average of states in a given spin system, which evolves in time under the influence of the Hamiltonian. It can be expressed in two distinct formalisms, namely, Hilbert space and Liouville space, in which it is represented as a 2^N × 2^N matrix or a 4^N vector, respectively, when all N spins have a spin quantum number I = 1/2. Generally, Hilbert space calculations are considerably faster to compute and are thus applied nearly exclusively throughout the literature.² The notable exception to this rule concerns cases dealing with large, complex spin systems since Liouville space lends itself more readily to state space reductions and approximations.^56–58 Probably the most important example of this is the widely used product operator formalism.⁵⁹ In Liouville space, the basis vectors (⁠$B^r$⁠) that make up the density matrix are products of single spin operators,

where these operators (Î_i,r) can take the following definitions:

and q_r corresponds to the spin order of a given basis vector. The spin order of a basis vector is the number of spins that are not represented by identity (Ê_i) in the product operator. The density matrix is then the sum of these different basis vectors, and their respective time-dependent coefficients

The density matrix can be propagated using the Liouville-von-Neumann equation,

to calculate the density matrix at a later time point after the application of a particular Liouvillian (⁠$L^^$⁠),

In cases where the Liouvillian is time-dependent, such as MAS, the propagation can be performed using small time-independent increments,

where the time from 0 to t is separated into P time-independent segments.

This organizational scheme of the density matrix has two important features that can be exploited. First, since we generally only excite coherences with low orders and can only observe single-quantum, single-spin operators experimentally, Liouville space provides a means of removing the high-coherence-order operators that are expensive to calculate and contribute little to the outcome of the calculation. Second, the presence of an identity operator (Ê_i) in the product operator basis set means that, aside from coherence order, specific basis operators can also be expunged depending on their spin orders or the indices of the spins in the product.³⁸ 

Using these concepts Dumez and co-workers were able to drastically reduce the size of the density matrix and quantitatively model spin diffusion in rotating solids.^47–52 They defined their basis set as operators with a coherence order of 0 (which is exact in the absence of pulses) and low q_r values; the maximum q_r (q_max) being typically less than or equal to 5. We later used the same formalism to simulate MAS-DNP in moderate spin systems ab initio.²² In the case of solution-state NMR, Kuprov and co-workers also found success using far more drastic state space reductions.^38,39 For their basis, they selected operators having low q_r values consisting of spins that are close to one another in space. They also recently demonstrated that good agreement between experimental and theoretical ¹H MAS linewidths could be obtained using this approach.⁵⁴ 

In part due to our interest in simulating DNP processes in large spin systems, for which spin diffusion plays a crucial role, we were keen to determine to what extent the basis set could be further reduced without sacrificing computational accuracy with regard to simulating spin diffusion. To this aim, we will briefly discuss a simplistic kinetic 1D spin diffusion model. This model consists of a linear chain of n equally spaced spins. As the magnetization traverses through the chain from one end to the other, it can do so through a number of different pathways involving magnetization exchange between various numbers of different spins. In the two extreme cases, magnetization can be exchanged directly between the spins at each end of the chain, or it can hop sequentially from one neighboring spin to the next until it reaches the end. The right basis set for spin diffusion calculations should then be able to capture the most important magnetization transfer pathways while discarding those that are the slowest.

It is known from previous work that the classical exchange rate between two sites depends on the inverse sixth power of the distance between these spins.^31,32 To a reasonable approximation, the magnetization is then transferred between two spins as follows:

where A is a phenomenological constant that encapsulates the MAS dependence and other effects, t is time, and r is the internuclear distance. In the case of a transfer occurring in n equal steps r, for instance, in a linear chain of spins, the overall magnetization transferred obeys the following series:⁶⁰ 

This expression can be rewritten using incomplete gamma functions as follows:

In this model, the fastest transfer pathway will always be the one that does the shortest jumps (n = N − 1), whereas the slowest will be the direct pathway (n = 1) due to the aggressive n⁶ dependence. Determining at which point the longer, more direct, transfers become insignificant should yield the most efficient basis set for the model. For this comparison, we use Eq. (9) to calculate at which point in time the transfer through the fastest pathway is 90% complete. We can then determine how much magnetization could be transferred during the same time through any of the competing pathways. For instance, during the time it takes for 2 successive jumps to transfer 90% of the magnetization, only 4.5% can be transferred through the direct pathway and 0.5% when there are 2 intermediary spins (n = 3), assuming that each pathway consists of equal steps. This simplistic model thus suggests that the dipolar coupling between distant spins can be safely ignored when simulating spin diffusion since the fastest pathway will correspond to the one following a large number of short-range magnetization transfers. Note that this model disregards backwards transfers and treats each pathway as though there were no competing pathways present. This is necessary in order to distinguish and compare the polarization transferred through each pathway.

The above analysis demonstrates that Kuprov’s approach of restricting the state space to neighbors within a given distance from each other^38,39 should work well for modeling spin diffusion. An important caveat to be made here is that the relationships presented above for the spin diffusion are independent of the ¹H density. Thus, while a cutoff based on distance works well for the calculation of the NMR spectrum,^53,54 for spin diffusion, the nearest neighbors rather than distance should be used as the selection criterion for the basis set in order to properly describe sparsely populated or inhomogeneous ¹H networks. This proposed basis set can be considered an IK-1(1, q_max) basis set employing a distinct maximum neighbors’ cutoff (N_max) as opposed to a cutoff radius.⁵⁵ 

To simplify the organization of a very sporadically sparse density matrix, we opted for the creation of a neighbors’ list wherein the nearest N_max spins to each spin are listed. A basis operator is then kept within the basis set if (1) it has a q_r value smaller than a given cutoff, (2) it has a total coherence order of zero, and, additionally, (3) all the nonunitary operators in the basis operator belong to spins that are on each other’s respective neighbor list. An example is given in Fig. 1.

Although we,²² and others,^41,51 have used complex indexing schemes when performing LCL calculations in order to minimize the memory requirements, this is difficult to do in a general way for the locally restricted model described here. As such in our software, we have ordered the operators according to (1) their spin with the smallest index and (2) their operator type. The coefficients (b_r) are then stored within these groups using a hash table, which is a fast sparse storage format. The explicit manner in which this is achieved is given in the supplementary material.

The propagation of the density matrix can then be conveniently performed using the Suzuki-Trotter algorithm, also used by Dumez et al.,⁴⁹ by simply cycling through the nonredundant spin pairs in the neighbor list,

where N is the number of spins in the simulation and $L^i,j$ is the dipolar coupling Liouvillian describing the coupling between spins i and j. If necessary, a second (or higher) order Suzuki-Trotter scheme can also be applied⁴⁹ although this doubles the computational time and, in many cases, yields identical results,

Each of the $exp(L^^i,jt)$ operations in Eqs. (10) and (11) affects a large number of product operators. They are then separated into a series of 4-dimensional subspaces that are evaluated sequentially by multiplying the subspace with the following rotation matrices:

for subspaces of the form [Î_iz$B^r$⁠, Î_jz$B^r$⁠, I_i+Î_j−$B^r$⁠, Î_i−Î_j+$B^r$] and

for subspaces of the form [Î_i±$B^r$⁠, Î_j±$B^r$⁠, Î_i±Î_jz$B^r$⁠, Î_j±Î_iz$B^r$] as described previously, where $B^r$ is an operator excluding spins i and j.⁴⁹ The dipolar frequency, ω_D,i,j, is calculated using Wigner rotation matrices,

where R_DD,i,j corresponds to the dipolar coupling constant between spins i and j,

and three sets of Euler angles are used to relate the orientation of the crystal within the rotor frame $(0,θ,ϕ)$⁠, the dipolar coupling tensors to the crystal frame (α, β, γ), and the rotor frame to the lab frame (ω_rt, θ_r, 0). In Eq. (15), γ_i corresponds to the magnetogyric ratio of spin i, and r_i,j corresponds to the internuclear distance between the spins i and j. Due to its close relationship with the LCL method, we have termed this approach locally restricted LCL or LR-LCL.

Chemical shifts may also be included in the simulation with relative ease by cycling through the basis operators, as in Eq. (10), and applying the following operation to the coefficients of each pair of raising and lowering operators Î_i+ Î_j−:

In Eq. (16), ω_i corresponds to the resonance frequency of spin i in the rotating frame.

To validate the LR-LCL method, we chose to simulate the ¹H–¹H spin diffusion along a linear alkane, which is reminiscent of the polymer that was simulated by Kuprov et al.³⁸ Such a model system provides a simple means of following the spin diffusion across long distances while also minimizing the number of spins in the system. This is important given that the equivalent LCL calculations needed for benchmarking LR-LCL can be quite computationally demanding for larger spin systems. This is particularly the case when following spin diffusion over long distances since the evolution times can be quite long, and the propagators cannot be stored and reused due to their size. Note that the LCL method has been extensively benchmarked and shown to provide results that approach those of exact calculations.^47,48 As such we will not compare the results obtained from LR-LCL calculations to those of exact calculations and simply rely on the results from full LCL calculations.

We focus, in particular, on octadecane, which contains 38 ¹H spins, since this is the longest alkane for which spin diffusion could be modeled, within reasonable time, using LCL. The number of basis vectors for q_max = 4 and 5 LR-LCL basis sets, as well as the ratio between the two, is plotted in Fig. 2 as a function of N_max. Note that for octadecane, the full LCL basis set corresponds to N_max = 37. Given that the number of matrix-vector operations depends linearly on the number of basis vectors and that all matrix-vector operations have the same dimensionality when using Suzuki-Trotter propagation [Eqs. (11)–(13)], the computational time also depends linearly on the size of the basis set. As a result, several orders of magnitude in computational time savings are obtained when performing LR-LCL simulation. Furthermore, when we locally restrict the basis sets, the cost of increasing the accuracy of the simulation by increasing q_max is also dramatically reduced [see Fig. 2(b)] by nearly one order of magnitude.

The calculated spin diffusion curves for octadecane, when starting with the magnetization at one end of the chain and detecting the magnetization at the opposite end, are shown in Fig. 3 as a function of N_max. For the simulations, the spinning frequency was set to 10 kHz and time steps of 1 µs were applied up to a total simulation time of 20 ms. Powder averaging was performed using 66 orientations, selected using the REPULSION scheme,⁶¹ and 3 γ angles. As can be seen, the red curves, which correspond to the LR-LCL results, quickly converge to the LCL result as N_max is increased with little improvement obtained when using an N_max value greater than 15. Note that this is quite significant given that the N_max = 15 LR-LCL calculation took only 66 min to complete on an Intel E5-2620 12-core processor whereas the LCL calculation required 3.2 days of computing time.

The convergence can also be followed by plotting the average root mean squared deviation (RMSD) between the LR-LCL and LCL polarizations of all spins at all time points [Fig. 4(a)],

Interestingly, the convergence observed in Fig. 4(a) is not continuous but rather occurs in a series of discrete jumps every time N_max is increased by 4. This corresponds to the times when the basis set is expanded to include the nearest proton of a new carbon neighbor. The inclusion of the interactions between the farther protons on adjacent carbons appears to be less significant, given that thepolarization between the two ¹H spins on a CH₂ equilibrates very rapidly.

We have additionally repeated the calculations with all the atomic coordinates scaled by a factor of 2 [Fig. 4(b)] to demonstrate that indeed it is the number of neighbors which is important to consider when selecting the basis set, rather than the distances between the spins. As can be seen, the final jump in accuracy again occurs when N_max reaches a value of 15 and the spins are connected to ¹H spins situated within 4 C–C bonds away.

In order to rationalize the optimal basis set for this structure, we have applied the rate-based unidirectional spin diffusion model described in Eqs. (7)–(9). We treated both protons on a given carbon as a single spin, due to their rapid equilibrium, and set the spin diffusion rates according to the expected r⁻⁶ relationship. We then calculated the time it would take for the polarization transfer to reach 90% through the n-step process [using Eq. (9)], from a ¹H situated n C–C bonds away, and then calculated how much polarization could be transferred through the equivalent single-step transfer in that same amount of time. The results are plotted in Fig. 4(c). Interestingly, for n = 2, direct polarization transfer is faster than the relayed process. This is the case since the ¹H spins on adjacent carbons were 2.512 Å apart, while the nearest ¹H on the subsequent carbon is only 2.531 Å away due to the zig-zag structure of the molecule. There is then a sharp drop for n = 3 due to the significantly increased distance for the direct transfer and a lesser drop for n = 4. All direct pathways with n > 4 are insignificant. Note that this is in good agreement with the LR-LCL results since, for this model, n = 4 would correspond to an N_max value of 15–17, which we identified as the most time efficient, accurate, basis set. Note that the linear alkane studied here is a particular case and that N_max would need to be increased to 25–35 in order to get the same coverage in a 3D distribution of spins. In those cases, a similar analysis as was used here could be applied to aid in the selection of a basis set. Ideally, however, calculations should be performed with different N_max values to ensure that convergence is reached.

Finally, we performed a series of spin diffusion calculations on a variety of linear alkanes having between 10 and 50 carbons, the longest alkane having 102 ¹H spins, in order to test the scalability of the LR-LCL method. In Fig. 5(a), the average time required to perform a single time increment, on a single crystallite orientation, is plotted as a function of the number of ¹H spins. Clearly, the LR-LCL method scales linearly with the spin system size, which should enable spin diffusion calculations on far larger spin systems than was previously possible. Simulations performed on a 3D system with 100–700 spins are also given in the supplementary material showing that the linear scaling is indeed general. In Fig. 5, we also compare the size of the LR-LCL basis set, with N_max = 15 and q_max = 4, to the exact [Fig. 5(b)] and LCL ones [Fig. 5(c)]. Between 10 and 60 orders of magnitude in faster performance are obtained, with respect to exact calculations, for spin systems of this size and between 1 and 4 orders of magnitude when compared to LCL, with further improvements obtained in the case of larger spin systems. To put this into perspective, while the calculation of a single time step for a particular orientation took on average only 3 ms for the C50 molecule when using LR-LCL, the same calculation would require 13 s with LCL and 3 × 10⁴⁰ years if done using the full basis set.

We have tested the utility of locally restricted basis sets for low-order correlations in Liouville space (LR-LCL) calculations of spin diffusion. These basis sets are similar to those that have been proven effective for solution state NMR studies, in particular, the IK-1 basis set, albeit using a neighbors’ cutoff instead of a cutoff radius. The simulated spin diffusion curves were found to quickly converge to those obtained using LCL as the basis set size was increased to include interactions and product operators between a greater number of neighboring spins. Unlike exact and LCL calculations, which scale exponentially and polynomially, respectively, as the spin system size is increased, the LR-LCL method scales linearly and enables the ab initio calculation of spin diffusion in solids undergoing MAS with near-exact precision for spin systems containing upwards of 1000 spins.

All LR-LCL calculations were performed using an in-house code written in C/C++ and run on a local cluster composed of Intel E5-2620 12-core processors. Each calculation used a single processor, with the powder averaging being run in parallel on the different cores using OpenMP. The molecular polyethylene models were optimized using Avogadro.⁶² The MAS spinning frequency was set to 10 kHz, and the rotor period was separated into 100 time increments of 1 µs. A total of 200 rotor periods were simulated, equaling a spin diffusion time of 20 ms, with the values of N_max and q_max set in accordance to those described in the text.

Details pertaining to the storage format of the density matrix, as well as additional timing calculations on a three-dimensional model, are provided in the supplementary material.

Dr. Lin-Lin Wang is thanked for many insightful discussions. This research was supported by the U.S. Department of Energy (DOE), Office of Science, Basic Energy Sciences, Materials Science and Engineering Division. The Ames Laboratory is operated for the DOE by Iowa State University under Contract No. DE-AC02-07CH11358.