Exact NMR simulation of proteinsize spin systems using tensor train formalism
Abstract
We introduce a new method, based on alternating optimization, for compact representation of spin Hamiltonians and solution of linear systems of algebraic equations in the tensor train format. We demonstrate the method’s utility by simulating, without approximations, a \ce^15N NMR spectrum of ubiquitin — a protein containing several hundred interacting nuclear spins. Existing simulation algorithms for the spin system and the NMR experiment in question either require significant approximations or scale exponentially with the spin system size. We compare the proposed method to the Spinach package that uses heuristic restricted state space techniques to achieve polynomial complexity scaling. When the spin system topology is close to a linear chain (e.g. for the backbone of a protein), the tensor train representation is more compact and can be computed faster than the sparse representation using restricted state spaces.
Keywords: density matrix renormalization group, alternating minimal energy, tensor train, nuclear magnetic resonance, protein
1 Introduction
The amount of patience required to simulate exactly a nuclear magnetic resonance (NMR) spectrum of an spin system scales approximately as . That much is rarely available, and considerable thought has consequently been given over the last decade to more efficient methods [1, 8, 19, 28, 70], particularly those that promise to achieve that objective in polynomial time. Such algorithms do exist [8, 31], but they make significant a priori assumptions about the spin system evolution — it is usually assumed that the system stays weakly correlated for the duration of the experiment [31, 42].
Outside the NMR community, significant progress was recently made with the development of tensor structured methods [61, 39, 62, 36, 25], all of which descend broadly from the density matrix renormalization group (DMRG) [74, 75] as well as matrix product state (MPS) [23, 38] and matrix product operator (MPO) [69] formalisms. Typical applications of DMRG in condensed matter theory are 1D spin chains [61, 7, 65, 73] with recent extensions to 2D lattices [44, 46, 80, 66]. DMRG has also been put to good use in electronic [77, 11, 81, 43, 78, 63, 79] and nuclear [53, 58] structure theory, but magnetic resonance spectroscopy has so far received little attention — the spin systems encountered in the daily practice of NMR and EPR (proteins, radicals, polynucleotides, polysaccharides) are irregular threedimensional roomtemperature networks with multiple interlocking loops in the spin coupling graph and no identical couplings [10]. When the strict requirement for correct wavefunction phase during the very long (milliseconds to seconds) dissipative spin system trajectories is added to the list, timedomain DMRG methods are currently struggling.
There are some biologically relevant cases, however, that may still be treated as linear chains — for the purposes of simulating simple backbone NMR experiments, protein side chains may often be ignored. This makes the corresponding spin system a weakly branched linear chain that is amenable to DMRG type treatment. Simple NMR experiments can also be reformulated as a matrixinversetimesvector problem in the frequency domain, for which efficient algorithms in tensor product formats have recently emerged [30, 15, 16, 17]. We report in this communication the behavior of the AMEn algorithm [16, 17, 18], applied to the solution of the NMR simulation problem in the frequency domain, as well as to the technical task of adding together, without loss of accuracy, tensor train representations of thousands of spin Hamiltonian terms for a protein.
Having integrated the algorithms described below into Spinach (a largescale magnetic resonance simulation library [28]), we are reporting here the first exact quantum mechanical simulation of a liquidstate 1D NMR spectrum for a protein backbone spin system with several hundred coupled spins. Beyond the physical assumptions made by chemists at the problem formulation stage and the controllable numerical rounding error of the tensor train format itself [50], there are no approximations.
2 Tensor product formats
Tensor product expressions appear naturally in spin dynamics because the state space of a multispin system is a direct product of state spaces of individual spins [22]. A simple example is the nuclear Zeeman interaction Hamiltonian
(1) 
where is the number of spins, is the applied magnetic field, are nuclear chemical shielding tensors, and the sum runs over all nuclei. Cartesian components of nuclear spin operators have the following tensor product form
(2) 
where denotes an identity matrix of appropriate dimension and Pauli matrices occur at the th position in the tensor product sequence. This representation is known in numerical linear algebra as the canonical polyadic (CP) format [39]. Although CP representations have been known in magnetic resonance spectroscopy for a long time [6], they suffer in practice from rapid inflation — spin Hamiltonians encountered in NMR and ESR (electron spin resonance) systems can be complicated [22] and, even for simple initial conditions, the number of terms in the canonical decomposition increases rapidly during system evolution. More ominously, the number of CP terms can change dramatically after small perturbations of the Hamiltonian or the system state. A simple example is
(3) 
where The left hand side of this equation contains direct product terms, given by Eq. (2), but the expression approximating it on the right hand side has only two direct product terms, and one could be tempted to use it to reduce storage and CPU time. However, both terms of the approximation grow to infinity when and the accuracy is lost due to rounding errors. Such instabilities in the CP format make it difficult to use — in finite precision arithmetic the number of terms in the decomposition quickly becomes equal to the dimension of the full state space and any efficiency savings disappear.
Unlike the CP format, which is an open tensor network, closed tensor network formats are stable to small perturbations. The most popular closed tensor network format was repeatedly rediscovered and is currently known under three different names: DMRG in condensedmatter physics [74, 75], MPS [23, 38, 68, 52]/MPO [69, 12] in computational physics, and TT (tensor train) in numerical linear algebra [50]. A tensor train is defined, using the standard notation of numerical linear algebra [13], as follows:
(4) 
The TT representation of the total operator in Eq. (3) is similar to the high–dimensional Laplacian [35]:
(5) 
with and The number of terms in each summation (known as bond dimension, or TT rank) is two, and the number of entries of the decomposition is now bounded. The TT representation of in Eq. (5) has single–spin operators, each of which is either zero, or identity , or the Pauli matrix The CP representation of in Eq. (3) has such operators — the tensor train representation is clearly more memory–efficient.
Another notable example is the ZZ coupling Hamiltonian that often makes an appearance in models of simple linear spin chains:
(6) 
As written, this is a CP format with terms and singlespin operators entering direct products. The corresponding TT representation is
(7) 
Here each summation runs over three terms only, and the total number of singlespin operator matrices appearing in is , much fewer than the one of the CP format in Eq. (6).
Storage requirements of tensor structured representations (both CP and TT) stand in sharp contrast with the classical approach to magnetic resonance simulations [1, 70], where the Hamiltonian is represented as a sparse matrix with all nonzero entries stored in memory. As soon as the matrix is assembled, CPU and memory resources grow exponentially with the number of spins making the simulation prohibitively difficult for large systems. Tensor structured methods avoid this problem (it is known colloquially as the curse of dimensionality) by keeping all data in compressed formats of the form given in Eqs. (1) and (4) and manipulating it without ever opening up the Kronecker products.
A very considerable body of literature exists on manipulating expressions directly in tensor product formats [39, 36, 25]. In particular, a given matrix may be converted into the TT format using sequential singular value decompositions [61, 52, 50]. Given tensors in the TT format, one can perform linear or bilinear operations (addition, elementwise multiplication, matrixvector multiplication) [62, 50], Fourier transform [5, 14], and convolution [34] directly in the TT format, avoiding exponentially large arrays and computational costs.
These developments would have permitted large–scale magnetic resonance simulations entirely in the TT format, were it not for a significant obstacle — the summation operation in tensor train representations is an expensive procedure that carries a significant accuracy penalty due to the need to recompress the representation to keep the bond dimensions low. Spin Hamiltonians of practically interesting biological systems contain many thousands one– and two–spin terms of the kind shown in Eq. (2). Intermediate expressions in spin dynamics simulations also frequently involve large sums. We demonstrate below that in those circumstances the standard bundleandrecompress tensor network summation procedure leads either to the bond dimension expansion beyond the limits of modern computing hardware, or to a catastrophic accuracy loss. This problem also occurs with three–dimensional potentials encountered in electronic structure theory [37, 59, 24]. Here we propose an alternative algorithm for computing large sums, based on alternating tensor train optimization, and use it to enable NMR simulations on proteinsize spin systems.
3 Simulation setting and experimental context
Fully \ce^13C and \ce^15N–labelled protein human ubiquitin (PDB code 1D3Z, Figure 1) containing over a thousand magnetic nuclei in 76 amino acid residues was chosen for testing purposes with two types of spin subsystem selection: backbone (H, N, C, CA, HA) and extended backbone (H, N, C, CA, CB, HA, HB). Both cases involve a weakly branched continuous chain of spinspin couplings and are encountered in the simulation of a large class of protein backbone NMR experiments that map out the protein bonding network and thereby assist in molecular structure determination: HNCO [33], HNCOCA [3], HNCA [33], and HSQC [60]. The isotropic NMR Hamiltonian was assembled using chemical shift values from the BMRB database [67] and couplings from the literature data [9, 45, 51, 71, 72]. In the cases where an experimental value of a particular coupling was not available in the literature, it was estimated based on the known values for structurally similar substances [2, 27, 40] — for most NMR simulation purposes and certainly for the purpose of the demonstration of the performance of the tensor train algorithm the accuracy of such coupling estimates (about ) is sufficient. The raw data for the magnetic couplings used in this work is available in the example set supplied with the current public version of the Spinach library [28].
NMR experiments were performed at on a Varian Inova MHz ( Tesla) spectrometer equipped with a Z–gradient triple–resonance cryogenic probe using a mM sample of uniformly \ce^13C – and \ce^15N–labelled human ubiquitin in \ceD_2O. \ce^15N spectra were collected as 2D \ce^1H–\ce^15N HSQC [4] spectra incorporating gradient enhanced coherence selection [32] and water flipback. The spectra were recorded with acquisition times of (, \ce^15N) and (, \ce^1H). During the \ce^15N evolution period, \ce^1J_HN and \ce^1,2J_NC couplings were either allowed to evolve, or decoupled by insertion of a rectangular \ce^15N or a shaped \ce^13C inversion pulse using the central lobe of the sinc function. During \ce^1H acquisition \ce^15N nuclei were either evolved or decoupled using ppm broadband WURST sequence [41].
The liquid state NMR Hamiltonian of \ce^13C,\ce^15N–labelled ubiquitin is:
(8) 
where canonical NMR spectroscopy notation is used [22], index runs over all nuclei, and indices run over pairs of nuclei that belong to the same isotope, and run over pairs of nuclei that belong to different isotopes, and run over the nuclei influenced by radiofrequency pulses, and are time profiles of those pulses, are offset frequencies arising from the chemical shielding of the corresponding nuclei [56], are “strong” NMR couplings [57], are “weak” NMR couplings [22], and spin operators are defined by Eq. (2). In the case of extended ubiquitin backbone, the Hamiltonian in Eq. (8) contains shielding terms, coupling terms, and radiofrequency terms. All calculations reported below were performed by extending the functionality of Spinach library [28] to the tensor train formalism and interfacing it to TTToolbox [49] where appropriate.
Due to the abundance of complicated multipulse NMR experiments with timedependent Hamiltonians [22], magnetic resonance simulations are generally carried out in the time domain. They always require longterm evolution trajectories with accurate phases (at least ms, much longer than the reciprocal Hamiltonian norm) for the density operator under the Liouville–von Neumann equation:
(9) 
where is the relaxation superoperator ( model with literature values for relaxation times [10] was used in the present work), is the thermal equilibrium state, and is the observable operator, usually a sum of or operators on the spins of interest. In very simple cases where the Hamiltonian is not timedependent, the general solution to Eq. (9) can be written as:
(10) 
where is the Hamiltonian commutation superoperator.
Direct time domain evaluation of this equation in tensor train format, either using explicit operator exponentiation or Krylov type propagation techniques, does not appear to be possible — in all cases described by Eq. (8) the ranks in the tensor train expansion quickly grow beyond the capacity of modern computers. Increasing the singular value cutoff threshold at the representation compression stage leads to catastrophic loss of accuracy. Fortunately, there are simple cases (most notably pulseacquire 1D NMR spectroscopy) where amplitudes at only a few specific frequencies are actually required for the Fourier transform of Eq. (10), meaning that the problem can be reformulated in the frequency domain:
(11) 
That is, to compute the observable at the point in the frequency domain, we need to solve a linear system The problem formulation in Eq. (11) sacrifices a great deal of generality compared to Eq. (9) (simulation of arbitrary NMR pulse sequences is no longer possible), but it does serve as a stepping stone and enables the demonstration calculation presented below.
4 Tensor train algorithm for the summation and solution of linear systems
The DMRG algorithm was initially proposed [74, 75] to find the ground state of a Hermitian matrix by the minimization of the Rayleigh quotient The dynamical DMRG algorithm [30] was then developed to find the solution of a linear system with a Hermitian positive definite matrix by the minimization of the energy function Apart from the change of the minimization target function, the two algorithms are similar.
In DMRG formalism the solution is sought in the form of a tensor train introduced in Eq. (4), but the minimization over all cores simultaneously is a complicated nonlinear problem. To make the procedure feasible, it is replaced by a sequence of optimizations carried over one core at a time:
(12) 
The TT format is linear in all cores . This fact may be expressed as where the frame matrix maps the parameters of the TT core to the vector . The linearity allows to rewrite Eq. (12) as where is the energy function for the local problem with and . Using the nonuniqueness of the tensor train representation (4), one can always construct the representation with the unitary frame matrix , that guarantees the stability of the local problem. Such a choice is known as gauge condition in the MPS literature, and canonical form in the DMRG literature. After the solution is computed, we substitute in the tensor train, and continue for and then back and forth along the chain.
The convergence of the above described onesite DMRG procedure depends on the initial guess and in particular on the initial choice of the TT ranks because they remain the same during the sequence of updates defined by Eq. (12). This is a severe restriction and additional measures are therefore taken to adapt the TT ranks during the computations. One way to do that is to replace the optimization over single cores by the optimization over pairs of neighboring cores, and then to adapt the TT rank between them. Another possibility is to expand the search space by adding auxiliary directions. The first method of the latter type is the corrected onesite DMRG algorithm [76], which targets in addition to a surrogate of the next Krylov vector
For the solution of linear systems, the alternating minimum energy (AMEn) algorithm was recently proposed [16, 17], which also uses an additional direction to adapt tensor train ranks. The local optimization step in AMEn is carried over one site only. To adapt TT ranks and improve convergence, TT blocks are expanded by auxiliary information, The enrichment introduces new directions in the subspace spanned by . A good choice of the enrichment is the component of the TT representation (exact or approximate) of the residual AMEn algorithm is as fast as onesite methods, but as rank adaptive as the twosite DMRG algorithm, and demonstrates comparable or better convergence rates. For the solution of a linear system with a Hermitian positive definite matrix, it has a proven global bound on the geometrical convergence rate. Unlike the corrected onesite DMRG method [76], the AMEn algorithm is stable to perturbations and free from tuning parameters and heuristics [18]. The rank adaptation strategy in the enrichment phase of AMEn is determined by a single relative accuracy parameter.
In this work we use the AMEn algorithm for two purposes. First, we apply it to a system with a trivial matrix but a complicated righthand side , which is a sum of many elementary tensors like the one in Eq. (2). This allows us to compress a Hamiltonian returned by the Spinach package from the CP format given by Eq. (8) into the TT format Eq. (4). The Hamiltonian is stretched into a vector, and the target functional is a Frobeniusnorm distance between a given Hamiltonian and Hamiltonian sought in the tensor train format. The onesite optimization in Eq. (12) is effectively the solution of the overdetermined linear system using the least squares method. For the unitary frame matrix we have and therefore the local optimization step is obtained by contracting the frame matrix with the given Hamiltonian The enrichment step uses a lowrank approximation of the error which is obtained by onesite DMRG optimization.
After the Hamiltonian is compressed into the tensor train format, we compute 1D NMR spectra by solving the linear system in Eq. (11). Since the matrix is not expected to be Hermitian positive definite, we consider instead an equivalent symmetrized problem For demonstration purposes, we chose a simple nonselective damping relaxation model and the same operator for the initial and the detection state, where are the total spin operators of all \ce^15N nuclei in the system. This avoids explicit radiofrequency pulses and makes the Hamiltonian in Eq. (8) timeindependent and realvalued, those properties are also inherited by the commutation superoperator Since the detection state is also realvalued, the NMR spectrum in Eq. (11) can be computed from that we obtain as follows:
(13) 
This equation is solved by the AMEn algorithm at each point in the userspecified frequency interval.
5 Results
As discussed above, a major problem in the application of tensor train methods to magnetic resonance simulation of large systems is the calculation of lengthy sums involved in the construction of spin Hamiltonians and density matrices, and their compression into the TT format. Fig. 2 illustrates the performance of our proposed solution to this problem in the case of minimal (H, N, C, CA, HA) and extended (H, N, C, CA, HA, CB, HB) ubiquitin backbone spin systems. Storage requirements for the TT format in Eq. (4) depend on all TT ranks (bond dimensions) and are characterized by the effective TT rank defined by It is clear from the left panels of Fig. 2 that the primary showstopper — rapid growth in the tensor train rank — has been removed by the AMEn method: the effective ranks stay below for the extended backbone and below for the minimal backbone, well within the capability of modern desktop workstations. Since is smaller than the number of terms in the CP representation, the TT format with operators provides more compact storage than the CP format.
The alternative to AMEn is binary summation, which adds up Hamiltonian terms pairwise and recompresses the representation after each addition. As demonstrated in Fig. 2, binary summation drives tensor train ranks up to several hundred and thereby makes the solution of the linear system in Eq. (13) exceedingly difficult. It is clear from the right panels of Fig. 2 that the CPU time requirements of AMEn summation compared to binary summation are essentially the same, making AMEn procedure clearly superior for all practical purposes. The resulting representation of the ubiquitin backbone spin Hamiltonian matrix is, up to the rounding error of the complex double precision arithmetic, exact. In magnetic resonance spectroscopy this is an unprecedented development — ubiquitin NMR simulation is currently just about feasible [21], with significant approximations and colossal computational resources. Tensor train representation is therefore a large step forward, even though Eq. (11) is not in general applicable to arbitrary NMR experiments.
After the Hamiltonian is compressed, we compute \ce^15N pulseacquire NMR spectra using Eq. (13) with the AMEn algorithm and compare it to the simulation produced by the restricted state space (RSS) approximation [28], which is currently the only other method that is capable of handling NMR systems of this size. As demonstrated in Fig. 3 (top), When the basis set used by RSS is increased, its result converges to the one produced by AMEn, and the relative deviation between two methods falls below across the frequency interval.
It is instructive to compare the results of AMEn simulations with those produced by the dynamical DMRG [30] technique. As shown in Fig. 3 (bottom), the NMR spectrum computed by AMEn matches the reference spectrum returned by RSS with only minor deviations, while the accuracy of the result computed by the dynamical DMRG algorithm at the same relative accuracy parameter is unacceptable. DMRG does of course produce the right answer if a much tighter accuracy parameter is specified, but the simulation time goes up by several orders of magnitude. AMEn does therefore appear to have a better accuracytoeffort ratio. This is also confirmed by the convergence graph of AMEn and DMRG, given in the same figure, where the relative deviation between the computed and the reference values is shown during the iterations (sweeps) for both DMRG and AMEn. Note also that the inexact values of the spectrum, computed by AMEn and DMRG are always below the reference values; this was first noted by Jeckelmann [30]. The comparison in Fig. 3 is made using to visually emphasize the observed difference between the two methods; the same conclusion also holds for more accurate calculations using
Due to the intrinsically low sensitivity of liquid state \ce^15N protein NMR spectroscopy, it is not possible to record the experimental equivalent of Fig. 3 directly with a sufficient signaltonoise ratio; we have therefore taken a somewhat longer route to the experimental validation of the tensor train simulation — Fig. 4 shows experimental protondetected \ce^15N–\ce^1H HSQC spectra of ubiquitin, compared to the simulations obtained at the basis set limit of the RSS formalism [31]. Perfect agreement is apparent in both cases. This provides an experimental evidence to the accuracy of the restricted state space method. The tensor train results in Fig. 3 can now be justified by comparison to the RSS results — it is clear that the TT formalism performs as intended.
6 Discussion
The successful 1D NMR simulation notwithstanding, very significant obstacles remain on the path to practical applications of the tensor train formalism to NMR spectroscopy. The following issues should be addressed in future work to fully uncover the potential of the DMRG/MPS/TT formalism for spin dynamics simulations:
(a) The requirement for the spin system to be a chain or a tree should be lifted. Biological magnetic resonance spin systems are irregular polycyclic interaction networks with multiple interlocking loops in the coupling graph, particularly in solid state NMR, where internuclear dipolar couplings form very dense meshes. A generalization of tensor train algorithms to general contraction networks that fully mimic the molecular structure is therefore required.
(b) Rank explosion problem for timedomain simulations should be solved. It is clear from the success of the restricted state space approximation [42, 28] that the order of spin correlation in many evolving magnetic resonance spin systems either is or may safely be assumed to be quite low. This suggests data sparsity and separability, and indicates that some kind of lowrank decomposition is possible. One likely direction is through the enforcement of symmetries and conservation laws within the tensor train format itself during time evolution.
(c) Our experience indicates that tensor train objects are very far from being dropin replacements for their matrix counterparts in standard simulation algorithms and software — it does actually appear that nearly everything in the very considerable body of magnetic resonance simulation methods needs to be adapted to the realities of DMRG. Current implementation of tensor product methods still requires a number of tuning parameters (approximation accuracies, TT ranks of the enrichment, etc.). Broad adoption of tensor network algorithms would require basic linear algebra operations to be handled transparently and seamlessly by the existing simulation software packages, in the same way as sparse matrices currently are.
(d) Transparent and clear tensor train approximation accuracy criteria, rank control and a priori error bounds should be developed in order to estimate the influence of the representation compression errors on the accuracy of the final result. This problem is particularly acute for the state vector phase in time domain simulations: magnetic resonance experiments rely critically on the phase being correctly predicted.
All of that having been said, we are very optimistic about the future of lowrank tensor product DMRG/MPS/TT methods, having also found them useful in Fokker–Planck type formalisms related to NMR and EPR spectroscopy [20]. Their primary strength is the lack of heuristic assumptions and the controllable nature of the representation accuracy. An experimental implementation of tensor train magnetic resonance simulation paths, via an interface to the TTToolbox [49], is available in version 1.3.1980 of our Spinach library [28].
7 Conclusions and future work
Even with their welldocumented limitations (the requirement for the spin system to be close to a chain, difficulty with longrange timedomain simulations, code implementation challenges, etc.), the ability of tensor network formalisms to simulate simple liquid state NMR spectra of large spin systems essentially without approximations is impressive. They cannot yet match the highly optimized dedicated methods developed by the magnetic resonance community [21], but if some of the limitations are lifted by the subsequent research, DMRG methods would have the potential to become a very useful formalism in NMR research.
Having solved in this paper the last purely technical problem on the way to the broad adoption of tensor train formalism in magnetic resonance spectroscopy, we are quite optimistic about its potential. In particular, the following avenues appear promising:

