When simulating nonadiabatic molecular dynamics, choosing an electronic representation requires consideration of well-known trade-offs. The uniqueness and spatially local couplings of the adiabatic representation come at the expense of an electronic wave function that changes discontinuously with nuclear motion and associated singularities in the nonadiabatic coupling matrix elements. The quasi-diabatic representation offers a smoothly varying wave function and finite couplings, but identification of a globally well-behaved quasi-diabatic representation is a system-specific challenge. In this work, we introduce the diabatized Gaussians on adiabatic surfaces (DGAS) approximation, a variant of the *ab initio* multiple spawning (AIMS) method that preserves the advantages of both electronic representations while avoiding their respective pitfalls. The DGAS wave function is expanded in a basis of vibronic functions that are continuous in both electronic and nuclear coordinates, but potentially discontinuous in time. Because the time-dependent Schrödinger equation contains only first-order derivatives with respect to time, singularities in the second-derivative nonadiabatic coupling terms (i.e., diagonal Born-Oppenheimer correction; DBOC) at conical intersections are rigorously absent, though singular time-derivative couplings remain. Interpolation of the electronic wave function allows the accurate prediction of population transfer probabilities even in the presence of the remaining singularities. We compare DGAS calculations of the dynamics of photoexcited ethene to AIMS calculations performed in the adiabatic representation, including the DBOC. The 28 fs excited state lifetime observed in DGAS simulations is considerably shorter than the 50 fs lifetime observed in the adiabatic simulations. The slower decay in the adiabatic representation is attributable to the large, repulsive DBOC in the neighborhood of conical intersections. These repulsive DBOC terms are artifacts of the discontinuities in the individual adiabatic vibronic basis functions and therefore cannot reflect the behavior of the exact molecular wave function, which must be continuous.

## I. INTRODUCTION

Many important physical processes that are initiated by the creation of an electronic excitation involve nonadiabatic dynamics. As such, researchers have spent considerable effort over the last several decades developing theoretical and computational methods to model nonadiabatic processes.^{1–22} When developing a nonadiabatic molecular dynamics simulation protocol, a researcher must choose a representation of the electronic wave functions. With a small number of exceptions,^{23–25} both the efficiency and accuracy of simulations depend on this choice.^{8,26,27} The two most commonly used representations are the adiabatic and quasi-diabatic representations. A less widely adopted but very interesting representation arises from the exact factorization of the total molecular wave function into a single nuclear wave function and a single electronic wave function.^{28–30}

The adiabatic representation offers important advantages, particularly in the context of *ab initio* molecular dynamics simulations. The adiabatic electronic states are those that are obtained straightforwardly by diagonalization of the electronic Hamiltonian. Obviously no coupling between adiabatic states arises from the electronic Hamiltonian, therefore all coupling arises from non-Born–Oppenheimer terms—the nonadiabatic couplings. A convenient feature of these couplings is that they are negligibly small at the majority of the points in nuclear configuration space. In polyatomic molecules, they typically only become large near conical intersections between potential energy surfaces (PESs), at which the adiabatic electronic states become degenerate.^{31–35} These intersections exist as *N _{DOF}* − 2 dimensional seams in nuclear configuration space (where

*N*is the total number of nuclear degrees of freedom). Thus, to a reasonable approximation one need only to consider quantum mechanical transitions between adiabatic electronic states in a finite region surrounding these intersection seams, and adiabatic propagation is sufficient in the remainder of nuclear configuration space.

_{DOF}However, expansion of the total molecular wave function as a sum over vibronic wave functions in the adiabatic electronic basis,

is not without problems. (Here *ψ*_{ad,J}(**r**; **R**) is the adiabatic electronic wave function of state *J*, $XJR,t$ is the associated nuclear wave function, **r** and **R** are the electronic and nuclear coordinates, respectively, and *t* is the time.) Operation of the nuclear kinetic energy operator on the vibronic basis functions summed over in Eq. (1) results in the well-known first-derivative nonadiabatic coupling,

and diagonal Born-Oppenheimer correction (DBOC),

where *η* indexes nuclear degrees of freedom and *m*_{η} is the associated mass. Because the adiabatic electronic wave functions change discontinuously with respect to **R** at conical intersections, their derivatives with respect to **R** are singular at these points; thus, both the first-derivative coupling and the DBOC become singular. The singularities in the first-derivative couplings are integrable, and various strategies for handling them gracefully in numerical simulations of nonadiabatic dynamics are gaining popularity.^{36–41} However, the singularities in the second-derivative couplings are not integrable and therefore more difficult to incorporate in practice.

