Coupled wave-packets for non-adiabatic molecular dynamics: a generalization of Gaussian wave-packet dynamics to multiple potential energy surfaces.

Accurate simulation of the non-adiabatic dynamics of molecules in excited electronic states is key to understanding molecular photo-physical processes. Here we present a novel method, based on a semi-classical approximation, that is as efficient as the commonly used mean field Ehrenfest or ad hoc surface hopping methods and properly accounts for interference and decoherence effects. This novel method is an extension of Heller's thawed Gaussian wave-packet dynamics that includes coupling between potential energy surfaces. By studying several standard test problems we demonstrate that the accuracy of the method can be systematically improved while maintaining high efficiency. The method is suitable for investigating the role of quantum coherence in the non-adiabatic dynamics of many-atom molecules.

Introduction

First principles-based molecular dynamics (MD) is becoming an important tool for understanding properties of complex molecular systems.1–3 Unfortunately, the cost of exact dynamics, by direct calculation of the time-dependent Schrödinger equation (TDSE), scales exponentially with the dimensionality (i.e. number of atoms) of the system.4–10 Thus, for large systems one often approximates that the nuclei of a molecule propagate via classical equations of motion and calculates forces (due to coulombic interaction) via quantum chemistry methods. In doing so one typically relies on the Born–Oppenheimer approximation, where electrons remain in the same electronic quantum state |n(x))〉 with energy, with energy, E((x), that parametrically depends on nuclear coordinates, x = (x1, …, x).11 Thus, the nuclei propagate on a single potential energy surface (PES). For molecules in the ground electronic state, and at low temperatures, this situation often holds due to a sufficiently wide gap between the PES of the ground and excited electronic states. However, for certain nuclear configurations, common when the molecule is in an excited electronic state due to absorption of energy (e.g. a photon), the gaps can become small or even vanish. In these regions, where the nuclear-electronic coupling is the same order as the energy gap, non-adiabatic behavior is expected.12 This creates a superposition of electronic states, with different forces acting on the nuclei.

Since the full TDSE is numerically intractable for high dimensions, approximations for the non-adiabatic molecular dynamics (NAMD) must be made. The simplest approximation is to average over the electronic degree of freedom (DOF), a mean-filed approximation, to determine the force on the nuclei.13,14 This is known as the Ehrenfest approximation. Like any mean-field approximation, it breaks down when there is non-negligible correlation between the dynamical DOF (the nuclear) and the averaged DOF (the electronic), i.e. if the components of the nuclear wavefunction separate depending on which PES they propagate. In an attempt to correct for this problem, while maintaining efficiency and simplicity, Tully proposed the surface hopping method,15 most commonly used with the fewest-switching procedure (FSSH).16 In this method a swarm of classical trajectories propagate on an initial PES, with a finite probability to hop to a coupled PES in regions of non-adiabatic coupling. This method is ad hoc, and is only strictly accurate in the same limit as the Ehrenfest approximation.17,18 This incomplete treatment of the nuclear-electron correlation has two well known symptoms: the interference problem, where the incorrect phase of the nuclear wavefunction leads to incorrect levels of constructive/deconstructive interference, and the decoherence problem, where the separation of the nuclear wavefunction is improperly accounted for. Both problems were pointed out by Tully himself.16 These two methods, Ehrenfest and FSSH, are by far the most commonly used in the simulation of NAMD.19–31

Many attempts have been made to improve upon the basic foundation of these two methods, while retaining the independence of trajectories.17,18,32–40 Quantum coherence effects require first-principles mixed quantum-classical or semi-classical methods, which allow interference between trajectories. These methods are typically applied only in small, or reduced, systems due to inefficiency and/or complexity.41–54 Wave-packet methods, such as ab initio multiple spawning,55 non-adiabatic Herman–Kluk frozen Gaussians,56 and non-adiabatic matching pursuit/split-operator Fourier transform,57 to name a few, provide a more accurate alternative to surface hopping and Ehrenfest methods. Additional complexities and computational costs are involved, but they also benefit from direct calculation of the wavefunction, and thus a clear, unambiguous, route to expectation values.

An ideal NAMD method would have certain properties. It should (1) be based on localized dynamics, i.e. based on real-space trajectories, (2) use only local parameters easily calculated from common electronic structure methods, i.e. PES and electronic wavefunction, (3) require no empirical or ad hoc treatments, (4) include proper treatment of electron-nuclear coupling, (5) be at least as efficient as surface hopping, and (6) be systematically improvable. The coupled wave-packets for non-adiabatic dynamics method, presented in this paper, satisfies each of these conditions. Similar to ab initio multiple spawning, it is built on a framework of branching, due to non-adiabatic coupling, Gaussian wave-packets. These wave-packets form a Gaussian basis for the nuclear wavefunction. The method is systematically improvable, capable of converging to the exact solution, and accurately accounts for decoherence and interference effects. Efficiency is retained, relative to Ehrenfest and FSSH, even when very high accuracy is required.

Extension of thawed Gaussian approximation to multiple PES

To build a new NAMD method, which satisfies the first and second conditions, one must start from a sound foundation for these real-spaced trajectories. The closest analog to a classical particle, and thus real local trajectory, for a quantum system is a localized wave-packet, or superposition of wave-packets.58

The use of complex multi-dimensional Gaussian wave-packets (GWP):as approximations, or basis sets, for nuclei wavefunctions is well studied for semi-classical dynamics on a single potential energy surfaces.56,59–65 In 1975, Heller derived the equations of motion for the four parameters (position x0, momentum p0, complex width matrix α̂0, complex phase γ0) of the GWP, assuming the PES is locally quadratic around x0, the thawed Gaussian approximation (TGA). The key of this method is that the phase-space center of the wave-packet moves by classical mechanics. That classical point is “dressed” in the semi-classical width and phase.59,66

In the adiabatic limit, the dynamics can be formally described in the framework of quantum mechanical description of the nuclei, Ψ(x,t) = eΨ(x,0) (here and in the following we set ħ = 1 unless stated otherwise),where H(x) is the Hamiltonian of the system. The potential V(x) is a parametric function of geometry x, m is the nuclear mass, and Ψ(x,t) is the nuclear wavefunction. TGA can be alternatively derived by splitting the evolution operator operator e– into slices e–e–… with an infinitesimally small time step ε (see ESI†).66 If for a single time slice one expands e to first order in ε, applies the same approximation as Heller, and re-exponentiates, one recovers a new Gaussian with parameters shifted by one time step using Heller's equations of motion.59,66 We seek to follow similar steps to generalize TGA for multiple electronic states.

The non-adiabatic dynamics can similarly be obtained from a quantum mechanical description of the nuclei, |Ψ(x,t))〉 = e = e|Ψ(x,0),0)〉. Now the nuclei's potential energy operator . Now the nuclei's potential energy operator V[combining circumflex](x) and the wavefunction |Ψ(x,t))〉 are are M × M Hermitian matrix and M component vector respectively, where M is the number of relevant electronic states. For simplicity we will consider a situation with two levels crossing, i.e., with M = 2. This is the most common situation, typically more complex problems with multiple PESs and crossings can be modeled as consecutive transitions through well separated regions of coupling between two locally adjacent PESs. Furthermore, an extension to the M > 2 situation is straightforward. We assume that the initial state is a single Gaussian localized on the first PES, |Ψ(x,0),0)〉 = = N(1)0g(1)(x;x(1)0,p(1)0,α̂(1)0)|1[x(1)0]]〉, where |1[, where |1[x(1)0]]〉 an eigenstate of an eigenstate of V[combining circumflex](x) corresponding to the first PES evaluated at x(1)0 and N(1)0 is the real amplitude of the otherwise normalized state |Ψ(x,0),0)〉. Here and in the following the superscripts indicate the electronic state or PES. Non-Gaussian states can be treated as linear superpositions of finite number of Gaussians due to the linearity of the TDSE.. Here and in the following the superscripts indicate the electronic state or PES. Non-Gaussian states can be treated as linear superpositions of finite number of Gaussians due to the linearity of the TDSE.

