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

Abstract

Geometric integrators of the Schrö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 orderintegrators 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íč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.

An implicit split-operator algorithm for the nonlinear time-dependent Schrödinger equation

“…Therefore, they can be further composed using symmetric composition methods 4,34,[45][46][47][48][50][51][52][53] in order to obtain integrators of arbitrary even orders of convergence. To this end, starting from an integrator Ûp of even order p, an integrator Ûp+2 of order p + 2 is generated using the symmetric composition…”

High-order geometric integrators for representation-free Ehrenfest dynamics

“…Initial energy of the system is E 0 = −0.2 a.u. To avoid clutter, only the higher-order integrators obtained using the optimal composition schemes are shown (the Suzukifractal scheme is the optimal fourth-order scheme 84 ).…”

How important are the residual nonadiabatic couplings for an accurate simulation of nonadiabatic quantum dynamics in a quasidiabatic representation?

“…For a valid comparison of the two wavepackets propagated with either the exact or approximate quasidiabatic Hamiltonian, the spatial and time discretization errors must be much smaller than the errors due to omitting the residual couplings. Owing to its exact symmetry, the implicit midpoint method can be composed using various schemes [44][45][46][47] to obtain integrators of arbitrary even orders of accuracy in the time step; 51,73 we composed the implicit midpoint method according to the optimal scheme 44 to obtain an eighth-order integrator.…”

“…This consideration led us to choose the optimal eighth-order 44 composition [45][46][47][48] of the implicit midpoint method, 47,49,50 which satisfies both requirements and, moreover, preserves exactly geometric properties of the exact solution. 47,49,51 Although all quasidiabatization schemes remove the numerically problematic singularity of the nonadiabatic couplings, the magnitude of the residual couplings depends on the employed scheme. 21 To find out how the errors due to ignoring the residual couplings depend on the sophistication of the quasidiabatization, we, therefore, compare the first-and second-order regularized diabatizations [39][40][41] in the cubic E ⊗ e Jahn-Teller model, in which even the strictly diabatic Hamiltonian exists and both quasidiabatic and strictly diabatic Hamiltonians can be obtained analytically.…”