In fact, singularities in the DBOC are reflections of the discontinuity of the adiabatic vibronic basis functions, *ψ*_{ad,J}X_{J}, in which the molecular wave function is expanded. We have recently shown that the representation of most well-behaved molecular wave functions via the expansion in Eq. (1) requires the adiabatic nuclear wave functions, $XJR,t$, to be discontinuous at the intersection points;^{42} continuity of the total molecular wave function is achieved by cancellation of compensating discontinuities in the nuclear and electronic wave functions of the intersecting states. Upon such compensation and only upon such compensation the singular DBOC contributions cancel, resulting in a finite molecular energy. Managing these discontinuities and singular matrix elements is impractical in the context of most approximate calculations.

This formal analysis of the continuity of the total molecular wave function suggests that more accurate results are obtained when the DBOC is neglected in certain classes of approximate quantum dynamical and mixed quantum classical simulations.^{42} After all, as per the postulates of quantum mechanics, the total molecular wave function must be continuous; the singular DBOC arises from discontinuities in the vibronic basis which must cancel one another. Izmaylov and co-workers have also argued for the neglect of the DBOC in many cases on the grounds that (a) errors due to the neglect of the DBOC are canceled by errors introduced by the neglect of geometric phase, and (b) in many cases simulations of dynamics at conical intersections which include the DBOC show greater errors than those that neglect it.^{43–45} Despite the mounting evidence against inclusion of the DBOC in most practical simulations of nonadiabatic dynamics, it is unsettling to approximate divergent matrix elements as zero.

In order to circumvent the difficulties of the adiabatic representation, many adopt a quasi-diabatic electronic representation.^{46–56} Though a nontrivial completely diabatic basis—i.e., an electronic representation in which the derivative couplings are exactly zero—does not exist for most polyatomic molecules, a quasi-diabatic representation eliminates the singularities at conical intersections, leaving behind a small non-removable nonadiabatic coupling.^{48} Unfortunately, there are infinitely many possible quasi-diabatic representations, leaving the user a difficult choice. The optimal recipe for a well-behaved global quasi-diabatic representation for a particular chemical problem may not be obvious, even when one has significant prior knowledge of the photodynamics of the system in question. Selection of a poorly behaved quasi-diabatic representation may have catastrophic consequences; for example, unexpected energy discontinuities and singularities in the nonadiabatic couplings may be introduced at geometries other than the conical intersection seam.^{57} In addition, quasi-diabatic states are often coupled over a broader swath of nuclear configuration space than adiabatic states, complicating the simulation of electronic transitions.

Thus, it would be ideal to develop an approach to model nonadiabatic dynamics that eliminates the disadvantages of the adiabatic representation—namely the discontinuities in the vibronic basis functions and their consequences—while maintaining its advantages—the relative locality of the coupling between electronic states and the absence of the need to choose a quasi-diabatic basis that is globally well-behaved. With these goals in mind, we present a new quantum dynamical ansatz: diabatized Gaussians on adiabatic surfaces (DGAS). DGAS is a variant of the *ab initio* multiple spawning (AIMS) method for modeling nonadiabatic molecular dynamics^{12,58} that allows us to simultaneously achieve the best features of both electronic representations.

In this work, we will focus on one of the most troubling consequences of the expansion of the molecular wave function in the adiabatic representation: the divergent DBOC matrix elements at conical intersections. To this end we compare simulation results obtained using three variants of AIMS: the newly introduced DGAS approach, traditional AIMS in the adiabatic representation (adiabatic AIMS) in which the DBOC is neglected, and a variant of adiabatic AIMS including the DBOC (adiabatic AIMS + DBOC). In Section II, we will introduce both the DGAS and adiabatic AIMS + DBOC methods in detail. In Section III, we will apply traditional adiabatic AIMS, DGAS, and adiabatic AIMS + DBOC to model the photochemistry of a simple molecule, ethene (ethylene), in order to quantify the degree to which wave function continuity and the DBOC influence the results of approximate simulations of nonadiabatic dynamics. In Section IV, we quantify the errors attributable to the approximations made in the adiabatic AIMS + DBOC approach. In Section V, we draw conclusions.

## II. THEORY AND METHODS

In this section, we introduce two novel variants of the AIMS method that will be employed in this work. In Subsection II A, we review the AIMS method as it is traditionally presented. In Subsections II B and II C, we introduce the DGAS formalism and discuss some important details of its implementation. In Subsection II D, we will introduce a second variant of AIMS: the traditional AIMS method in the adiabatic representation with the addition of the singular DBOC terms that arise from discontinuities in the vibronic basis.

### A. Review of AIMS

AIMS is a hierarchy of approximations that allows the solution of the time-dependent Schrödinger equation (TDSE) for the wave function of an electronically excited molecule, including all electronic and nuclear coordinates explicitly. The AIMS molecular wave function is typically written as an expansion in vibronic states