We again split the evolution operator operator e– into slices e–e–… and now introduce a basis resolution between the first and the second slices (the subscripts here and below indicate the time steps). The point x(1)1 is the location of the classical trajectory to be specified below. This is not to be confused with representation of unity used in the traditional Born–Oppenheimer expansion ).67 The use of the prior basis set, which is locally defined by the center of the wave-packet, as compared to the later, which is defined non-locally, by the variable x, is a subtle but crucial deviation from previously derived path-integral GWP dynamics.68,69 The use of a local eigenfunction is in line with a Gaussian wave-packet or trajectory based approach, since the Gaussian wave-packet is fully defined by localized quantities. Physically, introduction of the basis resolution corresponds to projecting the wave-packet on the new, slightly shifted basis of the eigenstates of V[combining circumflex] at the new average position of the wave-packet at time ε. The new wave-packet will mostly remain in the electronic state |1[x(1)1]]〉 with a small (∝ with a small (∝ε) transfer to |2[x(1)1]]〉. After some calculation one gets. After some calculation one gets66

The change in the wave-packet g(1)(x) in eqn (3) (i.e., after a single time step) is infinitesimal with the same form as the Heller GWP dynamics, leading to equations of motion for the multistate case: ṗ(1)0 = – = −〈1[1[x(1)0]|∂V[combining circumflex](x0)|1[x(1)0]]〉,,

