Auxiliary matrix formalism for interaction representation transformations, optimal control and related theories D.L. Goodwin1, L.J. Edwards2, Ilya Kuprov1,* 1
School of Chemistry, University of Southampton, Highfield Campus, Southampton, SO17 1BJ, UK. 2
Wellcome Trust Centre for Neuroimaging, 12 Queen Square, London, WC1N 3BG, UK.
*
Corresponding author (
[email protected] ) 1
Abstract Auxiliary matrix exponential method is used to derive simple and numerically efficient general expressions for the following, historically rather cumbersome and hard to compute, theoretical methods: (1) average Hamiltonian theory following interaction representation transformations; (2) Bloch-Redfield-Wangsness theory of nuclear and electron relaxation; (3) gradient ascent pulse engineering version of quantum optimal control theory; (4) pulsed field gradient treatment in NMR spectroscopy. In the context of spin dynamics, the auxiliary matrix exponential method is more efficient than methods based on matrix factorizations and also exhibits more favourable complexity scaling with the dimension of the Hamiltonian matrix.
2
Introduction In the bestiary of complicated functions that harrow every magnetic resonance theorist, chained exponential integrals involving some square matrices A k and B k occur particularly often: t
tn2
t1
dt1 dt2 ... dtn1 e 0
0
0
A1 t t1
B1e
A 2 t1 t2
B 2 B n 1e Antn1 .
(1)
Examples include perturbative relaxation theories [1-3], reaction yield expressions in radical pair dynamics [4-6], average Hamiltonian theory [7], fidelity functional derivatives in optimal control theory [8,9] and pulsed field gradient propagators in nuclear magnetic resonance [10]. Their common feature is the complexity of evaluation: numerically expensive matrix factorizations are required [3,11,12]. This makes the application of the associated theories difficult when matrix dimension exceeds 103, i.e. for ten spins or more. Consider the example of spin relaxation theory. The currently used techniques for the evaluation of Redfield's integral, which involves a static Hamiltonian H , a rotational correlation function
G and an irreducible spherical component Q of the stochastic Hamiltonian, either require the diagonalization of H followed by the evaluation of a large number of Fourier transforms [1,3]:
iH iD † iH † † iD † G e Q e d G Ve V Q Ve V d 0 nk 0 nk G Ve iD V †Q† V eiD V † d G Vnr e iDrr V †Q† V eiDss Vks* d rs 0 nk rs 0
Vnr V Q V V rs †
rs
†
* ks
G e
i Drr Dss
d ;
(2)
H =VDV †
0
or a matrix-valued numerical quadrature with a large number of time steps [13]. The latter method scales better because matrix sparsity is preserved at every stage and diagonalization is avoided, but the evaluation is still computationally expensive. Such situations are ubiquitous in magnetic resonance and the integrals in Equation (1) are the bottleneck of many practically important simulations (e.g. of NOESY spectroscopy [14]). Ideally, their evaluation should be an elementary function that does not involve expensive matrix operations or numerical quadrature grids. In this communication we propose a solution to this problem, based on the observation that matrix exponentiation, when used judiciously, does not require factorizations and preserves spin operator sparsity [13,15] and on the auxiliary matrix technique [16-18] for the evaluation of the integrals given in Equation (1). The result is a simplification, automation and numerical acceleration of some of the oldest and most useful magnetic resonance theories, particularly for systems involving 3
large numbers of interacting spins. This work is motivated in particular by the practical needs arising during the development of Spinach [15], which is a large-scale spin dynamics simulation library that cannot, as a matter of general policy, rely on expensive matrix factorizations or expect the user to pre-process the Hamiltonian by hand. The methods described below were found, in our practical simulation work, to be user- and developer-friendly as well as numerically efficient. Exponentials of auxiliary matrices A method for computing some of the integrals of the general type shown in Equation (1) was proposed by Van Loan in 1978 [16]. He noted that the integrals in question are solutions to linear block matrix differential equations and suggested that block matrix exponentials are used to compute them. In the simplest case of a single integral: A B e At exp t 0 C 0
At1 Ct1 e B e dt 1 0 Ct e t
e
At
(3)
Van Loan's method was subsequently refined by Carbonell et al., who derived a convenient expression using the exponential of a block-bidiagonal auxiliary matrix [18]: A11 0 M 0 0 0
A12
0
A 22
A 23
0 0
A 33 0
0
0
0 0 0 , A k 1,k 0 A kk 0
0
B11 B12 0 B 22 exp Mt 0 0 0 0 0 0
B1k B2k B k 1,k 0 B kk
B13 B 23 B 33 0 0
(4)
The integrals in question populate the rows of the resulting block matrix, e.g.: t
t1
tk 2
0
0
0
B1k dt1 dt2 ... dtk 1 e A11 t t1 A12 e A22 t1 t2 A 23 A k 1,k e Akk tk 1
(5)
Equations (4) and (5) are particularly appealing in the context of magnetic resonance spectroscopy – spin Hamiltonians are guaranteed to be sparse in the Pauli basis [19-21] and their exponential propagators are also sparse when iHt 1 , if due care is taken to eliminate insignificant elements after each multiplication [15] in the "scaled and squared" Taylor series procedure: e
Mt
k 0
Mt k!
k
;
e
Mt
eMt
2n
2
2 ...
2
(6)
Out of the multitude of "dubious ways" [22] of computing matrix exponentials, the Taylor series with scaling and squaring is preferred in practice because it is compatible with dissipative dynamics (Chebyshev method is not [23]), only involves matrix multiplications (Padé method requires a
4
costly and perilous inverse [16]), uses minimal memory resources (Newton polynomials are more expensive [24]), does not require precise scaling (Newton and Chebyshev methods do [23,24]) and is resistant, in terms of stability and numerical accuracy [25], to inexperienced chemists with noisy and poorly conditioned matrices. For a publicly available simulation package that must preserve sparsity (matrix dimensions are often in the millions) and run reliably in a large variety of settings [15,20], the latter consideration is actually decisive – in our experience, the last few decimal places in the double precision arithmetic are rarely important, but the real world ruggedness always is. We demonstrate below that Equations (4)-(6) remove the problems associated with the calculation of nested exponential integrals from several widely used magnetic resonance simulation methods. The rest of this paper goes through those we could identify, reformulates them to use the auxiliary matrix formalism and benchmarks the results. Rotating frame transformations A common first stage of spin dynamics simulations is the interaction representation transformation wherein the Hamiltonian matrix is split into the "large and simple" part H 0 and the "small and complicated" part H 1 , and the following unitary transformation is applied [11,12]:
σ t eiH0t ρ t e iH0t
H1R t eiH0t H1eiH0t
d ρ t i H 0 H1 , ρ t dt
d σ t i H1R t , σ t dt
(7)
An important next step is to make H1R t effectively time-independent by averaging its exponential propagator over the period of exp iH 0t . This is achieved by combining the solution of the laboratory frame Liouville - von Neumann equation ρ t e
i H 0 H1 t
ρ 0 e
i H 0 H1 t
(8)
with the definition of the interaction representation (aka "rotating frame") σ t eiH0t ρ t e iH0t
(9)
and observing that σ 0 ρ 0 : σ t e iH 0 t e
i H 0 H1 t
σ 0 e
i H 0 H1 t iH 0 t
e
U t σ 0 U† t
(10)
The effective propagator in the rotating frame is therefore U t e iH 0 t e
i H 0 H1 t
(11)
and the effective rotating frame Hamiltonian over a time interval 0,T is the logarithm of U T
H
i iH0T i H0 H1 T ln e e T 5
(12)
If T is now chosen to be the period of exp iH 0t , the propagator exp iH 0T becomes equal to the unit matrix and Equation (12) acquires a form similar to that seen in the theory of generalized cumulant expansions [2,26]:
H
i i H0 H1 T ln e T
(13)
in which the reader should resist the temptation to cancel the exponential and the logarithm – matrix logarithms are defined in the sense of the principal value which need not be equal to the matrix that has gone into the exponential [27]. With the above mentioned condition on T , Equation (13) is a neat formulation of the exact rotating frame average Hamiltonian and a useful generating function for the perturbative expansion that we are about to derive. Under the typical magnetic resonance rotating frame assumptions H1 H 0 and therefore the
H1 term under the exponential in Equation (13) is a small correction to the H 0 term. The corresponding Taylor series with respect to a small parameter is: H
i i H 0 H1 T i ln e T T
n i H H T ln e 0 1 n n 1
n 0
(14)
n!
The logarithm disappears after the first differentiation, leaving remarkably simple expressions for the perturbative corrections to the effective Hamiltonian in which nested commutators do not occur, the number of matrix terms is linear with respect to the approximation order and further terms may be obtained easily by repeated application of the product rule: H 0 i 1 H D†0 D1 T i 2 H D1† D1 D†0 D2 2T i 3 H D†2 D1 2D1† D2 D†0 D3 6T 0
Dk
k e
i H 0 H1 T
k
(15) 0
The derivatives D k of exp i H 0 H1 T with respect to at 0 are known [17] to have the form that matches the auxiliary matrix integrals in Equation (5): T
Dk T keTA e tA BDk 1 t dt , 0
D0 t e , tA
A iH 0 ,
(16)
B iH1
This makes them easy to compute – the derivatives of exp A B t with respect to at 0 may be obtained by constructing and exponentiating a very sparse block-bidiagonal matrix M prescribed by Equation (4): 6
D0 0 exp tM 0 0 0
A B 0 0 0 0 A B 0 0 M0 0 A 0 0 0 0 B 0 0 0 0 A
D1 1! D2 2! Dk k ! D0 D1 1! Dk 1 k 1 ! 0 D0 0 0 D1 1! 0 0 0 D0
(17)
and extracting the first block row from the result. Because T is a period of exp iH 0t , and H1 H 0 the 2-norm of tM is approximately 2 , meaning that the exponential of tM is also sparse [21]. For the same reason, the D0 terms in Equations (15)-(17) are actually unit matrices. The primary advantage of this method over the standard average Hamiltonian theory [7] is the simplicity of implementation and the possibility of automation – perturbative corrections to the rotating frame transformation used to be an arduous manual process that had to be endured afresh for each new class of spin systems. The procedure described above is also an improvement in the sense that its numerical implementations require no knowledge of the content of H 0 and H1 – it may be incorporated into highly general simulation codes without the logistical overhead of storing every significant magnetic resonance spectroscopy assumption and rotating frame layout. Spin relaxation theories Exponential integrals of the form discussed above also appear in the context of Bloch-RedfieldWangsness spin relaxation theory [1,3], in which the relaxation superoperator R arising from the stochastic rotational modulation of spin interaction anisotropies has the following general form:
R Q km Gkmpq e iH0 Q †pq eiH0 d kmpq
(18)
0
where H 0 is the time-independent part of the spin Hamiltonian commutation superoperator, Q km are the 25 irreducible spherical components of its anisotropic part and the rotational correlation functions Gkmpq t are defined as ensemble averages (denoted by angular brackets) of products of second-rank Wigner D functions of molecular orientation [13]: Gkmpq t Dkm t Dpq 0 2
2*
(19)
In systems undergoing rotational diffusion these functions may be expressed as linear combinations of decaying exponentials [28,29]: G t an e nt
(20)
n
that are scalars and therefore commute with all matrices. With this observation in place, we can conclude that individual matrix-valued integrals in Equation (18) have the following general form: 7
e
iH 0 t
Qei ( H 0 i 1) t dt
(21)
0
These integrals used to be evaluated by diagonalizing H 0 and computing the Fourier transform analytically [1,3,11,12] as shown in Equation (2). Because eigenvector arrays of sparse matrices are usually dense, this has limited the applicability of BRW theory to systems with fewer than about five spins. A useful workaround was found by Kuprov in 2011 – he has pointed out that exponentiating H 0 is cheaper than diagonalizing it and suggested a numerical route to the same integral [13] that extended the range of matrix dimensions into hundreds of thousands [14]. In this communication we point out that Equation (21) is a case of an auxiliary exponential relationship: T
e
iH 0 t
Qei ( H0 i 1)t dt A †B
0
iH 0 A B exp 0 0 C
T iH 0 1 Q
(22)
in which we made use of the fact that H 0 is Hermitian. The upper integration limit should be set according to the accuracy goal for the relaxation superoperator as suggested by Kuprov [13]:
e
T
t
dt e t dt
T ln 1
(23)
0
where is the desired relative accuracy of the integral. Equation (22) is easier to implement than a numerical quadrature and cheaper to compute than the diagonalization pathway. As of version 1.5, this procedure is the default in Spinach [15] – matrix dimensions that are routinely processed by the current version of Spinach on desktop computers are in the millions. Another aspect of semi-classical spin relaxation theories that may be cast into the auxiliary matrix form is the thermal equilibrium problem. Because spatial degrees of freedom are treated classically in the BRW theory and the back-action by the spins on those degrees of freedom is ignored [1,3,12], the resulting relaxation superoperator drives the system state vector t to the infinite temperature state rather than the thermal equilibrium: d t iH t t R t dt
lim t 0 t
(24)
where H t is the Hamiltonian commutation superoperator (a Hermitian matrix) and R is the relaxation superoperator (a symmetric negative-definite matrix). A simple empirical workaround for this problem is to switch the relaxation target to the thermal equilibrium state [1,3,12]: d t iH t t R t eq dt
8
(25)
where eq is the thermal equilibrium state vector. This transformation has the unfortunate effect
of making the equation of motion inhomogeneous, but the problem may be remedied by using the following auxiliary matrix transformation: d t i H t R R t dt eq 0 0 eq
(26)
which returns the equation of motion into the homogeneous form at the cost of doubling the matrix dimension. Equation (26) appears to be a known trick in magnetic resonance for which we claim no credit – in fact, more computationally efficient thermalization schemes exist [30,31] – but it should be mentioned in an auxiliary matrix paper for the sake of completeness. Optimal control theories
The task of taking a spin system from one state to another to a specified accuracy with minimal expenditure of time and energy, with the emphasis on the word minimal, is increasingly important in magnetic resonance [32-35]. Optimal solutions are usually found numerically, by maximizing "fidelity" – the overlap between the final state of the system and the desired destination state [9]:
T ˆˆ ˆ f Re ˆ ˆ T Re ˆ exp O i H t iRˆ dt ˆ 0 (27) 0 ˆ where ˆ 0 is the initial condition, ˆ is the desired destination state, Hˆ t is the Hamiltonian ˆ commutation superoperator, Rˆ is the relaxation superoperator and exp O denotes a time-ordered ˆ ˆ ˆ exponential. The Liouvillian superoperator of the spin system, defined as Lˆ t Hˆ t iRˆ , typiˆ cally contains the "drift" part Lˆ0 that cannot be influenced and the controllable part that the magnetic resonance instrument can vary within certain limits: ˆ ˆ ˆ k Lˆ t Lˆ0 c t Lˆk
(28)
k
ˆ where Lˆk are the control operators (usually radiofrequency and microwave Zeeman operators) and
t are their time-dependent coefficients. The maxima of the fidelity functional f with respect k to the control sequences c t , subject to experimental constraints on feasible amplitudes and c
k
frequencies, are the subject of the mathematical branch of the optimal control theory [36,37]. The standard modus operandi in magnetic resonance simulations is to discretize time and solve Equation (27) using thin slice propagators [8,33]. On a finite grid of time points 0 tn T , the control sequences c
k
t
become vectors:
c
k
t cn k , 9
tn 1 t tn
(29)
and the task of minimizing the fidelity functional becomes a high-dimensional non-linear optimization problem. Numerical optimization is remarkably well researched [38] and, in this context, gradient descent, quasi-Newton and Newton-Raphson families of methods are generally appropriate. From the computational efficiency point of view, the central problem is therefore the calculation of first and second derivatives of the fidelity functional:
J j cn cm 2
k
ˆ f ˆˆ ˆˆ Pˆn ˆˆ ˆ Re ˆ PN Pn 1 k Pn 1 Pˆ1 ˆ 0 k cn cn ˆ ˆ ˆˆ ˆˆ Pˆn ˆˆ ˆ Pˆm ˆˆ ˆ ˆ Re | PN Pn 1 k Pn 1 Pm 1 j Pm 1 Pˆ1 0 cn cm
(30)
where time slice propagators (assuming a fixed time grid step t ) are defined as:
ˆ ˆ k ˆ Pˆn exp i Lˆ0 cn Lˆk t k
(31)
The primary task therefore is to calculate first and second directional derivatives of this matrix exponential. This problem is very well researched [17] and the solution has already been mentioned in Equation (17). Restricting our attention to first and second derivatives specifically, we obtain the following auxiliary matrix relations:
Pn 0 0
1 2 Pn k k k j Lk 0 cn 2 cn cn L 0 cn L k k Pn k t (32) Pn 0 L 0 cn L k Lj exp i j cn k 0 Pn 0 0 L 0 cn k L k k ˆ Because each individual propagator Pˆn in Equation (31) only depends on a small number of coefPn
ficients cn k (the number of control channels enumerated by k is usually two in magnetic resonance experiments), the number of auxiliary matrix exponentials that need to be computed to obtain the fidelity functional Hessian is linear with respect to the number of time steps in the control sequence – by the usual standards of the Newton-Raphson method, the Hessian is extraordinarily cheap for the same reason that makes the gradient cheap [9,33]. Not using second order information in such a situation would be rather wasteful and we therefore recommend that this is always done – an exploration of the convergence properties of the resulting optimization algorithm is outside the scope of this work and will be published separately. Pulsed field gradients 10
The application of a B0 field gradient to the NMR sample makes the Hamiltonian of the spin system depend on spatial location of the molecule in the sample: k Hˆ r Hˆ 0 g T r Lˆ , Lˆ k LˆZ ,
(33)
k
where g is the B0 field gradient vector, r is the position vector, k are magnetogyric ratios and the index k runs over all spins in the sample. The propagator associated with this Hamiltonian is
Pˆ r , t exp i Hˆ 0 g T r Lˆ t .
(34)
Because the radiofrequency coils of the NMR instrument detect the total sample magnetization, the spatial coordinates must be integrated over in the expression for the detected signal:
ˆ ˆ f t w r ˆ exp i Hˆ 0 g T r Lˆ t ˆ 0 d 3r ˆ ˆ ˆ w r exp i Hˆ 0 g T r Lˆ t d 3r ˆ 0
(35)
The task is therefore reduced again to computing an integral involving a matrix exponential with a weight w r that corresponds to the sample probability density. A convenient feature of the high-field approximation in NMR is that the Hamiltonian Hˆ 0 commutes with Lˆ , therefore:
w r exp i
Hˆˆ g r Lˆˆ t d r T
3
0
(36) T ˆˆ 3 ˆˆ exp iH 0t w r exp i g r L t d r The probability density w r in a box of size aX , aY , aZ may be expanded in plane waves: w(r ) wpqs e ipx aX e iqy aY e isz aZ (37) pqs
This leads to the possibility of separating the integral in Equation (36):
T r Lˆˆ t d 3 r w Uˆ X Uˆ Y Uˆ Z pqs p q s pqs
w r exp i g
(38)
with the elementary integrals now compliant with the general form that allows the auxiliary matrix shortcut in their calculation:
p a g Lˆˆ t x dx ˆ exp i q a g Lˆ t y dy ˆ exp i s a g Lˆ t z dz
Uˆ pX
aX
exp i
X
X
Y
Y
0
Uˆ q Y
aY
0
Z Uˆ s
aZ
Z
0
11
Z
(39)
The following auxiliary matrix relation returns the integrals in question: 0 1 exp 0 i p aX g X Lˆˆ t
X Uˆ p 1 a X 0 exp i p aX g X Lˆˆ t aX
(40)
Y and similarly for Uˆ q and Uˆ r Z . In the most common (in NMR spectroscopy) case of a homospoil
gradient and a uniform concentration distribution, this was recently shown to accelerate the corresponding simulations considerably [10].
Radical pair dynamics Haberkorn [6] and Jones-Hore [5] models of radical pair recombination stipulate the following equations of motion for the spin density operator ˆ t : d ˆ t k k ˆ i Hˆ , ˆ t S ˆ t PˆS PˆS ˆ t T ˆ t PˆT PˆT ˆ t Rˆ dt 2 2 ˆ d t ˆ i Hˆ , ˆ t kS kT ˆ t kS PˆT ˆ t PˆT kT PˆS ˆ t PˆS Rˆ dt
(41)
ˆ where Hˆ is the spin Hamiltonian, Rˆ is the spin relaxation superoperator, PˆS,T are two-electron singlet and triplet projection operators, kS,T are singlet and triplet radical pair recombination rates and the assumption of first-order chemical kinetics unencumbered by spatial diffusion effects is made for the radical pair. Convenient solutions are possible in the adjoint representation (aka Liouville space [11]), where both models acquire the same general form d ˆ ˆ t iLˆ ˆ t dt
(42)
in which ˆ t denotes a column-wise vectorization of ˆ t and
i ˆ ˆ ˆ ˆ ˆ LˆH Hˆ kS PˆS kT PˆT iRˆ 2 ˆ ˆ ˆ ˆ LˆJH Hˆ i kS kT 1ˆ ikS PˆT PˆTT ikT PˆS PˆST iRˆ
(43)
for Haberkorn and Jones-Hore model respectively. In Equation (43), the upper T index indicates matrix transpose operation. Commutation and anti-commutation superoperators are defined as
ˆ Oˆ Oˆ 1ˆ 1ˆ Oˆ T
(44)
The target properties in radical pair simulations are singlet and triplet product yields YS,T t , defined as time integrals over the singlet and triplet populations: t
YS t kS
t
ˆ
ˆ
ˆ YT t kT PˆT e iLt ˆ 0 dt
ˆ PˆS e iLt ˆ 0 dt
0
0
12
(45)
where ˆ 0 is the density matrix for the initial state. Evaluation of these integrals is the slowest stage of radical pair dynamics simulations because either a diagonalization [4] or a matrix inverse-timesvector operation [39] is normally involved. It is clear, however, that both may be avoided because the yields in Equation (45) are matrix elements of an exponential integral that may be computed using a special case of Equation (3) with 0 denoting an all-zero vector of the same size as ˆ : t
YS t kS
t
ˆ
ˆ PˆS e iLt ˆ 0 dt kS PˆS
ˆˆ
iLt e dt ˆ 0 kS
0
0
Pˆ
S
0 1 0 0 exp t 0 iLˆˆ ˆ 0
(46)
and similarly for the triplet yield. The importance of using Krylov propagation for exp tA v type operations [40] should be stressed – it is not necessary to compute the entire matrix exponential when only one element is required by Equation (46). Similar solutions are possible for the simplified kinetic models in which radical pair recombination is assumed to be a first-order process with a time-dependent recombination probability [4]: t
t
ˆ ˆ YS t Tr e iHt ˆ 0 eiHt PˆS f t dt
ˆ ˆ YT t Tr e iHt ˆ 0 eiHt PˆT f t dt
0
(47)
0
The recombination probability function f t may be expanded as f t k n e kn t
(48)
n
and the integrals in Equation (47) transformed into a form that permits the application of the auxiliary matrix exponential technique. For the singlet yield case: t e iHtˆ ˆ 0 eiHtˆ PˆS ke kt dt kTr e iHtˆ ˆ 0 eiHtˆ PˆSe kt dt Tr 0 0 t
t ˆ ˆ ˆ kTr PˆS e iHt ˆ 0 e(iH kE )t dt 0
(49)
where Eˆ is an identity matrix of the same dimension as the Hermitian operator Hˆ . The integral under the trace is another instance of the problem already treated in Equation (22): t
e
ˆ iHt
ˆ
ˆ
ˆ 0ei ( H ikE )t dt A †B
0
iHˆ A B exp 0 0 C
t iHˆ kEˆ
ˆ 0
(50)
The advantage of this method over the diagonalization technique is that matrix sparsity is preserved in the Hamiltonian exponentiation procedure [15], leading to large memory and CPU time
13
savings [4]. Another practically relevant reminder is that Frobenius scalar product may be computed very efficiently, particularly for sparse matrices, by element-wise multiplication: A B =Tr A †B Tot A* B
(51)
where A † indicates the conjugate-transpose of A , A* indicates the element-wise complex conjugate of A , denotes element-wise matrix product and Tot stands for the total sum. Conclusions The primary advantage of auxiliary matrix techniques over other methods of computing the integrals involving matrix exponentials is quite specific to spin dynamics, where Hamiltonians and their small step exponentials are guaranteed to be very sparse [21]. Their various factorizations, however, are usually dense and must therefore be avoided whenever possible. Recognising this fact and recasting the relevant mathematical tools in a way that specifically prefers matrix exponentiation to factorization leads to significant improvements in numerical efficiency and scaling. Acknowledgements We are grateful to Sophie Schirmer for useful discussions. This work was made possible by EPSRC through a grant to IK group (EP/H003789/1) and a CDT studentship to DG. The European Commission has made a negative contribution to the progress of this research by demanding a surreal amount of paperwork and refusing to pay for investigator time on grant 297861/QUAINT. References [1] [2] [3] [4]
[5] [6] [7] [8]
A.G. Redfield, On the theory of relaxation processes, IBM Journal of Research and Development 1 (1957) 19-31. R. Kubo, Generalized Cumulant Expansion Method, Journal of the Physical Society of Japan 17 (1962) 1100. M. Goldman, Formal theory of spin-lattice relaxation, Journal of Magnetic Resonance 149 (2001) 160-187. C.R. Timmel, U. Till, B. Brocklehurst, K.A. McLauchlan, P.J. Hore, Effects of weak magnetic fields on free radical recombination reactions, Molecular Physics 95 (1998) 7189. J.A. Jones, P.J. Hore, Spin-selective reactions of radical pairs act as quantum measurements, Chemical Physics Letters 488 (2010) 90-93. R. Haberkorn, Density matrix description of spin-selective radical pair reactions, Molecular Physics 32 (1976) 1491-1493. U. Haeberlen, J.S. Waugh, Coherent Averaging Effects in Magnetic Resonance, Physical Review 175 (1968) 453-467. Z. Tošner, T. Vosegaard, C. Kehlet, N. Khaneja, S.J. Glaser, N.C. Nielsen, Optimal control in NMR spectroscopy: Numerical implementation in SIMPSON, Journal of Magnetic Resonance 197 (2009) 120-134. 14
[9] [10]
[11] [12] [13] [14]
[15]
[16] [17] [18] [19] [20] [21] [22] [23]
[24]
[25] [26]
[27] [28]
P.d. Fouquieres, S.G. Schirmer, S.J. Glaser, I. Kuprov, Second order gradient ascent pulse engineering, Journal of Magnetic Resonance 212 (2011) 412 - 417. L.J. Edwards, Simulation of coherence selection by pulsed field gradients in liquid-state NMR using an auxiliary matrix formalism, Journal of Magnetic Resonance 240 (2014) 95101. R.R. Ernst, G. Bodenhausen, A. Wokaun, Principles of nuclear magnetic resonance in one and two dimensions, Oxford University Press, 1987. A. Abragam, The principles of nuclear magnetism, Clarendon Press, 1961. I. Kuprov, Diagonalization-free implementation of spin relaxation theory for large spin systems, Journal of Magnetic Resonance 209 (2011) 31-38. L.J. Edwards, D.V. Savostyanov, Z.T. Welderufael, D.H. Lee, I. Kuprov, Quantum mechanical NMR simulation algorithm for protein-size spin systems, Journal of Magnetic Resonance 243 (2014) 107-113. H.J. Hogben, M. Krzystyniak, G.T. Charnock, P.J. Hore, I. Kuprov, Spinach - a software library for simulation of spin dynamics in large spin systems, Journal of Magnetic Resonance 208 (2011) 179-194. C.F. Van Loan, Computing integrals involving the matrix exponential, IEEE Transactions on Automatic Control 23 (1978) 395-404. I. Najfeld, T.F. Havel, Derivatives of the Matrix Exponential and Their Computation, Advances in Applied Mathematics 16 (1995) 321-375. F. Carbonell, J.C. Jimenez, L.M. Pedroso, Computing multiple integrals involving matrix exponentials, Journal of Computational and Applied Mathematics 213 (2008) 300-305. R.S. Dumont, S. Jain, A. Bain, Simulation of many-spin system dynamics via sparse matrix methodology, Journal of Chemical Physics 106 (1997) 5928-5936. M. Veshtort, R.G. Griffin, SPINEVOLUTION: A powerful tool for the simulation of solid and liquid state NMR experiments, Journal of Magnetic Resonance 178 (2006) 248-282. L.J. Edwards, I. Kuprov, Parallel density matrix propagation in spin dynamics simulations, Journal of Chemical Physics 136 (2012). C. Moler, C.F. Van Loan, Nineteen Dubious Ways to Compute the Exponential of a Matrix, SIAM Review 20 (1978) 801-836. M. Ndong, H. Tal-Ezer, R. Kosloff, C.P. Koch, A Chebychev propagator with iterative time ordering for explicitly time-dependent Hamiltonians, Journal of Chemical Physics 132 (2010). W. Huisinga, L. Pesce, R. Kosloff, P. Saalfrank, Faber and Newton polynomial integrators for open-system density matrix propagation, Journal of Chemical Physics 110 (1999) 55385547. C.F. Van Loan, The Sensitivity of the Matrix Exponential, SIAM Journal on Numerical Analysis 14 (1977) 971-981. B. Yoon, J.M. Deutch, J.H. Freed, Comparison of Generalized Cumulant and Projection Operator Methods in Spin-Relaxation Theory, Journal of Chemical Physics 62 (1975) 46874696. A.H. Al-Mohy, N.J. Higham, Improved Inverse Scaling and Squaring Algorithms for the Matrix Logarithm, Siam Journal on Scientific Computing 34 (2012) C153-C169. G. Lipari, A. Szabo, Model-Free Approach to the Interpretation of Nuclear MagneticResonance Relaxation in Macromolecules: 1. Theory and Range of Validity, Journal of the American Chemical Society 104 (1982) 4546-4559.
15
[29] G. Lipari, A. Szabo, Model-Free Approach to the Interpretation of Nuclear MagneticResonance Relaxation in Macromolecules: 2. Analysis of Experimental Results, Journal of the American Chemical Society 104 (1982) 4559-4570. [30] M.H. Levitt, L. Di Bari, The homogeneous master equation and the manipulation of relaxation networks, Bulletin of Magnetic Resonance 16 (1994) 94-114. [31] T. Levante, R. Ernst, Homogeneous versus inhomogeneous quantum-mechanical master equations, Chemical Physics Letters 241 (1995) 73-78. [32] K. Kobzar, T.E. Skinner, N. Khaneja, S.J. Glaser, B. Luy, Exploring the limits of broadband excitation and inversion pulses, Journal of Magnetic Resonance 170 (2004) 236-243. [33] N. Khaneja, T. Reiss, C. Kehlet, T. Schulte-Herbruggen, S.J. Glaser, Optimal control of coupled spin dynamics: design of NMR pulse sequences by gradient ascent algorithms, Journal of Magnetic Resonance 172 (2005) 296-305. [34] Z. Tosner, T. Vosegaard, C. Kehlet, N. Khaneja, S.J. Glaser, N.C. Nielsen, Optimal control in NMR spectroscopy: Numerical implementation in SIMPSON, Journal of Magnetic Resonance 197 (2009) 120-134. [35] M. Nimbalkar, R. Zeier, J.L. Neves, S.B. Elavarasi, H. Yuan, N. Khaneja, K. Dorai, S.J. Glaser, Multiple-spin coherence transfer in linear Ising spin chains and beyond: Numerically optimized pulses and experiments, Physical Review A 85 (2012) 012325. [36] L.S. Pontriagin, The mathematical theory of optimal processes, Pergamon, 1964. [37] F.L. Lewis, D.L. Vrabie, V.L. Syrmos, Optimal control, Wiley, 2012. [38] J. Nocedal, S.J. Wright, Numerical optimization, Springer, 2006. [39] H.J. Hogben, P.J. Hore, I. Kuprov, Strategies for state space restriction in densely coupled spin systems with applications to spin chemistry, Journal of Chemical Physics 132 (2010) 174101. [40] R.B. Sidje, Expokit: A software package for computing matrix exponentials, Acm Transactions on Mathematical Software 24 (1998) 130-156.
16