where $\psi Jr;R$ may be chosen to be either an adiabatic or quasi-diabatic electronic wave function. Most often the adiabatic basis is chosen for the reasons described in the Introduction, giving the expansion in Eq. (1). Hereafter we will refer to this variant of AIMS as adiabatic AIMS. The expansion in Eq. (4) is exact, of course, in both the adiabatic and quasi-diabatic electronic bases; approximations are made in the definitions of $\psi Jr;R$ and $XJR,t$.

The electronic wave functions, $\psi Jr;R$, are computed on the fly via *ab initio* electronic structure calculations. Each nuclear wave function, $XJR,t$, is expanded as a linear combination of Gaussian trajectory basis functions (TBFs),

with

and

Here *N _{J}*(

*t*) is the time-dependent number of TBFs on electronic state

*J*and

*j*indexes TBFs. The user-chosen Gaussian widths, $\alpha j\eta J$, are constant in time (“frozen”

^{59}). The average position and momentum vectors of each TBF—$R\u0304jJ(t)$ and $P\u0304jJ(t)$, respectively—are propagated via the classical equations of motion

and

Here

where $H\u02c6el$ is the electronic Hamiltonian, and the parametric dependence of *ψ _{I}*,

*ψ*, and $H\u02c6el$ on

_{J}**R**is implied. The time-dependent phase, $\gamma \u0304jJ(t)$, is propagated semiclassically

The time-dependent expansion coefficients, $CjJ(t)$, are solved by integration of the TDSE,

Here **S** is the overlap matrix,

**H** is the full molecular Hamiltonian matrix,

and $S\u0307$ is the right-acting time-derivative matrix,

Note that the second equivalences in Equations (13) and (15) arise from the orthonormality of the electronic basis, regardless of whether it is adiabatic or quasi-diabatic.

The Hamiltonian matrix elements can be split into a nuclear kinetic energy term and a nuclear potential energy term that arises from the electronic Hamiltonian

where $T\u02c6R$ is the nuclear kinetic energy operator. Integration over electronic degrees of freedom gives

Because TBFs are non-zero at all **R**, the exact evaluation of matrix elements like this requires knowledge of the electronic wave functions of states *I* and *J* at all **R**. Such infinite knowledge is obviously not possible in practical applications. One can take advantage of the locality of the TBFs, however, to approximate the matrix elements using only local information about the electronic structure. Specifically, matrix elements over the Hamiltonian are typically computed in the zeroth-order saddle-point approximation (SPA),

The centroid position, $R\u0304\u0304ijIJ$, is the maximum of the absolute overlap density, $\chi iI(R,t)*\chi jJ(R,t)$. It is trivially defined as

in the common case that $\alpha iI=\alpha jJ$. As can be seen, computing $Hel,ijIJ$ in the zeroth-order SPA requires knowledge of the value of $VIJR$, and therefore the electronic structure, at only a single point in nuclear configuration space: $R\u0304\u0304ijIJ$. In AIMS, the SPA is employed to approximate integrals over all operators of the form $O(R\u02c6)$ (e.g., nonadiabatic coupling matrix elements (NACMEs) or other electronic observables).

In practice, it is advantageous to minimize the number of TBFs so as to reduce computational cost. With this goal in mind, the initial nuclear basis is typically chosen to be small (often a single TBF), and then is adaptively expanded during the course of the simulation. In an adiabatic AIMS simulation, the NACMEs coupling each TBF to other electronic states are tracked as a function of time during the simulation. When the norm of the NACME vector rises above a user defined threshold, a new TBF is created (*spawned*) on the coupled state to accept population transfer. In this way, TBFs are only added to the nuclear basis set when the algorithm determines that they are required for the accurate solution of the TDSE.

For the sake of brevity, we do not present the specifics of the spawning procedure here, nor do we present simplified expressions for all matrix elements. These and other details of the AIMS algorithm can be found in the literature.^{58,60}

### B. Diabatized Gaussians on adiabatic surfaces

In order to introduce the DGAS approximation, we first generalize the AIMS formalism. Rather than associating each TBF with a set of global adiabatic or quasi-diabatic electronic states, as expressed in Eqs. (4) and (5), we associate a unique, **R**- and *t*-dependent electronic wave function, $\psi jr;R,t$, with each TBF. The AIMS wave function is then defined as

The adiabatic AIMS ansatz described above is a special case of this generalized form of AIMS in which the electronic wave function is the time-independent adiabatic electronic wave function

The recently developed *ab initio* multiple cloning algorithm is also encompassed by this generalization of AIMS.^{61}