The weight, N(1)1 = N(1)0, is unchanged.

The wave-packet g(2)(x) “hopped” to PES 2. It has the same classical position as the original wave-packet, x(2)1 = x(1)0, but different momentum: p(2)1 is such that p(2)1 – p(1)0 is parallel to the non-adiabatic coupling vector d12(x(1)0) = ) = 〈2[2[x(1)0]|∂|1[x(1)0]]〉, and its absolute value satisfies the energy conservation condition, , and its absolute value satisfies the energy conservation condition, .66,70 This rescaled momentum, and thus the idea of “hopping”, is a direct consequence of the projection onto local electronic basis functions. This can be seen in detail in the full derivation presented in the ESI.† 66 The parameters α(2) and N(2)1 are related to the coefficients of g(1)(x) as where v(1)0α = p(1)0α/m, v[combining macron](1)1 = (v(1)0 + v(2)1)/2 and Δv(1)0 = (v(1)0 – v(2)1)/2. Note that the parameters of the spawned wave-packet, e.g.eqn (5), change discontinuously at the moment of the hop. In practice we only keep the linear term in the expansion of V(x), affecting eqn (4) and (5), since for realistic systems calculation of the quadratic term can be very costly.

At the next time step each of the wave-packets propagates and spawns again, according to eqn (3)–(5) (with a replacement 1 → 2 for the wave-packet on the second PES). After each time step the total number of the wave-packets doubles. Such process can be viewed as branching on a tree, shown in Fig. 1a. This branching tree can be evaluated by a Monte-Carlo approach48,50,71 which becomes too computationally expensive in systems with multiple level crossings.

Fig. 1

(a) Branching tree solution to time dependent Schrödinger equation (sampled by Monte-Carlo). (b) Coupled wave-packets for non-adiabatic molecular dynamics (CW-NAMD) approximation to the branching tree. (c) CW-NAMD approximation with coarse branching.

An alternative way to represent the branching tree is to group terms with the same number of hops. In the continuous limit the resulting wavefunction can be written as

Here the first term in the right hand side is the wave-packet that did not hop, while the single integral term, is the sum over the wave-packets that hopped only once at time t1, etc. Note that all the integrals are time-ordered, which insures the convergence of the series for finite t. A full derivation is available in the ESI.† 66

