Efficient geometric integrators for nonadiabatic quantum dynamics. I. The adiabatic representation

Geometric integrators of the Schr\"{o}dinger equation conserve exactly many invariants of the exact solution. Among these integrators, the split-operator algorithm is explicit and easy to implement, but, unfortunately, is restricted to systems whose Hamiltonian is separable into a kinetic and potential terms. Here, we describe several implicit geometric integrators applicable to both separable and non-separable Hamiltonians, and, in particular, to the nonadiabatic molecular Hamiltonian in the adiabatic representation. These integrators combine the dynamic Fourier method with recursive symmetric composition of the trapezoidal rule or implicit midpoint method, which results in an arbitrary order of accuracy in the time step. Moreover, these integrators are exactly unitary, symplectic, symmetric, time-reversible, and stable, and, in contrast to the split-operator algorithm, conserve energy exactly, regardless of the accuracy of the solution. The order of convergence and conservation of geometric properties are proven analytically and demonstrated numerically on a two-surface NaI model in the adiabatic representation. Although each step of the higher order integrators is more costly, these algorithms become the most efficient ones if higher accuracy is desired; a thousand-fold speedup compared to the second-order trapezoidal rule (the Crank-Nicolson method) was observed for wavefunction convergence error of $10^{-10}$. In a companion paper [J. Roulet, S. Choi, and J. Van\'{\i}\v{c}ek (2019)], we discuss analogous, arbitrary-order compositions of the split-operator algorithm and apply both types of geometric integrators to a higher-dimensional system in the diabatic representation.