The generalized AIMS formalism defined in Eq. (20) is similar to that defined in Eqs. (4) and (5) in most respects. The classical and semiclassical propagation, as defined in Eqs. (8)–(11), is identical. The basis may still be expanded adaptively in the spirit of the spawning procedure described above. The TDSE (Eq. (12)) has the same form, and the matrix elements may still be approximated in the SPA. What differs are the definitions of the overlap and derivative matrix elements in Eqs. (13) and (15); the second equality in these equations does not hold because the electronic wave functions may not form an orthonormal set. Similarly, electronic overlaps that arise in the calculation of other matrix elements in the SPA may take on values intermediate between zero and one. The calculation of the elements of $S\u0307$ is further complicated by the time-dependence of the new electronic wave function

Now we introduce the diabatized Gaussians on adiabatic surfaces (DGAS) ansatz, which takes advantage of the time dependence of the generalized AIMS electronic wave function to eschew the primary difficulties of the adiabatic basis without resorting to a globally quasi-diabatic electronic basis. In DGAS, like in adiabatic AIMS, each TBF is associated with an adiabatic electronic state; thus, we will index each nuclear wave function with both *j* and *J*: $\chi jJ$. However, the electronic wave function, *ψ _{j}*, associated with $\chi jJ$ is not the adiabatic electronic wave function. Instead it is a time-dependent quasi-diabatic electronic wave function

We will call $\psi DGAS,jJr;R,R\u0304jJ(t)$ the *DGAS electronic wave function*.

The DGAS electronic wave function is anchored to the adiabatic electronic wave function at the time-dependent average position, $R\u0304jJ(t)$, according to

At all other **R**, the DGAS electronic wave function does not correspond to a single adiabatic electronic state but instead is defined via a truncated expansion over adiabatic electronic wave functions

where *N _{DGAS}* is the number of states included in the expansion. The coefficients in this expansion, which we term the

*DGAS coefficients*, are defined to maximize the overlap of the electronic wave function at all points

**R**with the adiabatic electronic wave function at $R\u0304jJ(t)$ of the associated TBF,

In short, $\psi DGAS,jJr;R,R\u0304jJ(t)$ is a quasi-diabatic wave function chosen to be equal to the adiabatic electronic wave function at $R\u0304jJ(t)$ and to vary continuously in the surrounding region of nuclear configuration space.

When reasonable simulation parameters are chosen, this DGAS electronic wave function will not exhibit discontinuities in **R** at conical intersections because the discontinuities in the adiabatic electronic wave functions, *ψ*_{ad,I}(**r**; **R**), at conical intersections will be canceled by compensating discontinuities in the **R**-dependent DGAS coefficients. Thus, singularities in the first- and second-derivative coupling matrix elements do not arise when the nuclear kinetic energy operator is applied to the DGAS vibronic basis, even in the presence of a conical intersection. Singular nonadiabatic matrix elements at conical intersections are the primary disadvantage of the adiabatic electronic basis, and their absence is the primary advantage of the DGAS ansatz. (A secondary advantage of the DGAS ansatz is that $\psi DGAS,jJr;R,R\u0304jJ(t)$ naturally carries with it geometric phase, which is practically challenging to include in simulations in the adiabatic representation.)

Though $\psi DGAS,jJr;R,R\u0304jJ(t)$ is continuous in **R** at conical intersections, the time dependence of the DGAS electronic wave function introduces a new complexity. As seen in Eq. (22), the time-derivative of the DGAS electronic wave function contributes to $S\u0307$. The application of the zeroth-order SPA to the term containing this derivative yields

When the centroid of a TBF passes directly over a conical intersection, the adiabatic electronic wave function at $R\u0304jJ(t)$ changes discontinuously in time. Because the DGAS electronic wave function is anchored to this adiabatic electronic wave function at $R\u0304jJ(t)$, it will also change discontinuously in time, and the rightmost bracket in Eq. (27) will become singular.

(This bracket is the *time-derivative coupling*, which is well known to be related to the first-derivative nonadiabatic coupling by the chain rule,

Hereafter we will use the term “first-derivative coupling” to refer to the coupling arising from the operation of the nuclear kinetic energy operator on the electronic wave function (Eq. (2)), while we will use the term “time-derivative coupling” to refer to that which arises from the operation of the time-derivative operator on the time-dependent DGAS electronic wave function.)

Thus, in moving from adiabatic AIMS to DGAS, we have replaced a basis that is discontinuous with respect to **R** with one that may be discontinuous with respect to *t*. This discontinuity in *t* is advantageous relative to a discontinuity in **R** because the TDSE contains derivatives of the wave function up to second order with respect to **R** but only first-order derivatives with respect to *t*. It is second-order terms (i.e., the DBOC) that exhibit non-integrable singularities at conical intersections. Like the first-derivative coupling, the singularities arising from the time-derivative coupling terms are integrable.