Coupled wave-packets

Here we propose a new approach based on the wave-packet reconstruction after each spawning event. The approach is schematically shown in Fig. 1b. That is, after two time steps, described in eqn (3)–(5), one creates two wave-packets, on each PES, which will give rise to four more, etc. We note, however, that if each pair of the wave-packets on the same surface has close coordinates and momenta, one can replace each pair by a single GWP, with slightly shifted parameters. We parameterize the new Gaussian by calculating the expectation values of x[combining circumflex], p[combining circumflex], x[combining circumflex]2, p[combining circumflex]2 of the superposition. of the superposition. 〈x[combining circumflex]〉 and 〈 and 〉 and 〈p[combining circumflex]〉 are taken as the position and momentum of the reconstructed wave-packet, while 〈 are taken as the position and momentum of the reconstructed wave-packet, while 〉 are taken as the position and momentum of the reconstructed wave-packet, while 〈x[combining circumflex]2〉 and 〈 and 〉 and 〈p[combining circumflex]2〉 directly give the new complex width. The new phase and weight, directly give the new complex width. The new phase and weight, γ and N, are determined by maximizing the overlap of the new wave-packet with the superposition, under the constraint that N2 is the same as the density of the superposition.66 Approximations are made in order to decouple the calculation of Approximations are made in order to decouple the calculation of 〈x[combining circumflex]〉 and 〈 and 〉 and 〈p[combining circumflex]〉 from the explicit form of the wavefunction, from the explicit form of the wavefunction, i.e. α̂.66 Thus, as with Heller's equations, the trajectories of the GWPs remains independent of the phase and width. At the next step the procedure is repeated, again we have only two GWP and so on. The process repeats until the overlap, O12, between the Gaussians within each pair becomes intolerable, O12 < Omin. At this point, or if the non-adiabatic coupling drops below its own threshold, the “coupling” between the GWPs stops and each is treated independently, thus new branching is allowed. This coarse branching, schematically shown in Fig. 1c, significantly reduces or eliminates the exponential growth of the number of wave-packets. We call this approximation coupled wave-packets for non-adiabatic molecular dynamics (CW-NAMD).

As the two wave-packets separate in position space, their electronic bases will become non-orthogonal. Formally this must be taken into account by considering the required basis rotations when reconstruction occurs. These rotations lead to a correction, but it is small and does not affect the results presented in this article.66

Results

Fig. 2 shows scattering probabilities for the standard Tully test problems II and III, Fig. 2d and e. These problems are frequently used to test new methods of non-adiabatic dynamics because they specifically probe the interference (Tully II) and decoherence (Tully III) questions directly. We compare the CW-NAMD results, with Omin = 0, to the standard fewest switching surface hopping (FSSH) and the mean-field Ehrenfest method as well as direct calculation of the time-dependent Schrödinger equation.16 When branching does not occur, the computational cost of the CW-NAMD method is similar to Ehrenfest (i.e. there is one force calculation per surface per time point), and is much lower than surface hopping. Fig. 2a and b demonstrates that for sufficiently high momentum the CW-NAMD method produces the correct scattering results. The CW-NAMD does not suffer from the interference or decoherence errors of Ehrenfest or FSSH. This can be observed by comparing the position of the peaks of the Stueckelberg oscillations72 in Fig. 2a and the lack of false oscillations in the reflected probabilities in Fig. 2b. Unlike Ehrenfest, CW-NAMD produces the correct momenta and positions of the wave-packets on the upper and lower surface (see Fig. 2f). However at low momenta, the total scattering probability is not conserved and may be poorly estimated. This is evident in both Tully-II (see Fig. 2a, 3b) and Tully-I (Fig. 3a). This can be corrected by allowing the coupled GWPs to branch, i.e. set Omin > 0.

Fig. 2