I. INTRODUCTION
Separating electronic from nuclear degrees of freedom leads to the Born-Oppenheimer approximation 1,2 and the intuitive picture of electronic potential energy surfaces.5][6][7][8] To address such processes, one can forget the Born-Oppenheimer picture and treat electrons and nuclei on the same footing, 9,10 use an exact factorization 11,12 of the molecular wavefunction, or, most commonly, determine which Born-Oppenheimer states are significantly coupled 13,14 and then solve the time-dependent Schrödinger equation with a molecular Hamiltonian that contains the nonadiabatic couplings.Below, we will only consider the third, yet the most traditional way to treat quantum nonadiabatic dynamics.
An approach particularly suited to study the nonadiabatic population dynamics of large chemical systems is the ab initio multiple spawning 15,16 and related methods, all of which represent the wavefunction by a superposition of time-dependent Gaussian basis functions moving along classical 17,18 or variational 19,20 trajectories.If high accuracy is required and especially if the Hamiltonian can be expressed as a sum of products of one-dimensional operators, a nonadiabatic algorithm of choice is the multiconfigurational time-dependent Hartree (MCTDH) method 21,22 or its multilayer extension, 23 which expand the state using orthogonal time-dependent basis functions.The power of the MCTDH method relies on the fact that only a small fraction of the tensor-product Hilbert space is typically accessible during the time of interest; sparse-grid methods 24,25 also take advantage of this phenomenon.
However, there are systems, in which the full Hilbert space is accessible, and then full grid or time-independent basis sets are preferable. 25,26ere also exist situations, where, in addition to prescribed accuracy, it pays off to conserve certain invariants of the exact solution exactly, regardless of the accuracy of the wavefunction.Because the above-mentioned methods typically conserve none or only some of these invariants, other methods, called geometric integrators, 27 are needed in this setting.
The geometric integrators acknowledge that the Schrödinger equation is special, and not just another general differential equation.Using these integrators can be likened to realizing that the Earth is not flat but round, and even approximate models of its surface should take this curvature into account.Geometric integrators are highly exploited in classical molecular dynamics, where the deceptively simple Verlet algorithm, 28,29 despite its only second-order accuracy, results in exact conservation of D invariants in a D-dimensional system, where D can easily reach thousands or millions in state-of-the art simulations of proteins.1][32] Symmetric compositions of the first-order split-operator algorithms, 25,32 including the standard second-order splitting, 31 are unitary, symplectic, stable, symmetric, and time-reversible regardless of the size of the time step.
4][35][36][37] In a companion paper 37 (below referred to as Paper II), we implement such higher-order compositions for the nonadiabatic quantum molecular dynamics in the diabatic representation.
Although the split-operator algorithms preserve numerous geometric properties of interest of the exact evolution operator, their use is limited to systems with Hamiltonians separable into a sum Ĥ = T (p) + V (q) of two terms, the first depending only on the momentum operator and the second only on the position operator.One must use a different time propagation scheme for systems with a non-separable Hamiltonian; for example, the nonadiabatic dynamics in the adiabatic representation or particles in crossed electric and magnetic fields.
The explicit Euler method is the simplest integrator applicable to non-separable Hamiltonians; it is, however, unstable. 27,38The implicit Euler method is stable regardless of the size of the time step but requires solving a large, although sparse, system of linear equations at every time step; furthermore, the method fails to preserve the unitarity, time reversibility, energy conservation, and other geometric properties of the exact evolution operator.The second-order differencing method [39][40][41] introduces symmetry by combining the forward and backward step of the explicit Euler method.It is explicit and stable for small enough time steps, but does not conserve the norm or energy exactly.
Another issue with the second-order differencing is that a much too small time step is required to obtain an accurate solution. 42This problem has been addressed by using the Chebyshev 43 and short iterative Lanczos algorithms; 41,44,45 both methods increase remarkably the efficiency of numerical integration by effectively approximating the exact evolution operator.However, these two methods are neither time-reversible nor symplectic, and the Chebyshev propagation scheme does not even conserve the norm.
To address either the low accuracy or nonconservation of geometric properties by various nonadiabatic integrators, we propose time propagation schemes based on symmetric compositions of the trapezoidal rule (also known as the Crank-Nicolson method 30,46 ) or implicit midpoint method.As we show below, because these elementary methods are unitary, symplectic, energy conserving, stable, symmetric, and time-reversible, so are their symmetric compositions.4][35] Methods with higher orders of accuracy are useful for obtaining highly accurate solutions because, for that purpose, they are more efficient than the second-order algorithms.Although each time step of a higher-order method costs more, the solution with the same accuracy can be obtained using a larger time step and, hence, a smaller total number of steps in comparison to lower-order methods.
The final benefit of the proposed geometric integrators is the simple, abstract, and general implementation of the compositions of the trapezoidal rule and implicit midpoint methods; indeed, even these "elementary" methods are, themselves, compositions of simpler explicit and implicit Euler methods.
In the adiabatic representation, the proposed integrators cannot be fully compared with the integrators based on the compositions of the split-operator algorithm, which are only applicable to separable Hamiltonians.Both types of integrators, however, can be used in the diabatic representation, which is the focus of Paper II. 37 We, therefore, compare the two integrators there, using a one-dimensional model 47 of NaI and a three-dimensional model 48 of pyrazine.
The remainder of this paper is organized as follows: In Section II, after defining geometric properties of the exact evolution operator, we discuss their breakdown in elementary methods and recovery in the proposed symmetric compositions of the trapezoidal rule and implicit midpoint methods.Next, we present the dynamic Fourier method for its ease of implementation and the exponential convergence with the number of grid points.Yet, the proposed integrators can be combined with any other basis or grid representation.We conclude Section III by discussing the relationship between the molecular Hamiltonians in the adiabatic and diabatic representations.In Section III, the convergence properties and conservation of geometric invariants by various methods are analyzed numerically on a two-surface NaI model 47 in the adiabatic representation.This system has a non-separable Hamiltonian due to an avoided crossing between its potential energy surfaces and a corresponding region of large nonadiabatic momentum coupling.Section IV concludes the paper.

II. THEORY
For a time-independent Hamiltonian Ĥ, the time-dependent Schrödinger equation has the formal solution ψ(t) = Û (t)ψ(0), where ψ(0) is the initial state and Û (t) the so-called evolution operator.The exact evolution operator is linear (in particular, independent of the initial state), reversible, stable, and, moreover, conserves both the norm and energy of the quantum state.Let us define and discuss these and other geometric properties of the exact evolution operator because they are also desirable in approximate numerical evolution operators Ûappr (t).