With such integrable singularities in the time-derivative coupling in mind, we have recently developed an efficient approach to integrate the TDSE based on a norm-preserving interpolation (NPI) of the electronic wave functions. Integrating the TDSE with NPI ensures the prediction of accurate population transfer probabilities even in the presence of such singularities.^{40,41} In Subsection II C, we describe the adaptation of this approach, which was originally formulated in the adiabatic representation, to the DGAS ansatz.

Before continuing, we briefly note an interesting relationship between DGAS and fewest switches surface hopping (FSSH). No second-derivative coupling matrix elements arise in Tully’s original formulation of FSSH,^{7,10,62} just as in DGAS, and the absence of these matrix elements arises for similar reasons; the derivative couplings in both DGAS and FSSH arise directly from the differentiation of the electronic wave function with respect to time, not with respect to **R**. Though the formal absence of the second-derivative couplings from FSSH is often viewed as an approximation, DGAS demonstrates that such absence can be achieved in a rigorous, fully quantum mechanical ansatz.

### C. Calculation of derivative couplings in DGAS

As noted above, rapid changes in the adiabatic electronic wave function can result in sharp spikes in the time-derivative nonadiabatic coupling (as reflected in the $S\u0307$ matrix) near conical intersections. Such spikes become singularities if the centroid of one of the interacting TBFs passes directly over the intersection itself. The integration of the TDSE in the presence of such spikes is numerically challenging; when implemented naively, errors in the predicted population transfer probabilities may be large if the coupling varies on a time scale shorter than the time step. Such errors can be reduced by the application of an adaptive or multiple-time-step integrator in most cases;^{63,64} however, this comes at a high computational cost, and in the worst case such a short time step is required that the calculation becomes intractable. NPI addresses this problem by interpolating the electronic wave function between time steps. With this interpolated electronic wave function one can compute the time-derivative coupling as a continuous (as opposed to discrete) function of time.

Here we describe the adaptation of the original formulation of NPI to compute the DGAS nonadiabatic coupling term (Eq. (27)). Consider a single simulation time step from time *t*_{0} to *t*_{0} + Δ*t*. The energies and electronic wave functions at the centroids of two TBFs are known at these two times but no intermediate time. Here we consider configuration interaction (CI) wave functions, but this formalism could trivially be generalized to other forms of electronic wave functions. The CI vector of adiabatic state *J* at position **R** will be represented $Xad,JR$. The molecular orbitals at time *t*_{0} + Δ*t* are rotated (diabatized) to maximize the overlaps with those at *t*_{0}, and we then approximate the overlap between orbital *m* at time *t*_{0} and orbital *n* at time *t*_{0} + Δ*t* as *δ _{mn}*. The electronic overlap matrix elements between the adiabatic electronic wave functions of states

*I*and

*J*at times

*t*

_{0}+ Δ

*t*and

*t*

_{0}, respectively, are thus well approximated as the inner products of CI vectors,

The CI vector of the DGAS electronic wave function at the centroid geometry can be defined in terms of the DGAS coefficients,

Now we define the interpolated DGAS electronic wave functions at this centroid geometry as a continuous function of time, *τ*, between *t*_{0} and *t*_{0} + Δ*t*,

where

and

In short, Eqs. (31)–(36) define a linear rotation of the CI vector from its value at time *t*_{0} to that at *t*_{0} + Δ*t*.

With the CI vectors now defined as a continuous function of time, the time-derivative coupling between DGAS electronic wave functions can also be computed as a continuous function of time

The derivative coupling at the midpoint of the time step is approximated as the average of this continuously defined coupling over the time step

In DGAS, it is this average value that is employed in Eq. (27). As with our previous formulation of NPI, the advantage of this approach is that the entire change in the electronic wave function is resolved at every time step. NPI does not incur error when “stepping over” a conical intersection or narrowly avoided crossing in a single time step, so long as the time step is shorter than the timescale for decoherence and dephasing.^{40}

This adaptation of NPI to DGAS shares the simplicity of implementation of the version reported in our previous publications.^{40,41} The differentiation and integration in Eq. (38) are performed analytically leading to a relatively simple set of expressions for the time-derivative coupling. These expressions, which are reported in the supplementary material, require minimal input from electronic structure calculations: CI vectors of the adiabatic electronic states at $R\u0304iI(t0)$, $R\u0304iI(t0+\Delta t)$, $R\u0304jJ(t0)$, $R\u0304jJ(t0+\Delta t)$, $R\u0304\u0304ijIJ(t0)$, and $R\u0304\u0304ijIJ(t0+\Delta t)$.

