Efficient geometric integrators for nonadiabatic quantum dynamics. II. The diabatic representation

Exact nonadiabatic quantum evolution preserves many geometric properties of the molecular Hilbert space. In a companion paper [S. Choi and J. Van\'{\i}\v{c}ek, 2019], we presented numerical integrators of arbitrary-order of accuracy that preserve these geometric properties exactly even in the adiabatic representation, in which the molecular Hamiltonian is not separable into a kinetic and potential terms. Here, we focus on the separable Hamiltonian in diabatic representation, where the split-operator algorithm provides a popular alternative because it is explicit and easy to implement, while preserving most geometric invariants. Whereas the standard version has only second-order accuracy, we implemented, in an automated fashion, its recursive symmetric compositions, using the same schemes as in the companion paper, and obtained integrators of arbitrary even order that still preserve the geometric properties exactly. Because the automatically generated splitting coefficients are redundant, we reduce the computational cost by pruning these coefficients and lower memory requirements by identifying unique coefficients. The order of convergence and preservation of geometric properties are justified analytically and confirmed numerically on a one-dimensional two-surface model of NaI and a three-dimensional three-surface model of pyrazine. As for efficiency, we find that to reach a convergence error of 10$^{-10}$, a 600-fold speedup in the case of NaI and a 900-fold speedup in the case of pyrazine are obtained with the higher-order compositions instead of the second-order split-operator algorithm. The pyrazine results suggest that the efficiency gain survives in higher dimensions.


I. INTRODUCTION
The celebrated Born-Oppenheimer approximation 1,2 assumes the separability of the nuclear and electronic motions in a molecule, and provides an appealing picture of independent electronic potential energy surfaces.5][6][7] To investigate such processes, one can abandon the Born-Oppenheimer representation and treat electrons and nuclei explicitly, [8][9][10] use an exact factorization 11,12 of the molecular wavefunction, or, determine which Born-Oppenheimer states are coupled strongly 13,14 and then solve the time-dependent Schrödinger equation with a nonadiabatically coupled molecular Hamiltonian; below, we will only consider the third and most common strategy.
In a companion paper 15 (which will be referred to as Paper I), we surveyed several algorithms for the nonadiabatic quantum dynamics, applicable to higher dimensions, including Gaussian basis methods, [16][17][18][19][20][21] variations of the multiconfigurational time-dependent Hartree (MCTDH) method, [22][23][24] and sparse-grid methods. 25,26There are situations, however, in which the wavepacket spreads over large parts of the available Hilbert space, and then time-independent basis sets or full-grid methods can become more efficient.
As for the molecular Hamiltonian used in nonadiabatic simulations, the ab initio electronic structure methods typically yield the adiabatic potential energy surfaces, which are nonadiabatically coupled via momentum couplings.However, in the regions of conical intersections, 27,28 the Born-Oppenheimer surfaces become degenerate, and the nonadiabatic couplings diverge.To avoid associated problems, it is convenient to use the diabatic representation, in which the divergent momentum couplings are replaced with well-behaved coordinate couplings.Although exact diabatization is only possible in systems with two electronic states and one nuclear degree of freedom, 29 there exist more general, approximate diabatization procedures, [30][31][32] starting with the vibronic coupling Hamiltonian model. 33Another benefit of the diabatic representation is that it separates the Hamiltonian into a sum of kinetic energy, depending only on nuclear momenta, and potential energy, depending only on nuclear coordinates, which makes it possible to propagate the molecular wavefunction with the split-operator algorithm. 26,34,35The split-operator algorithm is explicit, easy to implement, and, in addition, it is an example of a geometric integrator 36,37 because, similarly to the integrators discussed in Paper I, 15 it conserves exactly many invariants of the exact so-lution, regardless of the convergence error of the wavefunction itself.Geometric integrators in general acknowledge special properties of the Schrödinger equation which differentiate it from other differential equations.Using these integrators can be likened to using a wellfitting screw-driver instead of a hammer to attach a screw.Note that the integrators for nonseparable Hamiltonians, presented in Paper I, are also geometric and, clearly, still applicable to the separable Hamiltonian in the diabatic representation, but the split-operator algorithm is expected to be more efficient because it is explicit.
The standard, second-order split-operator algorithm 34 is unitary, symplectic, stable, symmetric, and time-reversible, regardless of the size of the time step.However, to obtain highly accurate results, the standard algorithm requires using a small time step, because it has only second-order accuracy.There exist much more efficient algorithms, such as the shortiterative Lanczos algorithm, [38][39][40] which has an exponential convergence with respect to the time step, and also conserves the norm and energy, but not the inner product (because it is nonlinear) and other geometric properties.
To address the low accuracy of the second-order split-operator algorithm and the nonconservation of geometric properties by other more accurate methods, various higher-order splitoperator integrators have been introduced, 41-44 some of which allow complex time steps [44][45][46] or commutators of the kinetic and potential energies in the exponent, [47][48][49] thus reducing the number of splitting steps.Here we explore one type of higher-order integrators, designed for nonadiabatic dynamics in the diabatic basis, which we have implemented using the recursive triple-jump 42,43 and Suzuki-fractal, 42 as well as several non-recursive, "optimal" compositions of the second-order split-operator algorithm.While the recursive compositions permit an automated generation of integrators of arbitrary even order in the time step, 36,37,42,43,50,51 the efficiency of higher-order algorithms is sometimes questioned because the number of splitting steps grows exponentially with the order of accuracy, and, consequently, so does the computational cost of a single time step.Motivated by this dilemma, we have explored the convergence and efficiency of the higher-order compositions using a one-and three-dimensional systems, concluding that, despite the increasing number of splittings, the higher-order methods become the most efficient if higher accuracy of the solution is required, and that this gain in efficiency survives in higher dimensions.We have also confirmed that all composed methods are unitary, symplectic, stable, symmetric, and timereversible.A final benefit of the higher-order methods is the simple, abstract, and general implementation of the compositions of the second-order split-operator algorithm; indeed, even this "elementary" method is a composition of simpler, first-order algorithms. 26,35e of the only challenges of implementing the split-operator algorithm for nonadiabatic dynamics in the diabatic representation is the exponentiation of the potential energy operator, which is nondiagonal in the electronic degrees of freedom (in contrast to the diagonal kinetic energy operator).We, therefore, explored several methods for the exponentiation of nondiagonal matrices.
The main disadvantage of the split-operator algorithm and its compositions is that their use is restricted to separable Hamiltonians.To compare them with the integrators from Paper I, we cannot use the adiabatic representation, but instead must perform the comparison in the diabatic representation, where the compositions of the explicit split-operator algorithm are, as expected, much more efficient than the more generally applicable compositions 15 of the implicit trapezoidal rule (the Crank-Nicolson method 52,53 ) from Paper I. Nevertheless, the comparison serves as a higher-dimensional test of integrators from Paper I and confirms that, in contrast to the split-operator compositions, the integrators from Paper I conserve also the energy exactly.
The remainder of this paper is organized as follows: In Sec.II, after reviewing the geometric properties of the exact evolution operator, we discuss the lack of symmetry and time-reversibility in the first-order split-operator algorithms and the recovery of these properties in the symmetric compositions.Next, we describe several strategies for reducing the computational cost and memory requirements by pruning redundant splitting coefficients generated automatically by the symmetric compositions.After presenting the dynamic Fourier method for its ease of implementation and the exponential convergence with the grid density, we briefly discuss the molecular Hamiltonian in diabatic representation.In Section III, the convergence properties and conservation of geometric invariants by various methods are analyzed numerically on a one-dimensional two-surface model 54 of NaI and a three-dimensional three-surface model of pyrazine, 55 both in the diabatic representation.
Section IV concludes the paper.