Development of reliable tensor train methods for solving linear systems of algebraic equations with indefinite matrices, and time evolution problems.
Elsewhere in magnetic resonance, benefits to electron spin resonance spectroscopy, with its starshaped spin interactions graphs, are likely to be harder to achieve, but may still be obtained by exploiting the direct product structure of combined spin and spatial dynamics appearing in Fokker–Planck type problems [20].
Acknowledgements
We are grateful to Garnet K.L. Chan for patiently explaining the DMRG formalism to IK during his visit to Cornell University and to Zenawi T. Welderufael for finding some of the less obvious ubiquitin couplings in the literature. We acknowledge the Iridis High Performance Computing Facility, and the associated support services at the University of Southampton. JMW would like to acknowledge the Wellcome Trust for support of the Southampton NMR centre and we would like to thank the Geoff Kelly at NIMR (Mill Hill) for lending us the ubiquitin sample. The project is supported by EPSRC (EP/H003789/2, EP/J013080/1).
References
 [1] M. Bak, J. T. Rasmussen, and N. C. Nielsen, SIMPSON: A general simulation program for solidstate NMR spectroscopy, J Magn. Reson., 147 (2000), pp. 296–330.
 [2] M. Barfield, R. J. Spear, and S. Sternhell, Allylic interproton spin–spin coupling, Chem. Rev., 76 (1976), pp. 593–624.
 [3] A. Bax and M. Ikura, An efficient 3D NMR technique for correlating the proton and backbone amide resonances with the carbon of the preceding residue in uniformly / enriched proteins, J. Biomolecular NMR, 1 (1991), pp. 99–104.
 [4] G. Bodenhausen and D. J. Ruben, Natural abundance nitrogen15 NMR by enhanced heteronuclear spectroscopy, Chem. Phys. Lett., 69 (1980), pp. 185–189.
 [5] D. E. Browne, Efficient classical simulation of the quantum Fourier transform, New J. Phys., 9 (2007), p. 146.
 [6] R. Brüschweiler and R. R. Ernst, A cogwheel model for nuclearspin propagation in solids, J. Magn. Reson., 124 (1997), pp. 122–126.
 [7] R. Bursill, T. Xiang, and G. Gehring, The density matrix renormalization group for a quantum spin chain at nonzero temperature, Journal of Physics: Condensed Matter, 8 (1996), p. L583.
 [8] M. C. Butler, J. N. Dumez, and L. Emsley, Dynamics of large nuclearspin systems from loworder correlations in Liouville space, Chem. Phys. Lett., 477 (2009), pp. 377–381.
 [9] D. A. Case, C. Scheurer, and R. Brüschweiler, Static and dynamic effects on vicinal scalar J couplings in proteins and peptides: a MD/DFT analysis, J. Am. Chem. Soc., 122 (2000), pp. 10390–10397.
 [10] J. Cavanagh, W. J. Fairbrother, A. G. Palmer III, and N. J. Skelton, Protein NMR spectroscopy: principles and practice, Elsevier Academic Press, 1995.
 [11] G. K.L. Chan and M. HeadGordon, Highly correlated calculations with a polynomial cost algorithm: A study of the density matrix renornalization group, J. Chem. Phys., 116 (2002), p. 4462.
 [12] G. M. Crosswhite, A. C. Doherty, and G. Vidal, Applying matrix product operators to model systems with long–range interactions, Phys. Rev. B, 78 (2008), p. 035116.
 [13] S. V. Dolgov, B. N. Khoromskij, I. V. Oseledets, and D. V. Savostyanov, Computation of extreme eigenvalues in higher dimensions using block tensor train format, Computer Phys. Comm., 185 (2014), pp. 1207–1216.
 [14] S. V. Dolgov, B. N. Khoromskij, and D. V. Savostyanov, Superfast Fourier transform using QTT approximation, J. Fourier Anal. Appl., 18 (2012), pp. 915–953.
 [15] S. V. Dolgov and I. V. Oseledets, Solution of linear systems and matrix inversion in the TTformat, SIAM J. Sci. Comput., 34 (2012), pp. A2718–A2739.
 [16] S. V. Dolgov and D. V. Savostyanov, Alternating minimal energy methods for linear systems in higher dimensions. Part I: SPD systems, arXiv preprint 1301.6068, 2013.
 [17] , Alternating minimal energy methods for linear systems in higher dimensions. Part II: Faster algorithm and application to nonsymmetric systems, arXiv preprint 1304.1222, 2013.
 [18] , Corrected onesite density matrix renormalization group and alternating minimal energy algorithm, in Proc. ENUMATH 2013, accepted, 2014.
 [19] R. S. Dumont, S. Jain, and A. Bain, Simulation of manyspin system dynamics via sparse matrix methodology, J Chem. Phys., 106 (1997), pp. 5928–5936.
 [20] L. E. Edwards, D. V. Savostyanov, A. A. Nevzorov, M. Concistrè, G. Pileo, and I. Kuprov, Gridfree powder averages: On the applications of the FokkerPlanck equation to solid state NMR, J. Magn. Reson., 235 (2013), pp. 121–129.
 [21] L. J. Edwards, D. Savostyanov, Z. Welderufael, D. Lee, and I. Kuprov, Quantum mechanical NMR simulation algorithm for protein–size spin systems, J. Magn. Reson., 243 (2014), pp. 107–113.
 [22] R. R. Ernst, G. Bodenhausen, and A. Wokaun, Principles of nuclear magnetic resonance in one and two dimensions, vol. 14, Clarendon Press Oxford, 1987.
 [23] M. Fannes, B. Nachtergaele, and R. Werner, Finitely correlated states on quantum spin chains, Comm. Math. Phys., 144 (1992), pp. 443–490.
 [24] S. A. Goreinov, I. V. Oseledets, and D. V. Savostyanov, Wedderburn rank reduction and Krylov subspace method for tensor approximation. Part 1: Tucker case, SIAM J. Sci. Comput., 34 (2012), pp. A1–A27.
 [25] W. Hackbusch, Tensor spaces and numerical tensor calculus, Springer–Verlag, Berlin, 2012.
 [26] W. Hackbusch and S. Kühn, A new scheme for the tensor representation, J. Fourier Anal. Appl., 15 (2009), pp. 706–722.
 [27] P. E. Hansen, Carbon–hydrogen spin–spin coupling constants, Progress in NMR Spectroscopy, 14 (1981), pp. 175–295.
 [28] H. J. Hogben, M. Krzystyniak, G. T. P. Charnock, P. J. Hore, and I. Kuprov, Spinach — A software library for simulation of spin dynamics in large spin systems, J Magn. Reson., 208 (2011), pp. 179–194.
 [29] T. K. Huckle, K. Waldherr, and T. SchulteHerbrüggen, Exploiting matrix symmetries and physical symmetries in matrix product states and tensor trains, Linear and Multilinear Algebra, 61 (2013), pp. 91–122.
 [30] E. Jeckelmann, Dynamical density–matrix renormalization–group method, Phys. Rev. B, 66 (2002), p. 045114.
 [31] A. Karabanov, I. Kuprov, G. T. P. Charnock, A. van der Drift, L. J. Edwards, and W. Köckenberger, On the accuracy of the state space restriction approximation for spin dynamics simulations, J Chem. Phys., 135 (2011), p. 084106.
 [32] L. Kay, P. Keifer, and T. Saarinen, Pure absorption gradient enhanced heteronuclear single quantum correlation spectroscopy with improved sensitivity, J. Am. Chem. Soc., 114 (1992), pp. 10663–10665.
 [33] L. E. Kay, M. Ikura, R. Tschudin, and A. Bax, Threedimensional tripleresonance NMR spectroscopy of isotopically enriched proteins, J. Magn. Reson. (1969), 89 (1990), pp. 496–514.
 [34] V. Kazeev, B. Khoromskij, and E. Tyrtyshnikov, Multilevel Toeplitz matrices generated by tensorstructured vectors and convolution with logarithmic complexity, SIAM J. Sci. Comp., 35 (2013), pp. A1511–A1536.
 [35] V. A. Kazeev and B. N. Khoromskij, Lowrank explicit QTT representation of the Laplace operator and its inverse, SIAM J. Matrix Anal. Appl., 33 (2012), pp. 742–758.
 [36] B. N. Khoromskij, Tensorstructured numerical methods in scientific computing: Survey on recent advances, Chemometr. Intell. Lab. Syst., 110 (2012), pp. 1–19.
 [37] B. N. Khoromskij and V. Khoromskaia, Multigrid accelerated tensor approximation of function related multidimensional arrays, SIAM J. Sci. Comput., 31 (2009), pp. 3002–3026.
 [38] A. Klümper, A. Schadschneider, and J. Zittartz, Matrix product ground states for onedimensional spin1 quantum antiferromagnets, Europhys. Lett., 24 (1993), pp. 293–297.
 [39] T. G. Kolda and B. W. Bader, Tensor decompositions and applications, SIAM Review, 51 (2009), pp. 455–500.
 [40] L. B. Krivdin and E. W. Della, Spin–spin coupling constants between carbons separated by more than one bond, Progress in NMR Spectroscopy, 23 (1991), pp. 301–610.
 [41] Ē. Kupce and R. Freeman, Stretched adiabatic pulses for broadband spin inversion, J Magn. Reson. A, 117 (1995), pp. 246–256.
 [42] I. Kuprov, N. WagnerRundell, and P. J. Hore, Polynomially scaling spin dynamics simulation algorithm based on adaptive statespace restriction, J Magn. Reson., 189 (2007), pp. 241–250.
 [43] Y. Kurashige and T. Yanai, Highperfomance ab initio density matrix renormalization group method: Applicability to largescale multireference problems for metal compounds, J. Chem. Phys., 130 (2009), p. 234114.
 [44] S. Liang and H. Pang, Approximate diagonalization using the density matrix renormalizationgroup method: A twodimensionalsystems perspective, Phys. Rev. B, 49 (1994), pp. 9214–9217.
 [45] F. Löhr, J. M. Schmidt, and H. Rüterjans, Simultaneous measurement of and coupling constants in –labeled proteins, J. Am. Chem. Soc., 121 (1999), pp. 11821–11826.
 [46] I. P. Mcculloch, Density matrix renormalization group algorithm and the twodimensional tJ model, Philosophical Magazine B, 81 (2001), pp. 1603–1613.
 [47] I. P. McCulloch and M. Gulácsi, The nonAbelian density matrix renormalization group algorithm, Europhys. Lett., 57 (2002), p. 852.
 [48] V. Murg, F. Verstraete, O. Legeza, and R. M. Noack, Simulating strongly correlated quantum systems with tree tensor networks, Phys. Rev. B, 82 (2010), p. 205105.
 [49] I. Oseledets et al., TTToolbox.
 [50] I. V. Oseledets, Tensortrain decomposition, SIAM J. Sci. Comput., 33 (2011), pp. 2295–2317.
 [51] C. Pérez, F. Löhr, H. Rüterjans, and J. M. Schmidt, Selfconsistent Karplus parametrization of couplings depending on the polypeptide sidechain torsion , J. Am. Chem. Soc., 123 (2001), pp. 7081–7093.
 [52] D. PerezGarcia, F. Verstraete, M. M. Wolf, and J. I. Cirac, Matrix product state representations, Quantum Info. Comput., 7 (2007), pp. 401–430.
 [53] S. Pittel and N. Sandulescu, Density matrix renormalization group and the nuclear shell model, Phys. Rev. C, 73 (2006), p. 014301.
 [54] I. Pižorn, F. Verstraete, and R. M. Konik, Tree tensor networks and entanglement spectra, Phys. Rev. B, 88 (2013), p. 195102.
 [55] S. Ramasesha, S. K. Pati, H. R. Krishnamurthy, Z. Shuai, and J. L. Brédas, Symmetrized densitymatrix renormalizationgroup method for excited states of hubbard models, Phys. Rev. B, 54 (1996), pp. 7598–7601.
 [56] N. F. Ramsey, Magnetic shielding of nuclei in molecules, Phys. Rev., 78 (1950), p. 699.
 [57] , Electron coupled interactions between nuclear spins in molecules, Phys. Rev., 91 (1953), p. 303.
 [58] J. Rotureau, N. Michel, W. Nazarewicz, M. Płoszajczak, and J. Dukelsky, Density matrix renormalization group approach for manybody open quantum systems, Phys. Rev. Lett., 97 (2006), p. 110603.
 [59] D. V. Savostyanov, Fast revealing of mode ranks of tensor in canonical form, Numer. Math. Theor. Meth. Appl., 2 (2009), pp. 439–444.
 [60] J. Schleucher, M. Schwendinger, M. Sattler, P. Schmidt, O. Schedletzky, S. Glaser, O. W. Sørensen, and C. Griesinger, A general enhancement scheme in heteronuclear multidimensional NMR employing pulsed field gradients, J. Biomolecular NMR, 4 (1994), pp. 301–306.
 [61] U. Schollwöck, The density–matrix renormalization group, Rev. Mod. Phys., 77 (2005), pp. 259–315.
 [62] , The densitymatrix renormalization group in the age of matrix product states, Annals of Physics, 326 (2011), pp. 96–192.
 [63] S. Sharma and G. K.L. Chan, Spinadapted density matrix renormalization group algorithms for quantum chemistry, J. Chem. Phys., 136 (2012), p. 124121.
 [64] Y.Y. Shi, L.M. Duan, and G. Vidal, Classical simulation of quantum manybody systems with a tree tensor network, Phys. Rev. A, 74 (2006), p. 022320.