As in adiabatic AIMS, nonadiabatic coupling terms also arise from the nuclear kinetic energy operator, $T\u02c6R=\u2211\eta 12m\eta \u22022\u2202R\eta 2$, in DGAS. However, due to the quasi-diabatic nature of the DGAS electronic wave functions, the singularities in the DBOC and the first-derivative couplings that exist at conical intersections in adiabatic AIMS are eliminated. What is left behind is the finite non-removable portion of the nonadiabatic couplings. This non-removable portion is small in general^{65} and, as such, we approximate the nonadiabatic coupling terms arising from the Hamiltonian to zero,

and

### D. Incorporation of DBOC into adiabatic AIMS

Separate from DGAS, in this work we introduce a second modification to the AIMS method: the addition of approximate DBOC terms to the PESs employed in adiabatic AIMS. We introduce this method to quantify the unphysical effects of the wave function discontinuities that arise at conical intersections when the total molecular wave function is expanded as a sum over vibronic states in the adiabatic representation. The only change to traditional adiabatic AIMS is that we add the DBOC to the PES of each state, i.e., we replace the Born-Oppenheimer surfaces with the Born-Huang surfaces,

Each term in the summation can be expanded as a sum over all electronic states,

Here *E*_{ad,I} is the energy of the adiabatic electronic state *I*. Despite the facts that the DBOC can be very important in thermochemistry and kinetics,^{66–80} analytic techniques for computing second-derivative coupling terms are known,^{67,69} and several studies of the influence of the DBOC on molecular dynamics have been performed using analytic PESs,^{10,18,44,62,81–83} the incorporation of the DBOC into *ab initio* molecular dynamics simulations (AIMD) has to date been impractical. This is for two reasons: (a) the summation over all electronic states in Eq. (42) is intractable and (b) the analytic differentiation of the Born-Huang PES requires analytic differentiation of this expression, which in turn requires computationally expensive second-derivative terms. Here we introduce approximations to circumvent these difficulties. First, we truncate the summation in Eq. (42) to include only those states that we directly simulate in the dynamics. This truncation is justified by the fact that only those terms in which *I* and *J* intersect become singular; thus, these few terms dominate the DBOC near conical intersections.

To circumvent the need to compute second-derivative terms, we approximate the gradient of the DBOC. Specifically, to approximate this gradient at position, **R**_{0}, we build an analytical *model cone* approximation to the PES in the region surrounding that point. The model cone PESs are the adiabatic surfaces obtained by the diagonalization of a two-state Hamiltonian in which all matrix elements are Taylor expanded to first order in **R**,

where

and

At **R**_{0}, this model Hamiltonian reproduces the exact Born-Oppenheimer energies and energy gradients of the intersecting PESs. The model and exact Hamiltonians also predict exactly the same contribution to the DBOC from a given adiabatic electronic state, *I*,

Here *ψ*_{model,J} and *E*_{model,J} are, respectively, the eigenstates and eigenvalues of **H**_{model}.

However, exactly computing the derivative of the DBOC would require differentiation of the model Hamiltonian with respect to the expansion point, **R**_{0}, which in turn requires the computationally expensive second-derivative terms (i.e., first derivatives of **g**_{ad,I}, **g**_{ad,J}, and **h**_{ad,IJ}). To avoid the calculation of such derivatives, we approximate the derivative of the DBOC of the exact Hamiltonian as the derivative of the model Hamiltonian DBOC at fixed **R**_{0}

Because the model Hamiltonian is trivially differentiable to arbitrary order in **R**, expensive higher derivative terms need not be computed. We quantify the error in this approximation in Section IV.

It was noted above that singularities in the DBOC are non-integrable, and thus one might expect Born-Huang potential energy matrix elements,

to be undefined in the presence of a conical intersection. However, in the SPA, these matrix elements remain finite so as long as $R\u0304\u0304ijJJ$ does not fall exactly on the intersection seam. Thus, in the SPA the divergent DBOC integrals are approximated as finite. This is arguably a harsh approximation, but still less harsh than the standard practice of completely neglecting the DBOC.

## III. THE EFFECT OF ENFORCING CONTINUITY

