- Email: [email protected]
Dynamics of large nuclear-spin systems from low-order correlations in Liouville space
Chemical Physics Letters 477 (2009) 377–381
Contents lists available at ScienceDirect
Chemical Physics Letters journal homepage: www.elsevier.com/locate/cplett
Dynamics of large nuclear-spin systems from low-order correlations in Liouville space Mark C. Butler, Jean-Nicolas Dumez, Lyndon Emsley * Université de Lyon (CNRS/ENS Lyon/UCB Lyon1), Centre de RMN à très hauts champs, 5 rue de la Doua, Villeurbanne, France
a r t i c l e
i n f o
Article history: Received 22 May 2009 In final form 3 July 2009 Available online 8 July 2009
a b s t r a c t Simulation of the coherent dynamics of a large lattice of strongly coupled nuclear spins is demonstrated. The evolution of N spins is simulated in a reduced Liouville space X which excludes coherences involving more than k spins, where k N is chosen to minimize the dimension of X without introducing unacceptable error. For polarization transfer within a spinning lattice of protons, comparison with exact simulations shows that k = 4 allows accurate simulation of the experimentally relevant powder-averaged curves. A simulation of polarization transfer among 144 protons is presented. Ó 2009 Elsevier B.V. All rights reserved.
1. Introduction It is astounding to consider that while sophisticated multidimensional nuclear magnetic resonance (NMR) methods are today routinely applied to solids, it is currently impossible to simulate exactly the dynamics of the coupled many-spin systems present in solids that generate the observables in these experiments. NMR methods play a central role in the characterization of solids, with applications ranging from studies of the atomic level structure and dynamics of biological macromolecules [1] to quantum computing [2], and in most of these areas, the inability to simulate the dynamics of a large lattice of coupled spins is a principal barrier to further progress [3,4]. Several approaches can be used to model such dynamics. Perturbation theory can be used to characterize the average effect of the network of couplings on experimental observables, yielding models in which the number of relevant variables is significantly reduced [5,6]. While such models allow simple analysis of some experiments, they require the use of empirical parameters and cannot account for effects that arise from the coherent dynamics of the lattice. Another approach to modelling involves numerical propagation of the full density matrix, using software that can simulate a wide range of NMR phenomena and methods [7–10]. Highly optimised algorithms [10,11], in combination with methods exploiting translational symmetry [4] in the solid state, have been used to perform exact simulations including up to 15 spins. Simulating systems of this size does yield insight into various fundamental phenomena [12–15] but is insufficient to capture large-scale effects, for which empirical models have to be added [5,12]. Simulations involving the full density matrix are limited by the fact that the dimension of Liouville space scales exponentially with * Corresponding author. Fax: +33 478896761. E-mail address: [email protected] (L. Emsley). 0009-2614/$ - see front matter Ó 2009 Elsevier B.V. All rights reserved. doi:10.1016/j.cplett.2009.07.017
sample size. For a system of N spins I = 1/2, for instance, Liouville space has dimension 4N. One method of circumventing this obstacle is to use a reduced Liouville space to simulate the coherent dynamics of large spin systems. An example of this approach was first introduced by Bruschweiler and Ernst [16], who used a ‘cogwheel’ model of spin motion in a solid to simulate the short-term dynamics of a chain of 10 spins. More recently, a scheme for reducing Liouville space for liquid-state NMR simulations has been applied to large spin systems [17,18]. In this scheme, small clusters of nearest-neighbor spins are defined, and spin i interacts directly with spin j in the simulation only if the two spins belong to the same cluster. However, in the case where each spin is strongly coupled to many other spins, e.g., a strongly coupled lattice of protons in a solid, the applicability of this approach has not been established. In this Letter we introduce a method of simulating large, densely coupled nuclear-spin systems in a reduced Liouville space. Using this approach we are able to simulate the medium-term coherent dynamics of a powder sample under magic-angle spinning for a crystalline lattice of 144 protons. This demonstration, which constitutes an order-of-magnitude increase in the number of spins for which such dynamics have been simulated, suggests that the structure and dynamics of solids may be probed by directly studying the dependence of experimental observables on coherent interactions involving many spins.
2. Simulation scheme 2.1. Motivation Exponential scaling of Liouville space with sample size can be understood as resulting from the vast number of ways in which correlations can develop within subcollections of particles. If the dimension of a 1-particle Liouville space is m, then a space of
378
M.C. Butler et al. / Chemical Physics Letters 477 (2009) 377–381
dimension mN is sufficient to characterize a system of N particles that are constrained to remain statistically independent, while the full Liouville space that includes all correlations between the particles has dimension mN. It is the requirement that all possible correlations between the particles be included in the computation that renders simulations of the full Liouville space of large systems prohibitively demanding. In general, however, only a small fraction of Liouville space is observed experimentally (either directly or indirectly), and the possibility exists that accurate simulations of large quantum systems could be performed in reduced Liouville spaces which include both the degrees of freedom monitored experimentally and some ‘surrounding space’ containing the degrees of freedom most strongly excited during the experiment. As an example, consider the problem of simulating spin diffusion (i.e., transfer of spin polarization) within a dipole-coupled lattice of protons. This problem is today central to the determination of structure in drug molecules [19–21] and solid proteins [22,23] by NMR crystallography. Since each spin i interacts strongly with a relatively large number of spins j, many distinct pathways for polarization transfer exist, and accurate simulation of spin diffusion requires a model in which the network of couplings is sufficiently dense to include the important transfer pathways. Note, however, that simultaneous transfer between these polarizations corresponds to motion along an even larger number of trajectories in Liouville space; it is conceivable that interference between the different trajectories in Liouville space could limit the net effect of high-order coherences, thereby enabling accurate simulation in a reduced space. An analogous effect was observed in multiple-quantum spectra [24] when phase interference among high-order coherences hindered the observation of their excitation [25]. Note as well that most solid-state NMR experiments involve magic-angle spinning (MAS) and are carried out on powder samples, and that these experimental conditions might shrink the region of Liouville space which effectively contributes to detected observables. During MAS, the modes (eigenoperators) of Liouville space are continually changing, and if the time scale of the sample rotation is sufficiently close to that of the spin dynamics, MAS could distribute the excitation which spreads through Liouville space among many modes whose contribution to the experimental signal could interfere destructively. Consistent with this hypothesis is the finding of Spiess and co-workers that sample rotation at frequencies of the order of the highest couplings often allows analysis of sideband intensities using only two-spin coherences [26]. Similarly, distinct crystallites in general have distinct sets of modes, and the evolution of an experimental observable for a given static crystallite can be written as a sum of oscillating terms associated with the modes of that crystallite. By summing contributions from distinct sets of modes, powder averaging could increase the efficiency of destructive interference among the modes excited under MAS.
the density operator while respecting the constraint that it remain within X (i.e., elements of the full Liouvillian that would couple operators within X to those lying outside of X could be considered formally set to zero). A first-order scheme was used to compute the evolution of the density matrix r under the action of LX. Given a basis fri g of the reduced space, the density operator can be expressed as a vector of coefficients si , where
r¼
X
s i ri :
The first-order change in r that occurs during a time step Dt is
Dr ¼ L X r Dt ¼
X ðsi DtÞLX ri :
ð1Þ
For simulations of polarization transfer within an N-proton system, each basis element ri was chosen as a product of N single-spin operators, with at most k < N of the single-spin operators different from the identity. Explicit formulas for the operators LX ri appearing in Eq. (1) were derived, and first-order numerical propagation of r was implemented by using these formulas to compute r þ Dr. Use of the derived formulas for LX ri eliminated the need to store LX as a matrix, and so the factor limiting the size of the system that could be simulated was the memory required to store r as a vector of coefficients, rather than the memory needed for an evolution matrix acting on that vector. Note as well that for fixed k, the memory required for this simulation scheme scales polynomially with N rather than exponentially. This improved scaling is a key feature for the simulation of large systems. For a solid sample containing protons in a field of a few Tesla or more, the dominant spin–spin interactions are governed by the secular homonuclear dipolar Hamiltonian:
HD ¼
X i
l0 c2 h ð1 3 cos2 hij Þð3Iiz Ijz Ii Ij Þ: 8pr3ij
Here Ii and Ij are nuclear-spin operators, rij is the internuclear distance, hij is the angle between the internuclear vector and the static magnetic field, and c is the gyromagnetic ratio. For simplicity, additional terms in the Hamiltonian were neglected in these simulations, although it would be straightforward to include them. Since solid-state NMR experiments are most often performed on powder samples using MAS, sample rotation and powder averaging were included in some simulations. Rotation of the sample was simulated by treating the Hamiltonian as time-dependent and updating it repeatedly during the numerical propagation of the density matrix. Powder averaging [27] was performed by summing observables over a set of simulations carried out with selected crystallites, using the ZCW scheme. Geometric structures were obtained from published organic crystal structures. For all simulations, the initial density matrix was set to r(0) = Iiz for a single proton i, and the evolution of hIjz(t)i was recorded for all protons j in the system.
2.2. Implementation
3. Results and discussion
We have tested these hypotheses by comparing exact simulations of polarization transfer within a lattice of N protons to simulations performed in a reduced Liouville space X that included only low-order correlations. All products of single-spin operators involving more than k correlated spins were excluded from X, where k < N is a small integer. Given an initial state defined by polarizations on certain spins, exact simulations of the evolution under a secular spin Hamiltonian can be performed within the zero-quantum subspace Z of the full Liouville space [16], and we defined reduced Liouville spaces by excluding from Z all coherences involving more than k correlated spins. The restriction LX of the full Liouvillian to a reduced space X was used to propagate
Exact simulations carried out with the SPINEVOLUTION software [10] on molecules containing up to 12 protons were compared with simulations performed in reduced Liouville spaces. Fig. 1 illustrates typical results obtained for a single orientation of a static molecule containing 12 protons, b-L-aspartyl-L-alanine. The simulations using only low-order correlations in Liouville space, which we dub LCL, reproduce the exact evolution for a very short initial time period, but after 100 ls, large errors develop in the simulations having k = 5 (dashed curves), i.e., in the simulations for which product operators involving more than 5 correlated spins were excluded. The simulations performed in a further reduced space, k = 4 (dotted curves), show pathological departure
379
M.C. Butler et al. / Chemical Physics Letters 477 (2009) 377–381
1 H H H H H H N O
12 spins, static
0.8
O
< Iiz > / || Iiz ||
O
H
H H
N H H O
O
H
0.6 0.4 0.2 0
-0.2 0
0.2
0.4 0.6 time (ms)
0.8
1
Fig. 1. Simulations of polarization transfer for a single molecule of b-L-aspartyl-Lalanine (CSD entry: fumtem). Exact simulations (solid) of this system of 12 protons are compared to simulations performed in reduced Liouville spaces that exclude coherences involving more than k = 4 spins (dotted) or k = 5 spins (dashed). The initial density matrix is r(0) = ICH2,z, and the observables hICH2,zi (red) and hICHasp,zi (blue) are plotted. The time step in the simulation is 1 ls. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
from the exact simulations (solid curves), and from the curves corresponding to k = 5. The relatively simple oscillations in the k = 4 curves are consistent with the hypothesis that the initial state of the spin system is dominated by the contribution from a few spurious modes in the corresponding reduced space. Similar short-term accuracy was reported using a cogwheel model to propagate the density matrix in a reduced Liouville space that excluded product operators involving more than 3 correlated spins [16]. The error present in these LCL simulations shows that the excluded regions of Liouville space are relevant for simulations of a single static crystal. The accuracy of the LCL simulations is significantly improved by adding spinning to the simulation, as shown in Fig. 2. While in the static case, agreement was obtained for about 50 ls, the agreement with exact simulations (solid curves) is pre-
served under 10 kHz MAS for longer time periods (around 250 ls), and the decay of the initial polarization, although not reproduced exactly, is qualitatively correct for both k = 4 (dotted curves) and k = 5 (dashed curves). The k = 4 curves appear to include contributions from many interfering modes, rather than showing the simple oscillations associated with the evolution of a small number of modes. The accuracy of the LCL simulations is further improved by the combination of spinning and powder averaging, as shown in Fig. 3, in which the exact curves are now quantitatively reproduced using the LCL scheme for times up to 1000 ls, corresponding to the whole of the experimentally relevant timescale. The Liouville space that includes 4-spin correlations (k = 4) is now sufficient to reproduce detailed features of the spin diffusion curves. LCL simulations were compared with exact simulations performed for a collection of 25 10-proton organic molecules and 15 12-proton organic molecules randomly selected from the Cambridge Structural Database (CSD) [28]. For a system containing N protons, an exact simulation was performed for each of the N distinct initial states r(0) = Iiz, and each of these simulations yielded N curves hIjz ðtÞi=kIjz k; i.e., 4660 exact curves were compared against LCL curves. The error in the LCL simulations was characterized by calculating a root-mean square (rms) error for each curve, and the distribution of rms errors for the set of 10-proton curves and the set of 12-proton curves were found to be to be similar. Among the k = 4 and k = 5 curves, respectively, 78% and 95% had an rms error of 0.01 or less, as shown in Fig. 4. (The Supplementary material includes all curves used for this test of the LCL simulation scheme.) These results support the hypothesis that for spinning organic powders, the dynamics of spin diffusion among protons can be accurately modelled using a small fraction of the full Liouville space. In the case of a 12-proton system, for instance, Table 1 shows that the dimension of the k = 4 reduced Liouville space is less than 0.1% as large as that of the full Liouville space. An additional finding is that reduced Liouville spaces X adequate for simulating these dynamics can be obtained simply by excluding all products of single-spin operators involving more than k = 4 or k = 5 correlated spins. Note that for simulations performed in such a space, excitation that develops in k-spin coherences can only be
1 1 12 spins, powder, MAS
0.8
12 spins, MAS < Iiz > / || Iiz ||
< Iiz > / || Iiz ||
0.8 0.6 0.4
0.4 0.2
0.2 0
0.6
0 0
0.2
0.4 0.6 time (ms)
0.8
1
Fig. 2. Simulations of polarization transfer for a single molecule of b-L-aspartyl-Lalanine under 10 kHz MAS. Exact simulations (solid) of the 12-proton system are compared to simulations performed in reduced Liouville spaces that exclude coherences involving more than k = 4 spins (dotted) or k = 5 spins (dashed). The initial state, observables, and time step are the same as in Fig. 1.
0
0.2
0.4 0.6 time (ms)
0.8
1
Fig. 3. Powder-averaged simulations of polarization transfer for a single molecule of b-L-aspartyl-L-alanine under 10 kHz MAS. Exact simulations (solid) of the 12proton system are compared to simulations performed in reduced Liouville spaces that exclude coherences involving more than k = 4 spins (dotted) or k = 5 spins (dashed). The initial state, observables, and time step are the same as in Fig. 1, and a ZCW50 set of 50 orientations is used for the powder average.
380
M.C. Butler et al. / Chemical Physics Letters 477 (2009) 377–381
4000
number of curves
k=4
k=5
2000
0
0.02
0.04
0.06
0.08 0 rms error
0.02
0.04
0.06
0.08
Fig. 4. The rms error for 4660 polarization curves obtained from 25 10-proton systems and 15 12-proton systems. For each curve, exact simulations were compared with simulations performed in reduced Liouville spaces that excluded coherences involving more than k = 4 or k = 5 spins.
Table 1 The dimensions of Liouville spaces that can be used for simulation of polarization transfer within 10-proton and 12-proton systems. The full Liouville space includes all possible spin correlations among the N protons. Given an initial state defined by polarizations on certain protons, exact simulations of the evolution under a secular spin Hamiltonian can be performed within the zero-quantum subspace. The reduced Liouville spaces are subspaces of the zero-quantum subspace, formed by excluding product operators involving more than k correlated spins. Full Liouville space
Zero-quantum subspace
Reduced spaces k=4
10 protons
1:0 10
12 protons
1:7 107
5
6
1:8 10
5:0 10
1:8 104
2:7 106
1:1 104
5:2 104
1 144 spins, powder, MAS
< Iiz > / || Iiz ||
0.8 0.6 0.4 0.2 0 0
0.2
0.4 0.6 time (ms)
0.8
k=5 3
1
Fig. 5. Powder-averaged simulations of polarization transfer in a system of 144 protons (solid) and a system of 12 protons (dotted) under 10 kHz MAS. The simulations are performed in a reduced Liouville space that excludes coherences involving more than k = 4 spins. The 12-proton system is a single molecule of b-Laspartyl-L-alanine (CSD entry: fumtem), while the 144-proton system is a rectangular block of three primitive unit cells of the same molecule. For the system of 144 spins, the minimum image convention [29] is used to define periodic boundary conditions, and the observables are summed over crystallographically equivalent spins. The initial density matrix is r(0) = ICH2,z, and the observables hICH2,zi (red) and hICHasp,zi (blue) are plotted. The time step in the simulation is 1 ls, and a ZCW50 set of 50 orientations is used for the powder average. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
reflected back to lower-order coherences, rather than propagated outward to coherences involving additional spins. We are investigating the way in which such errors introduced into the dynamics by the reduction scheme do not contribute significantly to the spin diffusion curves. A key feature of the LCL simulation scheme is that its favourable scaling properties allow simulation of large systems. As a demonstration, we have simulated powder-averaged spin diffusion curves for a spinning lattice of 144 protons, i.e., the 12 molecules contained in three primitive unit cells of b-L-aspartyl-L-alanine. Although a detailed comparison of the 12-proton simulation and the 144-proton simulation is outside the scope of this Letter, we note that Fig. 5 shows the expected qualitative similarity between the curves for the two systems, although quantitative differences are present. The full Liouville space for the 144-proton system has dimension 5 1086, which makes the full density matrix impossible to simulate. The reduced Liouville space has dimension 3 108, and the LCL simulation for a single orientation used 1 GB of random-access memory and required only 8 h of computing time on a single 2.8 GHz processor. 4. Conclusion We have found that exact powder-averaged simulations of polarization transfer within spinning organic crystals can be quantitatively reproduced over the full experimentally relevant timescale using reduced Liouville spaces that include only low-order correlations. For 10-proton and 12-proton organic structures, exact simulations of polarization transfer under the secular homonuclear dipolar Hamiltonian were quantitatively reproduced by simula-
M.C. Butler et al. / Chemical Physics Letters 477 (2009) 377–381
tions performed using reduced basis sets that excluded product operators involving more than 4 correlated spins. For the organic structures tested, the dimensions of the reduced Liouville spaces are orders of magnitude smaller than the dimension of the full Liouville spaces; however, the combination of magic-angle spinning and powder averaging enables accurate simulation of proton spin diffusion within the reduced spaces. Since the dimensions of the reduced space scale polynomially rather than exponentially with the number of spins, the simulation of large crystalline lattices is feasible, as demonstrated by a simulation of polarization transfer among 144 protons. This demonstration suggests that simulations including only low-order correlations in reduced Liouville spaces may enable new studies of the structure and dynamics of solids by revealing the dependence of experimental observables on coherent interactions involving many spins. Appendix A. Supplementary material Supplementary data associated with this article can be found, in the online version, at doi:10.1016/j.cplett.2009.07.017. References [1] A.E. McDermott, Curr. Opin. Struc. Biol. 14 (2004) 554. [2] L.M.K. Vandersypen, I.L. Chuang, Rev. Mod. Phys. 76 (2004) 1037. [3] W.X. Zhang, N. Konstantinidis, K.A. Al-Hassanieh, V.V. Dobrovitski, J. Phys.: Condens. Matter 19 (2007).
381
[4] P. Hodgkinson, D. Sakellariou, L. Emsley, Chem. Phys. Lett. 326 (2000) 515. [5] B.H. Meier, Adv. Opt. Magn. Reson. 18 (1994). [6] M. Ernst, B.H. Meier, in: I. Ando, T. Asakura (Eds.), Solid State NMR of Polymers, Elsevier, Amsterdam, 1998, p. 83. [7] M. Bak, J.T. Rasmussen, N.C. Nielsen, J. Magn. Reson. 147 (2000) 296. [8] F.S. Debouregas, J.S. Waugh, J. Magn. Reson. 96 (1992) 280. [9] P. Hodgkinson, L. Emsley, Prog. Nucl. Magn. Reson. Spectrosc. 36 (2000) 201. [10] M. Veshtort, R.G. Griffin, J. Magn. Reson. 178 (2006) 248. [11] V.E. Zorin, M. Ernst, S.P. Brown, P. Hodgkinson, J. Magn. Reson. 192 (2008) 183. [12] D. Sakellariou, P. Hodgkinson, L. Emsley, Chem. Phys. Lett. 293 (1998) 110. [13] V.E. Zorin, S.P. Brown, P. Hodgkinson, J. Chem. Phys. 125 (2006) 144508. [14] V.E. Zorin, B. Elena, A. Lesage, L. Emsley, P. Hodgkinson, Magn. Reson. Chem. 45 (2007) S93. [15] J.S. Waugh, Mol. Phys. 95 (1998) 731. [16] R. Bruschweiler, R.R. Ernst, J. Magn. Reson. 124 (1997) 122. [17] I. Kuprov, N. Wagner-Rundell, P.J. Hore, J. Magn. Reson. 189 (2007) 241. [18] I. Kuprov, J. Magn. Reson. 195 (2008) 45. [19] B. Elena, L. Emsley, J. Am. Chem. Soc. 127 (2005) 9140. [20] E. Salager, R.S. Stein, C.J. Pickard, B. Elena, L. Emsley, Phys. Chem. Chem. Phys. 11 (2009) 2610. [21] B. Elena, G. Pintacuda, N. Mifsud, L. Emsley, J. Am. Chem. Soc. 128 (2006) 9555. [22] F. Castellani, B. van Rossum, A. Diehl, M. Schubert, K. Rehbein, H. Oschkinat, Nature 420 (2002) 98. [23] C. Wasmer, A. Lange, H. Van Melckebeke, A.B. Siemer, R. Riek, B.H. Meier, Science 319 (2008) 1523. [24] M. Munowitz, A. Pines, Adv. Chem. Phys. 66 (1987) 1. [25] Y.S. Yen, A. Pines, J. Chem. Phys. 78 (1983) 3579. [26] C. Filip, S. Hafner, I. Schnell, D.E. Demco, H.W. Spiess, J. Chem. Phys. 110 (1999) 423. [27] M. Eden, Concepts Mag. Reson. A 18A (2003) 24. [28] F.H. Allen, Acta Crystallogr. Sect. B-Struct. Sci. 58 (2002) 380. [29] M.P. Allen, D.J. Tildesley, Computer Simulation of Liquids, Oxford Science, Oxford, 1987.
Contents lists available at ScienceDirect
Chemical Physics Letters journal homepage: www.elsevier.com/locate/cplett
Dynamics of large nuclear-spin systems from low-order correlations in Liouville space Mark C. Butler, Jean-Nicolas Dumez, Lyndon Emsley * Université de Lyon (CNRS/ENS Lyon/UCB Lyon1), Centre de RMN à très hauts champs, 5 rue de la Doua, Villeurbanne, France
a r t i c l e
i n f o
Article history: Received 22 May 2009 In final form 3 July 2009 Available online 8 July 2009
a b s t r a c t Simulation of the coherent dynamics of a large lattice of strongly coupled nuclear spins is demonstrated. The evolution of N spins is simulated in a reduced Liouville space X which excludes coherences involving more than k spins, where k N is chosen to minimize the dimension of X without introducing unacceptable error. For polarization transfer within a spinning lattice of protons, comparison with exact simulations shows that k = 4 allows accurate simulation of the experimentally relevant powder-averaged curves. A simulation of polarization transfer among 144 protons is presented. Ó 2009 Elsevier B.V. All rights reserved.
1. Introduction It is astounding to consider that while sophisticated multidimensional nuclear magnetic resonance (NMR) methods are today routinely applied to solids, it is currently impossible to simulate exactly the dynamics of the coupled many-spin systems present in solids that generate the observables in these experiments. NMR methods play a central role in the characterization of solids, with applications ranging from studies of the atomic level structure and dynamics of biological macromolecules [1] to quantum computing [2], and in most of these areas, the inability to simulate the dynamics of a large lattice of coupled spins is a principal barrier to further progress [3,4]. Several approaches can be used to model such dynamics. Perturbation theory can be used to characterize the average effect of the network of couplings on experimental observables, yielding models in which the number of relevant variables is significantly reduced [5,6]. While such models allow simple analysis of some experiments, they require the use of empirical parameters and cannot account for effects that arise from the coherent dynamics of the lattice. Another approach to modelling involves numerical propagation of the full density matrix, using software that can simulate a wide range of NMR phenomena and methods [7–10]. Highly optimised algorithms [10,11], in combination with methods exploiting translational symmetry [4] in the solid state, have been used to perform exact simulations including up to 15 spins. Simulating systems of this size does yield insight into various fundamental phenomena [12–15] but is insufficient to capture large-scale effects, for which empirical models have to be added [5,12]. Simulations involving the full density matrix are limited by the fact that the dimension of Liouville space scales exponentially with * Corresponding author. Fax: +33 478896761. E-mail address: [email protected] (L. Emsley). 0009-2614/$ - see front matter Ó 2009 Elsevier B.V. All rights reserved. doi:10.1016/j.cplett.2009.07.017
sample size. For a system of N spins I = 1/2, for instance, Liouville space has dimension 4N. One method of circumventing this obstacle is to use a reduced Liouville space to simulate the coherent dynamics of large spin systems. An example of this approach was first introduced by Bruschweiler and Ernst [16], who used a ‘cogwheel’ model of spin motion in a solid to simulate the short-term dynamics of a chain of 10 spins. More recently, a scheme for reducing Liouville space for liquid-state NMR simulations has been applied to large spin systems [17,18]. In this scheme, small clusters of nearest-neighbor spins are defined, and spin i interacts directly with spin j in the simulation only if the two spins belong to the same cluster. However, in the case where each spin is strongly coupled to many other spins, e.g., a strongly coupled lattice of protons in a solid, the applicability of this approach has not been established. In this Letter we introduce a method of simulating large, densely coupled nuclear-spin systems in a reduced Liouville space. Using this approach we are able to simulate the medium-term coherent dynamics of a powder sample under magic-angle spinning for a crystalline lattice of 144 protons. This demonstration, which constitutes an order-of-magnitude increase in the number of spins for which such dynamics have been simulated, suggests that the structure and dynamics of solids may be probed by directly studying the dependence of experimental observables on coherent interactions involving many spins.
2. Simulation scheme 2.1. Motivation Exponential scaling of Liouville space with sample size can be understood as resulting from the vast number of ways in which correlations can develop within subcollections of particles. If the dimension of a 1-particle Liouville space is m, then a space of
378
M.C. Butler et al. / Chemical Physics Letters 477 (2009) 377–381
dimension mN is sufficient to characterize a system of N particles that are constrained to remain statistically independent, while the full Liouville space that includes all correlations between the particles has dimension mN. It is the requirement that all possible correlations between the particles be included in the computation that renders simulations of the full Liouville space of large systems prohibitively demanding. In general, however, only a small fraction of Liouville space is observed experimentally (either directly or indirectly), and the possibility exists that accurate simulations of large quantum systems could be performed in reduced Liouville spaces which include both the degrees of freedom monitored experimentally and some ‘surrounding space’ containing the degrees of freedom most strongly excited during the experiment. As an example, consider the problem of simulating spin diffusion (i.e., transfer of spin polarization) within a dipole-coupled lattice of protons. This problem is today central to the determination of structure in drug molecules [19–21] and solid proteins [22,23] by NMR crystallography. Since each spin i interacts strongly with a relatively large number of spins j, many distinct pathways for polarization transfer exist, and accurate simulation of spin diffusion requires a model in which the network of couplings is sufficiently dense to include the important transfer pathways. Note, however, that simultaneous transfer between these polarizations corresponds to motion along an even larger number of trajectories in Liouville space; it is conceivable that interference between the different trajectories in Liouville space could limit the net effect of high-order coherences, thereby enabling accurate simulation in a reduced space. An analogous effect was observed in multiple-quantum spectra [24] when phase interference among high-order coherences hindered the observation of their excitation [25]. Note as well that most solid-state NMR experiments involve magic-angle spinning (MAS) and are carried out on powder samples, and that these experimental conditions might shrink the region of Liouville space which effectively contributes to detected observables. During MAS, the modes (eigenoperators) of Liouville space are continually changing, and if the time scale of the sample rotation is sufficiently close to that of the spin dynamics, MAS could distribute the excitation which spreads through Liouville space among many modes whose contribution to the experimental signal could interfere destructively. Consistent with this hypothesis is the finding of Spiess and co-workers that sample rotation at frequencies of the order of the highest couplings often allows analysis of sideband intensities using only two-spin coherences [26]. Similarly, distinct crystallites in general have distinct sets of modes, and the evolution of an experimental observable for a given static crystallite can be written as a sum of oscillating terms associated with the modes of that crystallite. By summing contributions from distinct sets of modes, powder averaging could increase the efficiency of destructive interference among the modes excited under MAS.
the density operator while respecting the constraint that it remain within X (i.e., elements of the full Liouvillian that would couple operators within X to those lying outside of X could be considered formally set to zero). A first-order scheme was used to compute the evolution of the density matrix r under the action of LX. Given a basis fri g of the reduced space, the density operator can be expressed as a vector of coefficients si , where
r¼
X
s i ri :
The first-order change in r that occurs during a time step Dt is
Dr ¼ L X r Dt ¼
X ðsi DtÞLX ri :
ð1Þ
For simulations of polarization transfer within an N-proton system, each basis element ri was chosen as a product of N single-spin operators, with at most k < N of the single-spin operators different from the identity. Explicit formulas for the operators LX ri appearing in Eq. (1) were derived, and first-order numerical propagation of r was implemented by using these formulas to compute r þ Dr. Use of the derived formulas for LX ri eliminated the need to store LX as a matrix, and so the factor limiting the size of the system that could be simulated was the memory required to store r as a vector of coefficients, rather than the memory needed for an evolution matrix acting on that vector. Note as well that for fixed k, the memory required for this simulation scheme scales polynomially with N rather than exponentially. This improved scaling is a key feature for the simulation of large systems. For a solid sample containing protons in a field of a few Tesla or more, the dominant spin–spin interactions are governed by the secular homonuclear dipolar Hamiltonian:
HD ¼
X i
l0 c2 h ð1 3 cos2 hij Þð3Iiz Ijz Ii Ij Þ: 8pr3ij
Here Ii and Ij are nuclear-spin operators, rij is the internuclear distance, hij is the angle between the internuclear vector and the static magnetic field, and c is the gyromagnetic ratio. For simplicity, additional terms in the Hamiltonian were neglected in these simulations, although it would be straightforward to include them. Since solid-state NMR experiments are most often performed on powder samples using MAS, sample rotation and powder averaging were included in some simulations. Rotation of the sample was simulated by treating the Hamiltonian as time-dependent and updating it repeatedly during the numerical propagation of the density matrix. Powder averaging [27] was performed by summing observables over a set of simulations carried out with selected crystallites, using the ZCW scheme. Geometric structures were obtained from published organic crystal structures. For all simulations, the initial density matrix was set to r(0) = Iiz for a single proton i, and the evolution of hIjz(t)i was recorded for all protons j in the system.
2.2. Implementation
3. Results and discussion
We have tested these hypotheses by comparing exact simulations of polarization transfer within a lattice of N protons to simulations performed in a reduced Liouville space X that included only low-order correlations. All products of single-spin operators involving more than k correlated spins were excluded from X, where k < N is a small integer. Given an initial state defined by polarizations on certain spins, exact simulations of the evolution under a secular spin Hamiltonian can be performed within the zero-quantum subspace Z of the full Liouville space [16], and we defined reduced Liouville spaces by excluding from Z all coherences involving more than k correlated spins. The restriction LX of the full Liouvillian to a reduced space X was used to propagate
Exact simulations carried out with the SPINEVOLUTION software [10] on molecules containing up to 12 protons were compared with simulations performed in reduced Liouville spaces. Fig. 1 illustrates typical results obtained for a single orientation of a static molecule containing 12 protons, b-L-aspartyl-L-alanine. The simulations using only low-order correlations in Liouville space, which we dub LCL, reproduce the exact evolution for a very short initial time period, but after 100 ls, large errors develop in the simulations having k = 5 (dashed curves), i.e., in the simulations for which product operators involving more than 5 correlated spins were excluded. The simulations performed in a further reduced space, k = 4 (dotted curves), show pathological departure
379
M.C. Butler et al. / Chemical Physics Letters 477 (2009) 377–381
1 H H H H H H N O
12 spins, static
0.8
O
< Iiz > / || Iiz ||
O
H
H H
N H H O
O
H
0.6 0.4 0.2 0
-0.2 0
0.2
0.4 0.6 time (ms)
0.8
1
Fig. 1. Simulations of polarization transfer for a single molecule of b-L-aspartyl-Lalanine (CSD entry: fumtem). Exact simulations (solid) of this system of 12 protons are compared to simulations performed in reduced Liouville spaces that exclude coherences involving more than k = 4 spins (dotted) or k = 5 spins (dashed). The initial density matrix is r(0) = ICH2,z, and the observables hICH2,zi (red) and hICHasp,zi (blue) are plotted. The time step in the simulation is 1 ls. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
from the exact simulations (solid curves), and from the curves corresponding to k = 5. The relatively simple oscillations in the k = 4 curves are consistent with the hypothesis that the initial state of the spin system is dominated by the contribution from a few spurious modes in the corresponding reduced space. Similar short-term accuracy was reported using a cogwheel model to propagate the density matrix in a reduced Liouville space that excluded product operators involving more than 3 correlated spins [16]. The error present in these LCL simulations shows that the excluded regions of Liouville space are relevant for simulations of a single static crystal. The accuracy of the LCL simulations is significantly improved by adding spinning to the simulation, as shown in Fig. 2. While in the static case, agreement was obtained for about 50 ls, the agreement with exact simulations (solid curves) is pre-
served under 10 kHz MAS for longer time periods (around 250 ls), and the decay of the initial polarization, although not reproduced exactly, is qualitatively correct for both k = 4 (dotted curves) and k = 5 (dashed curves). The k = 4 curves appear to include contributions from many interfering modes, rather than showing the simple oscillations associated with the evolution of a small number of modes. The accuracy of the LCL simulations is further improved by the combination of spinning and powder averaging, as shown in Fig. 3, in which the exact curves are now quantitatively reproduced using the LCL scheme for times up to 1000 ls, corresponding to the whole of the experimentally relevant timescale. The Liouville space that includes 4-spin correlations (k = 4) is now sufficient to reproduce detailed features of the spin diffusion curves. LCL simulations were compared with exact simulations performed for a collection of 25 10-proton organic molecules and 15 12-proton organic molecules randomly selected from the Cambridge Structural Database (CSD) [28]. For a system containing N protons, an exact simulation was performed for each of the N distinct initial states r(0) = Iiz, and each of these simulations yielded N curves hIjz ðtÞi=kIjz k; i.e., 4660 exact curves were compared against LCL curves. The error in the LCL simulations was characterized by calculating a root-mean square (rms) error for each curve, and the distribution of rms errors for the set of 10-proton curves and the set of 12-proton curves were found to be to be similar. Among the k = 4 and k = 5 curves, respectively, 78% and 95% had an rms error of 0.01 or less, as shown in Fig. 4. (The Supplementary material includes all curves used for this test of the LCL simulation scheme.) These results support the hypothesis that for spinning organic powders, the dynamics of spin diffusion among protons can be accurately modelled using a small fraction of the full Liouville space. In the case of a 12-proton system, for instance, Table 1 shows that the dimension of the k = 4 reduced Liouville space is less than 0.1% as large as that of the full Liouville space. An additional finding is that reduced Liouville spaces X adequate for simulating these dynamics can be obtained simply by excluding all products of single-spin operators involving more than k = 4 or k = 5 correlated spins. Note that for simulations performed in such a space, excitation that develops in k-spin coherences can only be
1 1 12 spins, powder, MAS
0.8
12 spins, MAS < Iiz > / || Iiz ||
< Iiz > / || Iiz ||
0.8 0.6 0.4
0.4 0.2
0.2 0
0.6
0 0
0.2
0.4 0.6 time (ms)
0.8
1
Fig. 2. Simulations of polarization transfer for a single molecule of b-L-aspartyl-Lalanine under 10 kHz MAS. Exact simulations (solid) of the 12-proton system are compared to simulations performed in reduced Liouville spaces that exclude coherences involving more than k = 4 spins (dotted) or k = 5 spins (dashed). The initial state, observables, and time step are the same as in Fig. 1.
0
0.2
0.4 0.6 time (ms)
0.8
1
Fig. 3. Powder-averaged simulations of polarization transfer for a single molecule of b-L-aspartyl-L-alanine under 10 kHz MAS. Exact simulations (solid) of the 12proton system are compared to simulations performed in reduced Liouville spaces that exclude coherences involving more than k = 4 spins (dotted) or k = 5 spins (dashed). The initial state, observables, and time step are the same as in Fig. 1, and a ZCW50 set of 50 orientations is used for the powder average.
380
M.C. Butler et al. / Chemical Physics Letters 477 (2009) 377–381
4000
number of curves
k=4
k=5
2000
0
0.02
0.04
0.06
0.08 0 rms error
0.02
0.04
0.06
0.08
Fig. 4. The rms error for 4660 polarization curves obtained from 25 10-proton systems and 15 12-proton systems. For each curve, exact simulations were compared with simulations performed in reduced Liouville spaces that excluded coherences involving more than k = 4 or k = 5 spins.
Table 1 The dimensions of Liouville spaces that can be used for simulation of polarization transfer within 10-proton and 12-proton systems. The full Liouville space includes all possible spin correlations among the N protons. Given an initial state defined by polarizations on certain protons, exact simulations of the evolution under a secular spin Hamiltonian can be performed within the zero-quantum subspace. The reduced Liouville spaces are subspaces of the zero-quantum subspace, formed by excluding product operators involving more than k correlated spins. Full Liouville space
Zero-quantum subspace
Reduced spaces k=4
10 protons
1:0 10
12 protons
1:7 107
5
6
1:8 10
5:0 10
1:8 104
2:7 106
1:1 104
5:2 104
1 144 spins, powder, MAS
< Iiz > / || Iiz ||
0.8 0.6 0.4 0.2 0 0
0.2
0.4 0.6 time (ms)
0.8
k=5 3
1
Fig. 5. Powder-averaged simulations of polarization transfer in a system of 144 protons (solid) and a system of 12 protons (dotted) under 10 kHz MAS. The simulations are performed in a reduced Liouville space that excludes coherences involving more than k = 4 spins. The 12-proton system is a single molecule of b-Laspartyl-L-alanine (CSD entry: fumtem), while the 144-proton system is a rectangular block of three primitive unit cells of the same molecule. For the system of 144 spins, the minimum image convention [29] is used to define periodic boundary conditions, and the observables are summed over crystallographically equivalent spins. The initial density matrix is r(0) = ICH2,z, and the observables hICH2,zi (red) and hICHasp,zi (blue) are plotted. The time step in the simulation is 1 ls, and a ZCW50 set of 50 orientations is used for the powder average. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
reflected back to lower-order coherences, rather than propagated outward to coherences involving additional spins. We are investigating the way in which such errors introduced into the dynamics by the reduction scheme do not contribute significantly to the spin diffusion curves. A key feature of the LCL simulation scheme is that its favourable scaling properties allow simulation of large systems. As a demonstration, we have simulated powder-averaged spin diffusion curves for a spinning lattice of 144 protons, i.e., the 12 molecules contained in three primitive unit cells of b-L-aspartyl-L-alanine. Although a detailed comparison of the 12-proton simulation and the 144-proton simulation is outside the scope of this Letter, we note that Fig. 5 shows the expected qualitative similarity between the curves for the two systems, although quantitative differences are present. The full Liouville space for the 144-proton system has dimension 5 1086, which makes the full density matrix impossible to simulate. The reduced Liouville space has dimension 3 108, and the LCL simulation for a single orientation used 1 GB of random-access memory and required only 8 h of computing time on a single 2.8 GHz processor. 4. Conclusion We have found that exact powder-averaged simulations of polarization transfer within spinning organic crystals can be quantitatively reproduced over the full experimentally relevant timescale using reduced Liouville spaces that include only low-order correlations. For 10-proton and 12-proton organic structures, exact simulations of polarization transfer under the secular homonuclear dipolar Hamiltonian were quantitatively reproduced by simula-
M.C. Butler et al. / Chemical Physics Letters 477 (2009) 377–381
tions performed using reduced basis sets that excluded product operators involving more than 4 correlated spins. For the organic structures tested, the dimensions of the reduced Liouville spaces are orders of magnitude smaller than the dimension of the full Liouville spaces; however, the combination of magic-angle spinning and powder averaging enables accurate simulation of proton spin diffusion within the reduced spaces. Since the dimensions of the reduced space scale polynomially rather than exponentially with the number of spins, the simulation of large crystalline lattices is feasible, as demonstrated by a simulation of polarization transfer among 144 protons. This demonstration suggests that simulations including only low-order correlations in reduced Liouville spaces may enable new studies of the structure and dynamics of solids by revealing the dependence of experimental observables on coherent interactions involving many spins. Appendix A. Supplementary material Supplementary data associated with this article can be found, in the online version, at doi:10.1016/j.cplett.2009.07.017. References [1] A.E. McDermott, Curr. Opin. Struc. Biol. 14 (2004) 554. [2] L.M.K. Vandersypen, I.L. Chuang, Rev. Mod. Phys. 76 (2004) 1037. [3] W.X. Zhang, N. Konstantinidis, K.A. Al-Hassanieh, V.V. Dobrovitski, J. Phys.: Condens. Matter 19 (2007).
381
[4] P. Hodgkinson, D. Sakellariou, L. Emsley, Chem. Phys. Lett. 326 (2000) 515. [5] B.H. Meier, Adv. Opt. Magn. Reson. 18 (1994). [6] M. Ernst, B.H. Meier, in: I. Ando, T. Asakura (Eds.), Solid State NMR of Polymers, Elsevier, Amsterdam, 1998, p. 83. [7] M. Bak, J.T. Rasmussen, N.C. Nielsen, J. Magn. Reson. 147 (2000) 296. [8] F.S. Debouregas, J.S. Waugh, J. Magn. Reson. 96 (1992) 280. [9] P. Hodgkinson, L. Emsley, Prog. Nucl. Magn. Reson. Spectrosc. 36 (2000) 201. [10] M. Veshtort, R.G. Griffin, J. Magn. Reson. 178 (2006) 248. [11] V.E. Zorin, M. Ernst, S.P. Brown, P. Hodgkinson, J. Magn. Reson. 192 (2008) 183. [12] D. Sakellariou, P. Hodgkinson, L. Emsley, Chem. Phys. Lett. 293 (1998) 110. [13] V.E. Zorin, S.P. Brown, P. Hodgkinson, J. Chem. Phys. 125 (2006) 144508. [14] V.E. Zorin, B. Elena, A. Lesage, L. Emsley, P. Hodgkinson, Magn. Reson. Chem. 45 (2007) S93. [15] J.S. Waugh, Mol. Phys. 95 (1998) 731. [16] R. Bruschweiler, R.R. Ernst, J. Magn. Reson. 124 (1997) 122. [17] I. Kuprov, N. Wagner-Rundell, P.J. Hore, J. Magn. Reson. 189 (2007) 241. [18] I. Kuprov, J. Magn. Reson. 195 (2008) 45. [19] B. Elena, L. Emsley, J. Am. Chem. Soc. 127 (2005) 9140. [20] E. Salager, R.S. Stein, C.J. Pickard, B. Elena, L. Emsley, Phys. Chem. Chem. Phys. 11 (2009) 2610. [21] B. Elena, G. Pintacuda, N. Mifsud, L. Emsley, J. Am. Chem. Soc. 128 (2006) 9555. [22] F. Castellani, B. van Rossum, A. Diehl, M. Schubert, K. Rehbein, H. Oschkinat, Nature 420 (2002) 98. [23] C. Wasmer, A. Lange, H. Van Melckebeke, A.B. Siemer, R. Riek, B.H. Meier, Science 319 (2008) 1523. [24] M. Munowitz, A. Pines, Adv. Chem. Phys. 66 (1987) 1. [25] Y.S. Yen, A. Pines, J. Chem. Phys. 78 (1983) 3579. [26] C. Filip, S. Hafner, I. Schnell, D.E. Demco, H.W. Spiess, J. Chem. Phys. 110 (1999) 423. [27] M. Eden, Concepts Mag. Reson. A 18A (2003) 24. [28] F.H. Allen, Acta Crystallogr. Sect. B-Struct. Sci. 58 (2002) 380. [29] M.P. Allen, D.J. Tildesley, Computer Simulation of Liquids, Oxford Science, Oxford, 1987.