II. THEORY A. Geometric properties of the exact evolution operator
The time-dependent Schrödinger equation with a time-independent Hamiltonian Ĥ and initial condition ψ(0) has the formal solution , where Û (t) is the evolution operator.While in Paper I, we considered general Hamiltonian operators Ĥ ≡ H(q, p), here we require that the Hamiltonian be separable as into a sum of kinetic and potential energies, which depend, respectively, only on the momentum p and position q operators.
The exact evolution operator is linear, unitary, symplectic, symmetric, time-reversible, stable, and conserves the norm, inner product, and energy.Because these properties are desirable also in approximate numerical evolution operator Ûappr (t), let us define them briefly.
An operator Û is said to preserve the norm if Û ψ = ψ for all ψ, and to preserve the inner product if Û ψ| Û φ = ψ|φ for all ψ and φ.For linear operators Û , these two properties are equivalent, whereas for general, possibly nonlinear operators, conservation of the inner product implies linearity 56 and hence the conservation of norm, but norm conservation implies neither linearity nor conservation of the inner product.An operator Û is said to be where Ûappr (∆t) is an approximate time evolution operator and ∆t the numerical time step.
Depending on the order of kinetic and potential propagations, the approximate evolution operator is in the VT split-operator algorithm and in the TV split-operator algorithm.Both ÛVT and ÛTV are unitary, symplectic, stable, but only first-order in the time step ∆t.Neither method conserves energy because neither evolution operator commutes with the Hamiltonian.Neither method is symmetric; in fact, they are adjoints of each other.Hence, neither method is time-reversible.These properties are justified in Appendix A and summarized in Table I.
Although the first-order split-operator algorithms are not time-reversible, composing them in a specific way leads to time-reversible integrators of arbitrary order of accuracy in the time step.

