Alexander
White
*ab,
Sergei
Tretiak
abc and
Dmitry
Mozyrsky
a
aTheoretical Division, Los Alamos National Laboratory, Los Alamos, NM, USA. E-mail: alwhite@lanl.gov; serg@lanl.gov; mozyrsky@lanl.gov
bCenter for Nonlinear Studies, USA
cCenter for Integrated Nanotechnologies, USA
First published on 25th April 2016
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.
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.
The use of complex multi-dimensional Gaussian wave-packets (GWP):
(1) |
In the adiabatic limit, the dynamics can be formally described in the framework of quantum mechanical description of the nuclei, Ψ(x,t) = eiH(x)tΨ(x,0) (here and in the following we set ħ = 1 unless stated otherwise),
(2) |
The non-adiabatic dynamics can similarly be obtained from a quantum mechanical description of the nuclei, |Ψ(x,t)〉 = eiĤ(x)t|Ψ(x,0)〉. Now the nuclei's potential energy operator (x) and the wavefunction |Ψ(x,t)〉 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)〉 = N(1)0g(1)(x;x(1)0,p(1)0,(1)0)|1[x(1)0]〉, where |1[x(1)0]〉 an eigenstate of (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)〉. 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−iĤt into slices e−iĤεe−iĤε… 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 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 (∝ε) transfer to |2[x(1)1]〉. After some calculation one gets66
|Ψ(x,ε)〉 = N(1)1g(1)(x)|1[x(1)1]〉 + εN(2)1g(2)(x)|2[x(1)1]〉. | (3) |
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:
ẋ0(1) = p0(1)−1 | (4) |
ṗ(1)0 = −〈1[x(1)0]|∂xV(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[x(1)0]|∂x|1[x(1)0]〉, 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
(5) |
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.
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
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
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
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.
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.
Footnote |
† Electronic supplementary information (ESI) available: A detailed derivation of the CW-NAMD method and a step-by-step description of the algorithm. See DOI: 10.1039/c6sc01319h |
This journal is © The Royal Society of Chemistry 2016 |