As discussed in the Introduction, the expansion of the total molecular wave function as a sum over vibronic states in the adiabatic representation is problematic; discontinuities in these vibronic wave functions result in non-integrable singularities in the DBOC. Most often simulations based on such an expansion neglect the divergent DBOC terms (as in adiabatic AIMS), though through careful approximation one can include them in a practical manner (as in the new adiabatic AIMS + DBOC approach introduced in Subsection II D). More ideal is to eliminate the discontinuities altogether, as we do in the DGAS approach (Subsection II B). In this section, we compare DGAS, adiabatic AIMS, and adiabatic AIMS + DBOC by application to a molecule with well-studied photochemistry: ethene.^{12,33,84,85} There are two goals of this comparison. First, in comparing DGAS and adiabatic AIMS + DBOC, we wish to quantify the effect of discontinuities in the vibronic basis and the associated energy spikes on nonadiabatic dynamics. Second, by the comparison of DGAS and adiabatic AIMS (without DBOC), we investigate two different methods of eliminating the singular DBOC energy terms; in DGAS the wave function is carefully constructed to eliminate these singularities, while in adiabatic AIMS the vibronic basis remains discontinuous, but the resulting singular DBOC terms are simply neglected.

In all three sets of simulations, the PES will be computed at the state-averaged complete active space self-consistent field (SA-CASSCF) level of theory. The state average includes three states, the active space contains two electrons in two orbitals (SA3-CAS(2/2)), and the 6-31G** basis is used. In each case, 50 simulations were initiated with a single TBF on the S_{1} (*π* → *π*^{∗}) state. The initial position and momentum vectors of each TBF are drawn at random from the ground state vibrational Wigner distribution computed in the harmonic approximation at the B3LYP/6-31G** level of theory. The same set of 50 initial conditions was employed in all three sets of simulations. The classical and quantum degrees of freedom are integrated using adaptive integrators (velocity Verlet and second-order Magnus expansion, respectively) with a minimum time step of 0.1 a.u. (0.002 fs). Time-derivative couplings were computed using NPI in all three sets of simulations. The DGAS simulations use NPI as defined in Subsection II C, while the two variants of adiabatic AIMS use NPI as presented in our past work.^{41} For computational convenience, each simulation is stopped when 90% of the population reaches the ground state. All electronic structure calculations were performed with the MolPro electronic structure software package.^{86–88} All AIMS calculations (including DGAS) were performed with the FMSMolPro package.^{60}

The average S_{0} populations of each set of simulations as a function of time after excitation are reported in Figure 1. The adiabatic + DBOC simulations predict a considerably slower decay to the ground state than DGAS. Fitting of the data in Figure 1 to a delayed exponential function gives an excited state lifetime of 50 fs for the adiabatic + DBOC calculations compared to 28 fs for DGAS. (Details of the delayed exponential fit are provided in the supplementary material.)

That inclusion of the DBOC slows population transfer can be rationalized intuitively. The DBOC energy is always positive, and as such the spikes in the DBOC at conical intersection points are purely repulsive. Thus, inclusion of the DBOC reduces the probability that a given TBF passes close to the conical intersection seam, where the coupling between adiabatic states is strongest. We emphasize that these spikes are a consequence of the fact that the adiabatic vibronic wave functions exhibit discontinuities. They do not exist in the DGAS approximation, which is designed to provide a well-behaved molecular wave function at conical intersection points.

The analysis of the couplings observed in our dynamics simulations shows that this intuitive picture is born out in practice. Figure 2 presents the DBOC experienced by the 50 initially excited TBFs as a function of time. Spikes in the DBOC above 0.25 eV are regularly seen, with the largest spike nearly reaching 1.5 eV. Clearly significant repulsive forces are attributable to the DBOC.

That this large repulsive potential reduces the probability of a TBF accessing a region of small energy gap (and thus large nonadiabatic coupling) can be seen in Figures 3 and 4. Figure 3 presents the logarithm of the S_{1}-S_{0} energy gap (in eV) as a function of time for each initially excited TBF. The DGAS trajectories are shown in panel (a) while the adiabatic + DBOC trajectories are shown in (b). The averages of the logarithms of these energy gaps are compared in Figure 4. It can be seen that TBFs in the DGAS simulations regularly explore regions of smaller energy gap than in the adiabatic + DBOC simulations and that the average energy gap is smaller in the DGAS simulations than in the adiabatic + DBOC ones.

Now we consider the adiabatic AIMS results without the DBOC. As noted above, adiabatic AIMS is based on the expansion of the total molecular wave function as a sum over vibronic states in the adiabatic representation, and therefore the individual vibronic states exhibit discontinuities at conical intersection points. The DBOC (which is therefore singular at intersection points) is neglected, however. It is noteworthy that the S_{0} population curves for DGAS and adiabatic AIMS (Figure 1) are effectively identical; the rate of decay is similar regardless of whether the elimination of singularities in the DBOC is achieved by careful construction of a continuous molecular wave function (DGAS) or by simple neglect (adiabatic AIMS). This supports the widespread practice of neglect of the DBOC in simulations in the adiabatic basis. Given the classical nature of the TBF basis, this result also supports the neglect (or formal absence) of the DBOC in mixed quantum classical simulations.

## IV. ANALYSIS OF ERRORS IN DBOC GRADIENTS