A. Geometric properties of the exact evolution operator
An operator Û on a Hilbert space is said to preserve the norm ψ := ψ|ψ 1/2 if Û ψ = ψ .For linear operators Û , preserving the norm is equivalent to preserving the inner product, where Û † is the Hermitian adjoint of Û .The preservation of inner product is, therefore, equivalent to the condition that Û † Û be the identity operator, or that Such an operator Û is said to be unitary.The exact evolution operator is unitary since An operator Û is said to be symplectic if it preserves the symplectic two-form ω(ψ, φ), i.e., a nondegenerate skew-symmetric bilinear form on the Hilbert space, if ω( Û ψ, Û φ) = ω(ψ, φ).In classical mechanics, conservation of the symplectic two-form has many farreaching consequences, one of which is Liouville's theorem-the conservation of phase space volume.In quantum mechanics, a symplectic two-form can be defined as 25 ω(ψ, φ) := −2 Im ψ|φ ; obviously, it is conserved if the inner product ψ|φ itself is.The exact evolution operator is therefore symplectic.
The expectation value of energy is conserved if the evolution operator is unitary and commutes with the Hamiltonian: The exact evolution operator is unitary, and because Û (t) = exp(−i Ĥt/ ) can be Taylor expanded into a convergent series in powers of Ĥ, Û (t) also commutes with Ĥ.As a result, the exact evolution conserves energy.
An adjoint Û (t) * of an evolution operator Û (t) is defined as its inverse taken with a reversed time step: An evolution operator is said to be symmetric if it is equal to its own adjoint: 27 An evolution is time-reversible if a forward propagation for time t is exactly cancelled by an immediately following backward propagation for the same time, i.e., if 27 Time reversibility in quantum dynamics is, therefore, a direct consequence of symmetry.The exact evolution operator is both symmetric and time-reversible because Û (t) * = exp(−i Ĥt/ ).
Finally, the time evolution is said to be: (i) stable 38,49,50 if for every > 0, there exists δ( ) > 0 such that (ii) attracting 49,50 if there exists a δ > 0 such that (iii) asymptotically stable if it is both stable and attracting.
These conditions are visualized in Fig. 1.The exact evolution operator is stable but not asymptotically stable because, due to norm conservation, Stable Asymptotically stable Unstable FIG. 1. Schematic representation of stability conditions in Euclidean space; the distance between corresponding points on the two curves (e.g., the tips of the arrows) is analogous to a metric ψ(t) − φ(t) in the Hilbert space; the dotted lines represent ψ(0) − φ(0) .