(a and b) Scattering probabilities Tully II (III) problems on the lower (upper) surface for different initial wave vectors k. Exact solution, FSSH (2000 trajectories), Ehrenfest and CW-NAMD are compared. Initial wave-packet position xinitial = –10 a.u. Initial width, αinitial = ik2/400 for all. (c–e) Potential energy surfaces (E) and Non-Adiabatic Coupling Vectors (NACV) for Tully I (II, III). (f) Average momentum for each surface after scattering (Tully I). Exact solution, FSSH, Ehrenfest and CW-NAMD are compared.

Fig. 3

(a and b) Comparison of low momentum transmission probabilities, on lower surface, with different values of Omin, compared to exact solution, for Tully I (II). Initial conditions as in Fig. 2 (except xinitial = –5 a.u. for Tully I). (c and d) Number of “effective” trajectories for Tully I (II) calculation with different values of Omin. Dynamics are run for a total time of a.u. for Tully I (II).

We compare the low momentum results for Tully I (II) with different values of Omin in Fig. 3a and b. The difference between exact and CW-NAMD solutions is systematically improved by increasing Omin. The increased cost can be seen in Fig. 3c and d. In direct dynamics simulations the bottleneck is typically the calculation of the PES gradients (forces). Trajectory methods like FSSH, require one force calculation per time step per trajectory. Thus we define an “effective” number of trajectories, by determining the total number of force calculations (summed over all branches) divided by the total number of time steps for the simulation, to compare the cost of a branching scheme to that of a trajectory based methods (i.e. FSSH). We see a growth of the number of trajectories required with increased Omin, however the cost of CW-NAMD is still lower than the 2000 trajectories used to calculate the FSSH result (Fig. 2b). In the limit Omin = 1 we recover the full branching tree (Fig. 1a). Lower values of Omin result in a coarse-grained tree (Fig. 1c). To prevent overgrowth of the tree, we place hard-limits on the spawning rate and utilize pruning procedures to discard irrelevant branches.66

To further test the capabilities of this new algorithm we apply the method to the three level, three crossing Model X problem (Fig. 4a–c).18 As is common in realistic systems there are multiple interfering pathways, reflection on multiple surfaces, and sharply changing potentials. The wave-packet is initialized on the middle surface, with an initial width, αinitial = ik2/1600, which is four times smaller (broader wave-packet in position space) than for the Tully test suite. There are four district regions of the momentum dependence: (1) when energy is too low for a classical particle to pass the peak in the middle PES (k < 11), (2) when the particle can pass but does not have energy to hop to the upper surface (k < 12), (3) when it can hop but cannot escape the well (k < 16), and (4) when it can transmit on the upper surface. The CW-NAMD method accurately describes scattering probabilities for all regions. The region between k = 12 and k = 16 is especially difficult to simulate for trajectories based methods, due to the many oscillations inside the upper well, and seems numerically infeasible for Monte-Carlo methods.50,71 FSSH results are also shown. The FSSH method shows breakdown at lower momentum due to over-coherence and incorrect phase interference.18

Fig. 4

(a) Model X PES's with non-adiabatic coupling vectors. (b) Transmission probabilities: lower, middle, and upper surface (top to bottom) as a function of initial momentum. (c) Reflection probabilities: lower, middle (top and bottom) as a function of initial momentum. Exact solution, CW-NAMD (Omin = 0.999), and FSSH (10 000 trajectories) results are compared. αinitial = ik2/1600 a.u., xinitial = –12 a.u. (d) Non-separable 2D Well PES's. (e) 2D Well non-adiabatic coupling vectors, x and y directions. (f) Transmission probabilities lower surface, exact solution, CW-NAMD (Omin = 0.99) and FSSH (10 000 trajectories) results are compared.