C. Recovery of geometric properties by composed methods
Composing the two first-order split-operator algorithms, each for a time step ∆t/2, yields a symmetric second-order method. 34Depending on the order of composition, one obtains either the VTV algorithm or TVT algorithm ÛTVT (∆t) := ÛTV (∆t/2) ÛVT (∆t/2).( 7) Both are explicit, unitary, symplectic, stable, symmetric, and time-reversible, regardless of the size of the time step.Neither evolution operator commutes with the Hamiltonian and, therefore, neither method conserves energy exactly.These properties are again justified in Appendix A and summarized in Table I.

D. Symmetric composition schemes for symmetric methods
As discussed in Paper I, composing any symmetric second-order method (such as one of those of Sec.II C) with appropriately chosen time steps leads to symmetric integrators of arbitrary order of accuracy. 36,37,42,43More precisely, there are a natural number M and and such that for any symmetric evolution operator Ûp (∆t) of an even order p, composing this symmetric evolution operator with coefficients γ n yields a symmetric integrator of order The simplest composition schemes (see Fig. n is the number of recursive compositions and C the total number of composition steps per time step (C = 3 n for the triple jump 42,43 , C = 5 n for Suzuki's fractal 42 ).+ or − denotes that the geometric property of the exact evolution operator is or is not preserved.
Method Order Unitary Symplectic Commutes Energy Symm-Time-Stable Cost with Ĥ cons.etric reversible this composition is sometimes more efficient than the triple-jump composition, despite requiring more composition steps (see Ref. 15 for a numerical example).For specific orders of convergence, more efficient non-recursive composition schemes exist and will be referred to as "optimal."These were implemented according to Kahan and Li 58 for the 6 th and 8 th orders, and according to Sofroniou and Spaletta 59 for the 10 th order (see Sec. II D of Paper I for more details about composition methods).

E. Compositions of split-operator algorithms
The split-operator algorithm is applicable if the Hamiltonian Ĥ can be written as a sum of operators Â and B with evolution operators, Û Â(t) = exp(−it Â/ ) and Û B (t) = exp(−it B/ ), whose actions on ψ can be evaluated exactly.A general split-operator evolution operator can be expressed as where N is the number of splitting steps, and a j and b j are the splitting coefficients associated with the operators Â and B. These coefficients in general satisfy the identity N j=1 a j = N j=1 b j = 1, and are a 1 = b 1 = 1 for the first-order VT and TV algorithms 60 and for the second-order VTV or TVT algorithms. 61cause the second-order split-operator algorithm 61 is symmetric, it can be composed by any of the composition schemes discussed in Sec.II D. For example, the splitting coefficients of a fourth-order method are with N = 6 if the triple-jump composition scheme is used, and with N = 10 if Suzuki's fractal is used instead.The remaining coefficients are obtained from symmetry as Both composition procedures can be applied recursively to obtain higher-order split-operator algorithms.These as well as the optimally composed algorithms of up to the tenth order are represented pictorially in Fig. 1. and nonrecursive "optimal" composition schemes shown in Fig. 2 of Paper I. 15 In other words, each elementary method Û (γ n ∆t) (solid line segment in Fig. 2 of Paper I) is replaced by a secondorder split-operator algorithm Û Â(γ n ∆t/2) Û B (γ n ∆t) Û Â(γ n ∆t/2), represented here by a triple of consecutive solid, dotted, and solid line segments.Solid line segments represent Û Â(γ n ∆t/2), whereas the dotted line segments represent Û B (γ n ∆t).N Ô is the number of actions of Û Ô on ψ.
All compositions of the second-order VTV or TVT split-operator algorithms are unitary, symplectic, and stable; all symmetric compositions are symmetric and, therefore, timereversible.The proof of this statement is a special case of the general proof of a corresponding theorem for the composition of geometric integrators in Paper I.

F. Pruning splitting coefficients
Many b j coefficients of the higher-order integrators obtained by recursive composition of the second-order split-operator algorithm are zero [for an example, see Eqs. (10) and (11)].
The computational time can be reduced by "pruning," i.e., removing the splitting steps corresponding to b j = 0 and merging the consecutive actions of Û Â(a j ∆t) and Û Â(a j+1 ∆t).
If b j = 0 and j = N , the splitting coefficients are modified as in order to merge the j th and (j + 1) th steps.The composed methods after the merge are exhibited in Fig. 2 and the reduction in the number N of splitting steps, which measures the computational cost, is summarized in Table II.
For a time-independent separable Hamiltonian, one can either precompute and store the evolution operators, Û Â(a j ∆t) and Û B (b j ∆t), or compute them on the fly.The former approach is more memory intensive than the latter, which does not store any evolution operators, but the computational time is reduced since the evolution operators are only computed once at initialization.To alleviate the memory requirement of the former approach, one can exploit the repetition of certain splitting coefficients, which is obvious from Eqs. ( 10) and (11) and Fig. 2. If either Â or B is time-dependent, it is always beneficial to compute the corresponding evolution operator pertaining to the time-dependent operator on the fly because no reduction in computational time is possible by precomputing the evolution operators.
The effort spent in searching for repeated coefficients is reduced if the symmetries of the composition scheme and of the elementary method are exploited [see Eq. ( 12)].The repeated coefficients are then identified from only half of the original coefficients a j and b j .
Once identified, only the unique evolution operators Û Â(a unq j ∆t) and Û B (b unq j ∆t) are stored in arrays of lengths N unq coefficient arrays, i.e., Exploiting the repeated coefficients, the number of stored evolution operators reduces from II).

G. Dynamic Fourier method
To propagate a wavepacket ψ(t) with any split-operator algorithm (see Secs.II B-II D), only the actions of the kinetic ( Û T ) and potential ( Û V ) evolution operators on ψ(t) are required, where Û T (∆t) := e −i∆tT (p)/ and Û V (∆t) := e −i∆tV (q)/ .
Since Û T and Û V are diagonal in the momentum and position representations, respectively, their action on ψ(t) is easy to evaluate in the appropriate representation.This is the main idea of the dynamic Fourier method, 34,35,62,63 in which the representation of ψ(t) is repeatedly changed, as needed, via the fast Fourier transform (for more details, see Sec.II E of Paper I).
In the numerical examples below, the Fourier transform was performed using the Fastest Fourier Transform in the West 3 (FFTW3) library. 64Although its accuracy is sufficient for a N Â = 2N B for order ≥ 2.
most applications, small deviations from unitarity, which were due to the high number of repeated application of the forward and backward Fourier transforms, affected the most converged calculations.To reduce the nonunitarity, we used the long-double instead of the default double precision version of FFTW3.

H. Molecular Hamiltonian in the diabatic basis
The molecular Hamiltonian in the diabatic basis can be expressed as where m is the diagonal D × D nuclear mass matrix, D the number of nuclear degrees of freedom, and V the potential energy.In Eq. ( 15), 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, where S is the number of included electronic states.Using the dynamic Fourier method, each evaluation of the action of the pair Û V(t V ) and ÛT (t T ) on a molecular wavepacket ψ(t), which now becomes an Scomponent vector of nuclear wavepackets (one on each surface), involves two changes of the wavepacket's representation.The above-mentioned nonunitarity of the solution, partially due to the numerical implementation of the FFT algorithm, was made worse by the matrix exponential required for evaluating the potential evolution operator Û V(t V ), which contains offdiagonal couplings between the electronic states.Although we tried different approaches for matrix exponentiation, including Padé approximants 65,66 and exponentiating a diagonal matrix obtained with the QR decomposition 65,67 or with the Jacobi method, 65 none of the three methods was better than the others in reducing the nonunitarity.Since both in the NaI and pyrazine models, only 2 × 2 matrices are relevant, and since for such matrices, the Jacobi method yields already after one iteration the analytically exact result for the exponential, we used the Jacobi method for all results in Sec.III.Note, however, that the other two methods (based on Padé approximants or QR decomposition), while not exact in the two models used in this paper, converge, in general, faster than the Jacobi method, and are, therefore, preferred in systems with more than two coupled electronic states.

I. Trapezoidal rule and implicit midpoint method
In addition to nonconservation of energy, the main disadvantage of the split-operator algorithms is that they can be applied to nonadiabatic dynamics only in the diabatic representation.Yet, there exist closely related, arbitrary-order geometric integrators, discussed in Paper I, which, in addition, conserve energy and are applicable both in the diabatic and adiabatic representations.These integrators are, like the higher-order split-operator algorithms, based on recursive symmetric composition (see Sec. II D) of the second-order trapezoidal rule (Crank-Nicolson method 52,53 ) or the implicit midpoint method, both of which are, themselves, compositions of the explicit and implicit Euler methods [see Eqs. ( 18), ( 19), (13), and ( 14) of Paper I]. Due to the presence of implicit steps, the trapezoidal rule, implicit midpoint method as well as their compositions require solving large, although sparse, linear systems iteratively, 15 and, as a result, in the diabatic representation are expected to be significantly less efficient than the explicit split-operator algorithms of the same order of accuracy.These integrators are, again, most naturally implemented in conjunction with the dynamic Fourier method described in Sec.II G; the only difference being that one must evaluate the operation ( T + V )ψ instead of Û T ψ and Û V ψ.More details about these higherorder integrators can be found in Paper I, 15 which discusses their geometric properties and studies their efficiency in applications to nonadiabatic quantum dynamics in the adiabatic representation, in which the molecular Hamiltonian is nonseparable.

III. NUMERICAL EXAMPLES
To test the geometric and convergence properties of the split-operator algorithms presented in Sections II B-II D, we used these integrators to simulate the nonadiabatic quantum dynamics in a one-and three-dimensional systems.

A. One-dimensional model of NaI
This model is a diabatized version of the one presented in Paper I, i.e., a one-dimensional two-surface model 54 of the NaI molecule.We used the same initial and final times, and the same approximations for the initial state and for the molecule-field interactions as in Paper I. For detailed calculation parameters, see Section III of Ref. 15.
The top panel of Fig. 3 shows the two diabatic potential energy surfaces as well as the initial wavepacket at t = 0 and the ground-and excited-state components of the final wavepacket at the final time t f = 10500 a.u.The population dynamics of NaI, displayed in the middle and bottom panels of Fig. 3, shows that after passing this crossing, most of the population jumps to the other diabatic state, while a small fraction remains in the original, dissociative diabatic state.On the scale visible in the figure, the converged populations obtained with the VTV and TVT split-operator algorithms agree with each other and also with the results of the trapezoidal and midpoint rule (middle panel).Moreover, the results of the triple-jump, Suzuki-fractal, and optimal compositions of the second-order VTV algorithm agree with each other (bottom panel).
For a quantitative comparison of various algorithms, it is necessary to compare their convergence errors at the final time t f .As in Paper I, the convergence error at time t f as a function of the time step ∆t is measured by the L 2 -norm error ψ ∆t (t f ) − ψ ∆t/2 (t f ) , where ψ τ (t f ) represents the wavepacket propagated with a time step τ .This error is shown in in the adiabatic basis, 15 the prefactor of the error is the largest for the triple-jump, 42,43 intermediate for the optimal, 58 and smallest for Suzuki-fractal composition.The figure also shows that for the smallest time steps, the error starts to increase again.This is due to the accumulating numerical error of the fast Fourier transform, which eventually outweighs the error due to time discretization.As a result, the predicted asymptotic order of convergence cannot be observed for some methods because it is only reached for very small time steps.
While the probability density has a classical analogue, the phase of the wavefunction is a purely quantum property.As a consequence, an accurate evaluation of the phase is very important in the calculation of electronic spectra and in other situations, where quantum effects play a role.To investigate the convergence of the phase as a function of the time step, we used the phase of wavefunction at the maximum of the probability density (for a 10 20 q (a.u.) 0.00 0.25 0.50 precise definition, see Paper I). Figure 5 displays the convergence of the error of the phase for the triple-jump compositions, and confirms that the order of convergence is the same as for the wavefunction itself (bottom left-hand panel of Fig. 4).
Because the number of composition steps depends on the composition scheme and increases with the order, the efficiency of an algorithm is not determined solely by the convergence error for a given time step ∆t.It is, therefore, essential to compare directly the    efficiency of the different algorithms.Figure 6 displays the wavefunction convergence error of each algorithm as a function of the computational (CPU) time.Comparison of the compositions of the VTV split-operator algorithm in the top panel of Fig. 6 shows that the fourth-order Suzuki composition already takes less CPU time to achieve convergence error 10 −2 than does the elementary VTV algorithm.To reach errors below 10 −2 , it is more efficient to use some of the fourth or higher-order integrators.Remarkably, the CPU time required to reach an error of 10 −10 is roughly 600 times longer for the basic VTV algorithm than for its optimal 6 th -order composition.The bottom right-hand panel of Fig. 6 confirms the prediction that the optimal compositions are the most efficient among composition methods of the same order.
Convergence curves in Figs.4-6 were obtained using the long-double precision for the FFTW3 algorithm, which lowered the error accumulation resulting from the nonunitarity of the FFTW3 Fourier transform.If high accuracy is not desired, the double precision of the FFTW3 algorithm can be used instead, resulting in much more efficient higher-order algorithms.This is shown for the NaI model in Fig. 7, which compares the efficiency of the optimal compositions of the VTV algorithm evaluated either with the double or longdouble implementation of the FFTW3, and also with the corresponding compositions of the trapezoidal rule (for which the double precision of FFTW3 was sufficient).Even the more expensive, long-double precision calculation with the compositions of VTV algorithm are faster than the corresponding double precision calculations with the trapezoidal rule, which requires an expensive iterative solution of a system of linear equations.In particular, the sixth-order optimal composition of the VTV algorithm reaches a convergence error of 10 −10 forty times faster than the same composition of the trapezoidal rule (see Fig. 7) and 30000 times faster than the elementary trapezoidal rule (see Figs. 6 and 7).
Note that the dependence of CPU time on the error in Fig. 7 is not monotonous for the compositions of the trapezoidal rule 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.
To check that the increased efficiency of higher-order compositions is not achieved by sacrificing the conservation of geometric invariants, we analyzed, using the NaI model, the conservation of norm, symplectic two-form, energy, and time reversibility.Conservation of the norm and symplectic two-form, and nonconservation of energy by all split-operator al-gorithms is demonstrated in panels (a)-(c) of Fig. 8.The tiny residual errors (< 10 −12 in all cases) result from accumulated numerical errors of the FFT and matrix exponentiation (see Sec. II G).Panels (d) and (e) confirm, on one hand, that the first-order split-operator algorithm is not time-reversible, and, on the other hand, that the second-order VTV algorithm together with all its compositions are exactly time-reversible; the tiny residual errors are again due to accumulated numerical errors of the FFT and matrix exponentiation.
The nonconservation of energy by the split-operator algorithms is further inspected in Fig. 9, showing the error of energy as a function of the time step.For the Suzuki-fractal compositions of the VTV algorithm, the energy is only conserved approximately; its conservation follows the order of convergence of the integrator, as indicated by the gray lines.In contrast, the trapezoidal rule conserves the energy to machine accuracy, regardless of the size of the time step.population dynamics obtained with sixth-order optimal compositions of the VTV algorithm and of the trapezoidal rule agree with each other.for convergence errors below 10 −2 for compositions of the VTV algorithms and, remarkably, already for errors below 10 −1 for compositions of the trapezoidal rule.In particular, to reach an error of 10 −10 , a 900-fold speedup over the second-order VTV algorithm and a 300fold speedup over the second-order trapezoidal rule are achieved by using their tenth-order optimal compositions.These results suggest that increasing the number of dimensions is either beneficial or, at the very least, not detrimental to the gain in efficiency from using the higher-order integrators.As in Fig. 7, the compositions of the VTV algorithms are much more efficient than the compositions of the trapezoidal rule, but this was expected, because the Hamiltonian (15) is separable.One must remember that the main purpose of the compositions of the trapezoidal rule is for nonseparable Hamiltonians, where the split-operator algorithms cannot be used at all.operator algorithm, both double and long-double precision versions are compared.As for the fourth-order algorithms, Suzuki's fractal is considered as the "optimal" composition scheme.

IV. CONCLUSION
We have described geometric integrators for nonadiabatic quantum dynamics in the diabatic representation, in which the Hamiltonian is separable into a kinetic term, depending only on momentum, and potential term, depending only on position.These integrators are based on recursive symmetric composition of the standard, second-order split operator algorithm, and as a result, are explicit, unconditionally stable and exactly unitary, symplectic, symmetric, and time-reversible.Unlike the original split-operator algorithm, which is only second-order, its recursive symmetric compositions can achieve accuracy of an arbitrary even order in the time step.These properties were justified analytically and demonstrated numerically on a diabatic two-surface model of NaI photodissociation.Indeed, the higherorder integrators sped up calculations by several orders of magnitude when higher accuracy 2 of Ref.15) are the triple jump[41][42][43]57 with M = 3, and Suzuki's fractal42 with M = 5. Boh are symmetric compositions, meaning that γ M +1−n = γ n .Because larger time steps can be used for calculations using Suzuki's fractal, TABLE I. Geometric properties and computational cost of the first-order and recursively composed second-order split-operator (SO) algorithms.Cost (here before speedup by pruning splitting coefficients) is measured by the number of fast Fourier transforms required per time step (see Sec. II G).
FIG.1.Split-operator algorithms composed by the recursive (triple jump and Suzuki's fractal)

Fig. 4 ,
Fig. 4, which confirms, for each algorithm, the asymptotic order of convergence predicted in Secs.II B-II D. For clarity, in this and all remaining figures, only the VT algorithm and compositions of the VTV algorithm are compared because the corresponding results of the TV algorithm and compositions of the TVT algorithms behave similarly.The top panel of Fig. 4 compares all methods, whereas the bottom left-hand panel compares only the different orders of the triple-jump composition and the bottom right-hand panel compares only different composition schemes with the sixth-order convergence.Similarly to the results

FIG. 3 .
FIG.3.Nonadiabatic dynamics of NaI.Top: Diabatic potential energy surfaces with the initial and final nuclear wavepacket components in the two diabatic electronic states (the inital groundstate component is not shown because it was zero: ψ 1 (q, 0) = 0).Middle: Populations of NaI in the two diabatic states computed with four different second-order methods.Bottom: Populations computed with three different sixth-order compositions of the VTV algorithm.Populations were propagated with a time step ∆t = 0.01 a.u.for the second-order methods and ∆t = 82.03125a.u.for the sixth-order methods, i.e., much more frequently than the markers suggest.The time step guaranteed wavepacket convergence errors below ≈ 10 −5 in all methods.

FIG. 4 .
FIG. 4. Convergence of the molecular wavefunction as a function of the time step.The wavefunction was propagated with the VT algorithm or with the compositions of the VTV algorithm.Gray straight lines indicate various predicted orders of convergence O(∆t n ).Top: all discussed methods, bottom left: methods composed with the triple-jump scheme, bottom right: sixth-order methods.

FIG. 5 .
FIG. 5. Convergence error of the phase of the wavepacket as a function of the time step for the triple-jump compositions.Gray straight lines indicate various predicted orders of convergence O(∆t n ).

FIG. 6 .
FIG.6.Efficiency of the VT algorithm and of various compositions of the VTV algorithm shown using the dependence of the convergence error on the computational (CPU) time.Top: all methods, bottom left: triple-jump compositions, bottom right: sixth-order methods.The reference wavefunction ψ 0 (t f ) was chosen as the most accurate point in Fig.4, i.e., the wavefunction obtained using the optimal eighth-order composition with a time step ∆t = t f /2 9 .

FIG. 9 .
FIG. 9. Energy conservation as a function of the time step in simulations of the nonadiabatic dynamics of NaI.Gray straight lines indicate various orders of convergence O(∆t n ).

FIG. 10 .
FIG.10.Population dynamics of pyrazine obtained using the sixth-order optimal compositions of the trapezoidal rule and VTV algorithm.The same time step t f /25600 was used for both calculations.

Figure 11 compares
Figure 11 compares the efficiency of different (yet always optimal) compositions of the VTV algorithm and trapezoidal rule.Higher-order integrators become more efficient already

FIG. 11 .
FIG. 11.Efficiency of the optimal compositions of the trapezoidal rule and VTV split-operator algorithm applied to the three-dimensional pyrazine model.For the trapezoidal rule, only the double precision version of the FFTW3 fast Fourier transform was used, while for the VTV split- Im ψ|φ , which is, obviously, conserved, if the inner product is.Û is said to conserve energy if Ĥ Û ψ = Ĥ ψ , where Â ψ := ψ| Â|ψ denotes the expectation value of operator Â in the state ψ.Finally, an adjoint Û (t) * of an evolution operator Û (t) is defined as Û (t) * := Û (−t) −1 .An evolution operator is said to be symmetric if 36,37 Û (t) * = Û (t) and

TABLE II .
Computational cost and memory requirement of the composed split-operator algorithms before and after pruning (i.e., removing zero coefficients and merging adjacent coefficients) and identifying repeated coefficients.The computational cost is measured by N Â + N B , where N Ô is the number of actions of Û Ô on the wavepacket.The memory requirement before and after pruning is N Â + N B , and after identifying repeated coefficients decreases to N unq a + N unq b .