B. Loss of geometric properties by approximate methods
In approximate propagation methods, the state ψ(t + ∆t) at time t + ∆t, where ∆t is the numerical time step, is obtained from the state ψ(t) at time t by applying an approximate time evolution operator Ûappr (∆t).This operator is in the explicit Euler method and Ûimpl (∆t in the implicit Euler method.Both Euler methods are of the first order in the time step ∆t, and both are neither unitary nor symplectic.Due to their lack of unitarity, the methods do not conserve energy, even though their evolution operators commute with the Hamiltonian.
Neither method is symmetric; in fact, they are adjoints of each other.Hence, neither method is time-reversible.The explicit Euler method is unstable with the distance between two wavefunctions diverging, whereas the implicit Euler method is asymptotically stable.
The second-order differencing method [39][40][41] recovers symmetry by combining a forward and backward steps of the explicit Euler method: The method can be also obtained directly from the time-dependent Schrödinger equation by using a finite-difference approximation While it is almost as simple as the explicit Euler method to implement, the second order differencing has a higher order of accuracy and, in contrast to the explicit Euler method, it is symmetric, time-reversible, and at least conditionally stable, meaning that it remains stable for sufficiently small time steps ∆t.The second order differencing does not conserve the inner product, norm, energy, or symplectic two-form exactly.Yet, it conserves quantities analogous to the inner product [see Eq.The properties of the different methods are summarized in Table I; a more thorough justification of these properties is provided in Appendix A. Although the explicit and implicit Euler methods are not geometric, composing them in a specific way leads to arbitraryorder integrators that preserve many important geometric properties of the exact solution.
Obviously, the compositions are applicable to systems with non-separable Hamiltonians just like the elementary methods themselves.

C. Recovery of geometric properties by composed methods
Composing the explicit and implicit Euler methods, each for a time step ∆t/2, yields a symmetric second-order method (see Proposition 7 of Appendix A).Depending on the order of composition, one obtains either the trapezoidal rule Ûtrap (∆t) := Ûimpl (∆t/2) Ûexpl (∆t/2), (18)   or implicit midpoint method The trapezoidal rule is also known as the Crank-Nicolson method, 46 although the latter It is necessary to stress that the geometric properties of the trapezoidal rule and implicit midpoint method are only preserved if the implicit step, which involves solving a set of linear equations, is executed exactly (or, in practice, to machine accuracy).We solved the system of equations using the generalized minimal residual method, [52][53][54] an iterative method based on Arnoldi process. 55,56It was an appropriate choice since the matrix being inverted was not positive-definite, symmetric, skew-symmetric, Hermitian, or skew-Hermitian, and therefore neither conjugate gradient nor minimal residual method was applicable. 54The initial guess for the implicit step was approximated with the explicit Euler method since for small time steps, the solutions from the explicit and implicit Euler methods differ only by (∆t/ ) 2 Ĥ2 |ψ(t) .

D. Symmetric composition schemes for symmetric methods
Recursively composing symmetric methods with appropriately chosen time steps leads to symmetric integrators of arbitrary orders. 27,33,35More precisely, there exist a natural number M and real numbers γ n , n = 1, . . ., M , called composition coefficients, satisfying M n=1 γ n = 1 and such that if Ûp (∆t) is any symmetric integrator of (necessarily even) order p, then is a symmetric integrator of order p+2.The most common composition schemes (see Fig. 2) are the triple jump 33,35,57,58 with M = 3, and Suzuki's fractals 35 with M = 5, The remaining coefficients are obtained from the relation γ M +1−n = γ n , which expresses that both of these are symmetric compositions.
Because each triple jump is formed of three steps while each Suzuki's fractal is composed of five steps, the p th -order integrator obtained using Suzuki's fractals has (5/3) p 2 −1 times more composition steps than the one obtained from the same symmetric second-order method using the triple jump.Therefore, the p th -order method obtained from Suzuki's fractals takes (5/3) p 2 −1 times longer to execute per time step ∆t than does the method of the same order achieved through the triple jump.Yet, the leading order error coefficient of the p th -order integrator based on Suzuki's fractal is smaller because the magnitude of each composition step is smaller in Suzuki's fractal.Consequently, to achieve the same accuracy at a final time t, larger time steps can be typically used for calculations using Suzuki's fractals compared to those based on the triple jump.
Non-recursive composition schemes, which require fewer composition steps and are also more efficient, have been obtained for various specific orders.We will refer to these as "optimal" schemes because they minimize the "magnitude" of composition steps.The magnitude of composition steps can be defined as either max n |γ n | or M n=1 |γ n |.With either definition, Suzuki's fractal is the optimal fourth-order scheme.The optimal sixth-and eighth-order schemes, 59 found by Kahan and Li by minimizing max n |γ n |, have two more composition steps (M = 9 and M = 17, respectively) than the minimum number possible for the respective order; the optimal tenth-order scheme, 60  Proof.We prove the theorem in much greater generality.Indeed, a composition of any unitary operators Û1 and Û2 is unitary since A composition of any symplectic operators is symplectic since However, a composition of two symmetric operators is, in general, not symmetric: It is symmetric if the two operators commute or if it is a symmetric composition, e.g., Finally, a composition of time-reversible operators is not necessarily time-reversible since The composition is time-reversible if the two operators commute or if it is a symmetric composition, e.g.,

E. Dynamic Fourier method
To propagate the wavepacket using the explicit or implicit Euler method, or one of their compositions (see Sec. II B-II D), only the action of the Hamiltonian operator Ĥ on ψ(t) is required provided that the implicit steps are solved iteratively.The dynamic Fourier method 31,32,40,61 is an efficient approach to compute f (x)ψ(t), where f (x) is an arbitrary function of x, which denotes either the nuclear position (q) or momentum (p) operator.
Each action of f (x) on ψ(t) is evaluated in the x-representation (in which x is diagonal) by a simple multiplication, after Fourier-transforming ψ(t) to change the representation if needed.On a grid of N points, f (x)ψ(t) is evaluated as f (x i )ψ(x i , t), 1 ≤ i ≤ N , where ψ(x, t) is the wavepacket in the x-representation and x i are either the position or momentum grid points.

F. Molecular Hamiltonian in the adiabatic basis
The molecular Hamiltonian in the adiabatic basis can be expressed as where m is the diagonal D × D nuclear mass matrix, D the number of nuclear degrees of freedom, V the diagonal S×S potential energy matrix, S the number of considered electronic states, and F the nonadiabatic coupling vector (more precisely, a D-vector of S×S matrices).
In Eq. ( 22), the dot • denotes the matrix product in nuclear D-dimensional vector space, the hat ˆrepresents a nuclear operator, and the bold font indicates an electronic operator, i.e., an S × S matrix.Using the dynamic Fourier method, each evaluation of the action of Ĥ on a molecular wavepacket ψ(t), which now becomes an S-component vector of nuclear wavepackets (one on each surface), involves 4D changes of the wavepacket's representation.
In two-state models (i.e., for S = 2), it is possible to obtain Ĥ in the adiabatic representation analytically from the one in the diabatic representation, 62-64 in which W(q) is the (real) diabatic potential energy matrix and in which the nonadiabatic vector couplings vanish.The adiabatic potential energy matrix V(q) is obtained by diagonalizing its diabatic analog W(q), V(q) = O(q) T W(q)O(q), ( and the molecular wavepacket in the adiabatic basis ψ(t) is obtained from its counterpart in the diabatic basis ψ diab (t) as using an orthogonal matrix with ].The two adiabatic energies are given by where ∆W := W 22 − W 11 and W := (W 11 + W 22 )/2.Finally, the transformed momentum operator is By comparing with Eq. ( 22), we see that, in the adiabatic basis, the nonadiabatic coupling vector satisfies F(q) = O(q) T O (q); in particular, F 11 (q) = F 22 (q) = 0, This allowed us to compare the proposed integrators, applied in the adiabatic basis, with the split-operator algorithm, which can only be used in the diabatic basis.Such a rigorous comparison would only be possible for a two-surface model potential because the split-operator algorithm requires the diabatization of the Hamiltonian formulated in the adiabatic representation and this cannot be done exactly for higher-dimensional ab initio potential energy surfaces with more electronic states.
Before the electronic excitation, the NaI molecule was assumed to be in the ground vibrational eigenstate of a harmonic fit to the ground-state potential energy surface at the equilibrium geometry.This vibrational wavepacket was then lifted to the excited-state surface, in order to obtain an initial Gaussian wavepacket (q 0 = 4.9889 a.u., p 0 = 0 a.u., σ 0 = 0.110436 a.u.) for the nonadiabatic dynamics.This use of the sudden approximation assumes an impulsive excitation, i.e., the simultaneous validity of the time-dependent perturbation theory and Condon and ultrashort pulse approximations during the excitation process.After that, the nonadiabatic dynamics was performed by solving the time-dependent Schrödinger equation using the dynamic Fourier method (see Sec. II E) on a uniform grid with 2048 points between q = 3.8 a.u. and q = 47.0 a.u.; Appendix B shows wavepacket represented on such a grid is converged for the duration of the dynamics.A long-enough propagation time, t f = 10500 a.u., was chosen so that the wavepacket traverses the avoided crossing and simultaneously witnesses the change of the nature of the excited adiabatic state from covalent to ionic.The top panel of Fig. 3 shows the two adiabatic potential energy surfaces as well as the initial wavepacket at t = 0 and the final wavepacket at t = t f .The population dynamics of NaI, displayed in the bottom panel of Fig. 3, shows that after crossing the region of highest nonadiabatic coupling, most of the wavepacket remains in the bound excited adiabatic state, while a small population transfer occurs to the dissociative ground state.The figure also confirms that the converged populations obtained with different integrators agree on the scale visible in the figure; in particular, the results obtained with integrators designed for the adiabatic basis agree with each other and also with the result of the split-operator algorithm in the diabatic basis.
To compare various integrators quantitatively, it is essential to "zoom in" and inspect the convergence error at the final time t f ; after all, the dynamic Fourier method 31,32,40 is expected to describe the wavepacket with a high degree of accuracy.In our setting, the convergence error at time t is defined as the L 2 -norm error ψ ∆t (t) − ψ ∆t/2 (t) , where ψ ∆t (t) denotes the wavepacket propagated with a time step ∆t.We omit the split-operator algorithm, which served as a benchmark in Fig. 3, from the following analysis because this algorithm is not applicable to time propagation in the adiabatic representation.Note, however, that 0.95 1.0 P 2

SOD Trapezoidal
Midpoint Split-operator 0 2500 5000 7500 10000 t (a.u.) 0.0 0.05 3. Nonadiabatic dynamics of NaI.Top: Adiabatic potential energy surfaces with the initial and final nuclear wavepacket components in the ground and excited adiabatic electronic states [Because the initial molecular wavepacket was in the excited state, its component ψ 1 (q, t = 0) ≡ 0 is not shown.].Bottom: Ground-and excited-state populations of NaI computed with four different second-order methods: SOD stands for the second-order differencing.The populations were propagated with a time step ∆t = 0.01 a.u., i.e., much more frequently than the markers suggest.The small time step guaranteed wavepacket convergence errors below ≈ 10 −6 in all methods.
for separable Hamiltonians, such as the nonadiabatic Hamiltonian in the diabatic basis, the split-operator algorithms are more efficient than the present integrators of the corresponding order (see Table I and Paper II).
Figure 4 plots the convergence error as a function of the time step and confirms, for each algorithm, the asymptotic order of convergence predicted in Sections II B-II D. Recall that the trapezoidal rule and implicit midpoint method are obtained by composing the explicit and implicit Euler methods, and that the higher order methods are obtained from the trapezoidal rule or implicit midpoint method using the triple jump, Suzuki's fractal, or optimal composition.The top panel of Fig. 4 compares the convergence of all methods, while, for clarity, the bottom left-hand panel only compares the different orders of composition for the Suzuki's fractal and the bottom right-hand panel compares different composition schemes with the sixth order of convergence.(In Fig. 4 and all following figures, we have omitted the results of the implicit midpoint method and of its compositions because they overlap almost perfectly with the corresponding results for the trapezoidal rule.)It is clear that, for a given order of convergence, the prefactor of the error is the largest for the triple jump composition, 33,35 intermediate for the optimally composed 59,60 method, and smallest for Suzuki's fractal 35 composition.To guarantee the correct order of convergence of all composed methods, the composed elementary second-order method must be exactly symmetric, which requires that the systems of linear equations arising from implicit steps must be solved numerically exactly.

Composition type
No comp Suzuki Optimal Triple J.In Section II B, we mentioned the instability of the explicit Euler method and the conditional stability of the second-order differencing.Both properties are reflected in the top panel of Fig. 4 in the divergence of the errors of the two methods for large time steps.The critical time step for the second-order differencing is ∆t ≈ 0.5 a.u., whereas the explicit Euler method is unstable regardless of ∆t but the effect of instability is more visible for larger time steps.In contrast, the trapezoidal rule, implicit midpoint method, and their compositions are stable, but implicit, and, therefore, require the solution of systems of linear equations.

Explicit E SOD Implicit E Trapezoidal
These methods could not be used beyond a certain time step (max n |γ n |∆t ≈ 100 a.u.for both the trapezoidal rule and implicit midpoint method) because the iterative generalized minimal residual algorithm did not converge for very large ∆t.
Convergence of the wavepacket's phase, which is very important, e.g., in the evaluation of spectra, is shown in Fig. 5.As a measure of the convergence error of the phase, we use , where ϕ ∆t := arg[ψ ∆t (q max , t f )] and q max := arg max q (|ψ ∆t (q, t f )|), i.e., ϕ ∆t is the phase of the wavefunction propagated with time step ∆t at the position q max , for which the amplitude of the wavefunction achieves its maximum.Note that the order of convergence is identical to that of the wavepacket.FIG. 5. Convergence of the phase as a function of the time step.All higher-order integrators were obtained through Suzuki's fractal composition.Gray lines were added to guide the eye.
The efficiency of an algorithm cannot be judged solely from the convergence error for a given time step ∆t because the number of composition steps depends on the composition scheme, and, indeed, grows exponentially for the triple-jump and Suzuki's fractal compositions.Figure 6, therefore, displays the convergence error ψ ∆t (t) − ψ 0 (t) , where ψ 0 (t) is the wavepacket propagated using the optimally composed 10 th -order trapezoidal rule with an infinitesimal time step (in practice, ∆t = 0.01 a.u.), as a function of the computational (CPU) time.Among the elementary first-and second-order methods, compared in the top panel of Fig. 6, the second-order differencing, which does not require the solution of a system of linear equations, is the most efficient.Comparison of the geometric integrators based on the trapezoidal rule in the middle and bottom panels of Fig. 6 shows that the fourth-order Suzuki composition takes less CPU time to achieve convergence error as high as 10 −2 than does the elementary trapezoidal rule (i.e., Crank-Nicolson method).To reach errors below 10 −2 , it is already more efficient to use the higher-order integrators.For a more dramatic example, note that the CPU time required to reach an error of 10 −10 is roughly 1000 times longer for the original trapezoidal rule than for its optimal 8 th -order composition.(In Paper II, we confirm that this gain in efficiency holds in higher dimensions by applying the compositions of the trapezoidal rule to the nonadiabatic dynamics in a three-dimensional model of pyrazine in the diabatic representation. 48) The bottom panel of Fig. 6 confirms the prediction that the optimal compositions are the most efficient among composition methods of the same order.Finally, note that the dependence of CPU time on the error in Fig. 6 is not monotonous for the integrators with implicit steps because the convergence of the numerical solution to the system of linear equations required more iterations for larger time steps; as a result, both the error and CPU time increased for time steps larger than a certain critical value (see Fig. 6).
Besides increased efficiency, another benefit of the algorithms based on the composition of the trapezoidal rule or implicit midpoint method is the conservation of the geometric properties of the exact evolution operator.Conservation of the energy, norm, symplectic twoform, and time reversibility by the trapezoidal rule and their compositions is demonstrated in Fig. 7. Time reversibility is measured by the distance of an initial state ψ(0) from Û (−t) Û (t)ψ(0), i.e., a state propagated first forward in time for time t and then backward in time for time t.The tiny residual errors (< 2 • 10 −12 in all cases) of the invariants result from accumulated numerical errors of the fast Fourier transform and generalized minimal residual algorithm.In contrast, the second-order differencing conserves energy, norm, and symplectic two-form only approximately with much larger, O(∆t 4 ) errors (see Propositions 5 and 6 of Appendix A).Although the second-order differencing is time-reversible in theory, its practical implementation is not.[For the second-order differencing to be exactly timereversible, the wavepackets at time t = 0 and t = −∆t would have to be known exactly before the start of the simulation.However, because only ψ(0) is typically available, ψ(−∆t) must be approximated with explicit methods such as the second-order Runge-Kutta scheme.

Composition type
No comp Suzuki Optimal Triple J.

Elementary method
Explicit E SOD Implicit E Trapezoidal FIG. 6. Efficiency of various methods shown using the dependence of the convergence error on the computational (CPU) time.Results of the elementary trapezoidal rule were extrapolated using the line of best fit to highlight the speedup achieved with higher-order compositions.Top: elementary methods; middle: trapezoidal rule and its compositions; bottom: detail of the middle panel.

Composition type
No comp Suzuki Optimal Triple J. by the Euler and second-order differencing methods (top panels), elementary second-order methods (middle panels), and composed methods (bottom panels).φ(0) is a complex Gaussian wavepacket with q 0 = 5.05 a.u., p 0 = 2.5 a.u., and σ 0 identical to that of ψ(0).As a consequence, ω(ψ(0), φ(0)) is nonzero.For the Euler and second-order differencing methods the time step was ∆t = 0.5 a.u. to ensure the stability of the second-order differencing.For all other methods, ten times larger time step (∆t = 5 a.u.) was used to highlight that the exact conservation of invariants is independent of the accuracy of the wavefunction itself.

Elementary method
addition, are exactly unitary, symplectic, and time-reversible.We have shown that unlike the original trapezoidal rule or implicit midpoint method, which are only second-order, their recursive symmetric compositions can achieve accuracy of arbitrary even order in the time step.
We have proven all these properties analytically as well as demonstrated them numerically on a two-surface model of NaI photodissociation.As expected, the higher-order integrators significantly sped up calculations when higher accuracy was required.For example, even to achieve a moderate wavefunction convergence error of 10 −5 , tenfold reduction in the computational time was observed by using higher-order methods compared to the elementary trapezoidal rule.It is probable that Chebyshev 43,65 and short iterative Lanczos schemes 44,45 the leading order local error of both methods is i( Ĥ) 3 /12.
Finally, the local error of the second-order differencing method is which is found by Taylor expanding ψ sod (t − ), assumed to be exact, in Eq. ( 16) to obtain Subtracting Eq. (A5) from Eq. (A10) gives the local error (A8).
The analysis of its geometric properties is simplified if the second-order differencing is found using det Ûsod ( ) = 1, shows that the second-order differencing is not unitary.
Although the second-order differencing is not unitary, a conserved quantity analogous to the inner product exists: (A21)], norm, 40,41 energy 40,41 [see Eq. (A32)], and symplectic two-form [see Eq. (A28)].The corresponding exact quantities are conserved only up to the fourth order in ∆t (see Propositions 5 and 6 of Appendix A).
obtained by Sofroniou and Spaletta by minimizing M n=1 |γ n |, has four more (M = 35).Theorem.All compositions of the trapezoidal rule or implicit midpoint method are unitary, symplectic, stable, energy-conserving, and their evolution operators commute with the Hamiltonian; all symmetric compositions are symmetric and therefore time-reversible.

FIG. 2 . 3 p 2 − 1 5 p 2 − 1
FIG. 2. Pictorial representation of recursive (triple jump and Suzuki's fractals) and non-recursive "optimal" composition schemes.The triple jump has 3 p 2 −1 composition steps per time step, where p is the order of the method, whereas Suzuki's fractal has 5 p 2 −1 composition steps.

)
III. NUMERICAL EXAMPLESTo test the geometric and convergence properties of the integrators presented in Sections II B-II D, we used these integrators to simulate the nonadiabatic quantum dynamics in a two-surface model47 of the NaI molecule.This one-dimensional model, motivated by the experiment by Mokhtari et al.,3 has two electronic states, and therefore an analytical transformation between diabatic and adiabatic representations is available (see Sec. II F).

FIG. 4 .
FIG. 4. Convergence of the molecular wavefunction as a function of the time step.Gray lines were added to guide the eye.Top: all discussed methods; bottom left: methods composed through Suzuki's fractals, bottom right: sixth-order methods.

represented by a 2 × 2
conjugate Ûsod ( ) † of Ûsod ( ) with its inverse, implies a second-order finite-difference approximation to the spatial derivative in the kinetic energy operator, whereas we use the dynamic Fourier method (see Sec. II E), which has exponential convergence with grid density (see Appendix B).Both the trapezoidal rule and implicit midpoint methods are Cayley transforms 51 of (i∆t/2 ) Ĥ and, therefore, are unitary; in addition, both are second-order, symplectic, symmetric, time-reversible, and stable regardless of the size of the time step.Both methods also commute with the Hamiltonian, are energy conserving, and can be further recursively composed to obtain arbitrary-order methods (see Sec. II D).The summary of the properties is given in TableIand a detailed justification provided in Appendix A.
35,35 I. Geometric properties and computational cost of various integrators.Cost is measured by the number of Fourier transforms required per time step (see Sec. II F).I is the number of iterations for the implicit step and n = 0, 1, 2, . . . is the number of recursive compositions.C is the total number of composition steps per time step (C = 3 n for the triple jump33,35, C = 5 n for Suzuki's fractals35).+denotes that the geometric property of the exact evolution operator is preserved and − denotes that it is not.SOD stands for the second-order differencing and SO for the split-operator algorithm.aStability holds for time steps that satisfy Eq. (A48).frequently that also the wavefunction will have an error increasing beyond any bound.As for the implicit Euler method, its error of the norm converges to −1 because ψ impl (t) → 0 as t → ∞ [see Fig.7(b), top panel].We have described geometric integrators for nonadiabatic quantum dynamics in the adiabatic representation, in which the popular split-operator algorithms cannot be used due to nonseparability of the Hamiltonian into a kinetic and potential terms.The proposed methods are based on the symmetric composition of the trapezoidal rule or implicit midpoint method, and as a result, are symmetric, stable, conserve the energy exactly and, in 40] None of the four geometric properties or analogous quantities is conserved by the Euler methods.The explicit Euler method is unstable regardless of ∆t, and will, for long enough times t f , result in a norm divergent to infinity [see Fig.7(b), top panel] even for very small ∆t, implying