The adiabatic AIMS + DBOC simulations are based on the approximation to the gradient of the DBOC defined in Eq. (48). Here we analyze the error in these approximate gradients by comparison to numerical gradients of the DBOC. Specifically, we have optimized the minimal energy conical intersection (MECI) of ethene at the CASSCF level of theory described in Sec. III. The difference gradient and nonadiabatic coupling vectors at the intersection point were computed, and a grid of geometries displaced from the MECI geometry in these two directions was constructed. The electronic structure was solved at the CASSCF level of theory at each grid point. The DBOC is then computed at each point with the expansion in Eq. (42) truncated to the two intersecting states. The derivative of the DBOC was then computed analytically according to Eq. (48) and numerically by central difference with a step size of 0.001 Å. This numerical derivative is taken as the reference, and the error in the model cone derivative is computed relative to it.

The DBOC itself is plotted in Figure 5(a). Figure 5(b) shows the relative error in the gradient vector, defined as

where **g**_{DBOC,num} is the numerical derivative of the DBOC with respect to **R**, **g**_{DBOC,model} is the analytical derivative of the DBOC in the model cone approximation (the right-hand side of Eq. (48)), and **g**_{ad,S1} is the gradient of the Born-Oppenheimer PES of the first excited state. In other words, this is the norm of the error in the model cone gradient vector divided by the norm of the total force on S_{1}. Not surprisingly, this relative error is largest in the region where the DBOC itself is large, exceeding 1% of the total gradient norm in roughly the same region where the DBOC exceeds 10^{−0.5} (0.32) eV. The maximum error is 6% (in the red region in the upper right quadrant of Figure 5(b)).

The error in the gradient results in the non-conservation of the classical energy of the 50 initial TBFs in the adiabatic AIMS + DBOC simulations described in Sec. III. The mean deviation (mean absolute deviation) of the energy at the end of the simulations relative to that at the beginning of the simulations is 0.06 (0.17) eV. Though not insignificant, this is a small fraction of the average total kinetic energy at the end of the simulations (2.99 eV); thus, we expect that it does not affect the conclusions of our study. Note that these averages include 49 of the 50 initial TBFs. The remaining TBF experienced a massive gain of energy (36.8 eV) due to an artifact of the CASSCF PES and was therefore excluded as an outlier.

## V. CONCLUSIONS

In this work, we have introduced DGAS, a new approximate quantum dynamical ansatz that simultaneously achieves the advantages of expanding in the adiabatic and quasi-diabatic bases. Like the quasi-diabatic representation, (A) the total molecular wave function is expanded in a basis of vibronic states that is well behaved at conical intersections, and (B) all Hamiltonian matrix elements are finite. However, similar to the adiabatic representation: (C) there is no need to choose a globally quasi-diabatic representation of the electronic wave function, and (D) the couplings between electronic states spike intermittently and are negligibly small the remainder of the time.

For comparison, we introduced a variant of the AIMS method in which the kinetic energy associated with discontinuities in the vibronic basis is explicitly included. This kinetic energy takes the form of infinite spikes in the DBOC at conical intersections. In practice these spikes repel population from the conical intersection seam, slowing population transfer to the ground state significantly in simulations of photoexcited ethene. Elimination of these spikes either by rigorously eliminating wave function discontinuities (DGAS) or by simple neglect (adiabatic AIMS) results in faster decay. Because spikes in the DBOC result from unphysical wave function discontinuities, we assert that the DGAS rate is more correct. That the adiabatic AIMS results agree with DGAS supports the widespread practice of neglecting the DBOC in the approximate simulation of nonadiabatic dynamics in the adiabatic representation.

Finally, we highlight the key idea behind DGAS, which we hope may stimulate developments in other areas: Many quantum dynamical ansatzes based on the adiabatic representation are designed such that the vibronic basis functions in which that wave function is expanded are discontinuous in **R** but vary continuously in time. In contrast, the DGAS vibronic basis functions are chosen to vary continuously in **R**, but may change discontinuously in time. In this way, all singularities in the derivative coupling terms are eliminated from the Hamiltonian. Instead, singular terms that are akin to first-derivative couplings arise on the time-derivative side of the TDSE. Because these singularities are integrable and occur only in off-diagonal matrix elements, they are far easier to cope with in practice than the divergent DBOC matrix elements that arise in the adiabatic representation.

## SUPPLEMENTARY MATERIAL

See the supplementary material for codable expression for the NPI couplings and details of the fitting of the population data in Figure 1.

## Acknowledgments

We gratefully acknowledge support from the National Science Foundation under Grant No. CHE-1565634, and we thank Artur Izmaylov and Todd Martínez for fruitful discussion.