Finally, accuracy in non-separable multidimensional PESs is key for application to realistic molecular systems. Such a model, introduced by Shenvi et al.,73 involves two coupled nuclear degrees of freedom in a non-trivial geometry (Fig. 4d) with significant non-adiabatic coupling in both directions (Fig. 4e). As in Tully II, strong Stueckelberg oscillations are present in the scattering probabilities. Fig. 4f shows the probability of transmission on the lower surface as a function of initial momentum in the x direction (k), there is no initial momentum in the y direction. Even in the free-thawed Gaussian approximation the oscillations are qualitatively more accurate with CW-NAMD than with FSSH. At very high momentums the FSSH result converges to the CW-NAMD result. In this regime the difference in integrated forces on the two PES are negligible compared to the initial momentum. Quantitative deviation from the exact result is due to the break down the Thawed Gaussian Approximation, the width of the wave-packet being of similar size to the PES well. This is further discussed in the following section.

Discussion and Conclusion

CW-NAMD is based on the thawed Gaussian approximation, which is known to be accurate only for short times, unless the PES is harmonic or linear. In practice we have approximated even further, by only keeping the linear expansion in V(x). The linear expansion, and indeed the full TGA is appropriate when imag[α̂0] ≫ ∂2V[combining circumflex](x0), essentially when the wave-packet is narrower in position than the variance of the potential. When this does not hold error will accumulate over time, thus limiting the procedure to short-time dynamics. However, TGA can be systematically improved through time-slicing procedures, which allow for the broadened wave-packets to be reconstructed from narrower, in position space, wave-packets.74,75 This can be done by Monte-Carlo sampling or gradient-descent,57,76 and adds an additional source of branching which interfaces well with the CW algorithm. By reducing the over-complete basis, the number of trajectories can be kept to the minimum necessary for an accurate description.44,57,75,76 Inclusion of this time-slicing procedure is a goal for the continued development of the CW-NAMD method.

The CW-NAMD method is similar in spirit to the ab initio multiple spawning (AIMS) method developed by Martinez et al.55,77,78 Both involve approximate solution to an infinitely branching tree, of GWPs. However, in practice AIMS is usually based on the so called independent first generation, where an initial sampling of non-interfering frozen Gaussian wave-packets are propagated. Note that subsequently spawned packets do interfere and full interference can be considered as well.55 CW-NAMD uses thawed GWPs, which can be expanded to improve accuracy, and considers the full superposition of GWPs. For AIMS the “spawning” procedure is based on well-reasoned but empirical considerations.79 The branching procedure in CW-NAMD has a simple numerical control parameter, Omin. In AIMS, coupling between spawned wave-packets results in an equation of motion for complex pre-factors (weight and phase), but, unlike CW-NAMD, does not result in shifts of the other Gaussian parameters.

In conclusion, the new CW-NAMD method is a highly efficient and accurate method of simulating non-adiabatic dynamics applicable to realistic molecular systems. CW-NAMD consistently accounts for decoherence and interference between different dynamical pathways. It can be as efficient as the Ehrenfest method in the high momentum limit, moreover it accurately describes the dynamics of branching wave-packets. In the low momentum limit the method can be systematically improved by increase the rate of allowed branching via the user controlled accuracy threshold, Omin. Combined with filtering of insignificant branches, the method is more accurate and more efficient than the standard FSSH. In our test problems we observe numerical cost of CW-NAMD ranging from about 2(M) to 1000 trajectories depending on initial momentum and desired accuracy. This needs to be compared with the number of effective trajectories in other methods: 2(M) (Ehrenfest), (2 – 10) × 103 (FSSH), (2 – 10) × 104 (Monte-Carlo approaches).

The development of CW-NAMD opens new avenues for future research. On the technical development side, this includes more advance branching criterion, manipulation of the electronic bases, optimization of the reconstruction and branch pruning procedures, as well as the introduction of time-slicing. In regards to application, the CW-NAMD opens a path towards investigation of the role of quantum coherent effects in the non-adiabatic dynamics of large scale molecular systems.

Click here for additional data file.