[65]
N. Shibata, Thermodynamics of the anisotropic heisenberg chain calculated
by the density matrix renormalization group method, J. Phys. Soc. Jpn., 66 (1997), pp. 2221–2223.  [66] E. M. Stoudenmire and S. R. White, Studying twodimensional systems with the density matrix renormalization group, Annual Review of Condensed Matter Physics, 3 (2012), pp. 111–128.
 [67] E. L. Ulrich, H. Akutsu, J. F. Doreleijers, Y. Harano, Y. E. Ioannidis, J. Lin, M. Livny, S. Mading, D. Maziuk, Z. Miller, E. Nakatani, C. F. Schulte, D. E. Tolmie, R. K. Wenger, H. Yao, and J. L. Markley, Biomagresbank, Nucleic acids research, 36 (2008), pp. D402–D408.
 [68] F. Verstraete and J. I. Cirac, Matrix product states represent ground states faithfully, Phys. Rev. B, 73 (2006), p. 094423.
 [69] F. Verstraete, J. J. GarcíaRipoll, and J. I. Cirac, Matrix product density operators: Simulation of finite–temperature and dissipative systems, Phys. Rev. Lett., 93 (2004), p. 207204.
 [70] M. Veshtort and R. G. Griffin, SPINEVOLUTION: A powerful tool for the simulation of solid and liquid state NMR experiments, J Magn. Reson., 178 (2006), pp. 248–282.
 [71] B. Vögeli, J. Ying, A. Grishaev, and A. Bax, Limits on variations in protein backbone dynamics from precise measurements of scalar couplings, J. Am. Chem. Soc., 129 (2007), pp. 9377–9385.
 [72] A. C. Wang and A. Bax, Determination of the backbone dihedral angles in human ubiquitin from reparametrized empirical karplus equations, J. Am. Chem. Soc., 118 (1996), pp. 2483–2494.
 [73] X. Wang and T. Xiang, Transfermatrix densitymatrix renormalizationgroup theory for thermodynamics of onedimensional quantum systems, Phys. Rev. B, 56 (1997), pp. 5061–5064.
 [74] S. R. White, Density matrix formulation for quantum renormalization groups, Phys. Rev. Lett., 69 (1992), pp. 2863–2866.
 [75] , Densitymatrix algorithms for quantum renormalization groups, Phys. Rev. B, 48 (1993), pp. 10345–10356.
 [76] , Density matrix renormalization group algorithms with a single center site, Phys. Rev. B, 72 (2005), p. 180403.
 [77] S. R. White and R. L. Martin, Ab initio quantum chemistry using the density matrix renormalization group, J. Chem. Phys., 110 (1999), pp. 4127–4130.
 [78] S. Wouters, P. A. Limacher, D. Van Neck, and P. W. Ayers, Longitudinal static optical properties of hydrogen chains: Finite field extrapolations of matrix product state calculations, J. Chem. Phys., 136 (2012), p. 134110.
 [79] S. Wouters, W. Poelmans, P. W. Ayers, and D. Van Neck, CheMPS2: A free opensource spinadapted implementation of the density matrix renormalization group for ab initio quantum chemistry, Computer Phys. Comm., (2014).
 [80] T. Xiang, J. Lou, and Z. Su, Twodimensional algorithm of the densitymatrix renormalization group, Phys. Rev. B, 64 (2001), p. 104414.
 [81] D. Zgid and M. Nooijen, On the spin and symmetry adaptation of the density matrix renormalization group method, J. Chem. Phys., 128 (2008), p. 014107.