J. Phys. Chem. 1996, 100, 7867-7872
7867
Multiple Time Scale Quantum Wavepacket Propagation: Electron-Nuclear Dynamics Seokmin Shin† and Horia Metiu* Department of Chemistry and Physics, UniVersity of California, Santa Barbara, California 93106-9510 ReceiVed: August 25, 1995; In Final Form: NoVember 22, 1995X
We have developed a multiple time scale (MTS) quantum wavepacket propagation algorithm for the electronnuclear dynamics in which the electron and nuclear motions are followed simultaneously. The large mass disparity between the electron and the nuclear degrees of freedom poses a problem when one propagates the total wavefunction using a single time step. The standard split operator propagation method requires a small time step to be used for stable propagation of the fast electronic degrees of freedom, which leads to an unnecessarily large number of propagations for slow nuclear motions. A new factorization of the time evolution operator based on the decomposition of the total Hamiltonian into the nuclear part and the electronic Hamiltonian is proposed. The MTS factorization allows one to choose the spatial grid size and the time step appropriate for each subsystem separately. The MTS algorithm is applied to a charge-transfer model for the evaluation of optical spectra and thermal rate constants. It is shown that the multiple time scale propagation gives correct results with moderate computational time savings.
I. Introduction The study of the dynamics of molecular systems has traditionally been done by following nuclear motions on a potential energy surface obtained a priori based on the BornOppenheimer approximation. The separation of nuclear and electronic degrees of freedom is a natural approximation for adiabatic processes in which the electrons closely follow the nuclei. However, a single adiabatic PES cannot describe processes involving electronic rearrangements. One can overcome that limitation by using multiple, coupled, potential energy surfaces. Recently, attempts have been made to treat such processes by solving coupled equations for the electronic and nuclear degrees of freedom without constructing a PES.1-4 By treating the dynamics of the electronic motions explicitly in oneelectron problems such as absorption spectra,5,6 charge-transfer rates,7 and electron-mediated surface phenomena,8,9 much detailed physical insights can be obtained. The treatment of the nuclear and the electronic dynamics on equal footing has also been done within the framework of the Born-Oppenheimer separation.10-13 The intrinsic problem one has to face when following the nuclear and the electronic motions simultaneously is the inherent time scale difference due to the large mass disparity. If one uses a single time scale for both the nuclear and the electronic motions, a very small step is necessary for the description of the fast movement of the electrons. As a result, such procedures for electron-nuclear dynamics are very time-consuming. The problem of disparate time scales has been studied in classical molecular dynamics simulations. In the classical simulations the separation in time scales occurs due to the disparity in masses, the difference in the range of forces, and the frequency mismatches among various types of motions. Several methods to solve the problem have been proposed. In the simplest case, one can use different time steps to integrate the equations of motion for the fast and the slow degrees of freedom separately.14-16 Recently, Berne and co-workers have developed multiple time scale algorithms for molecular dynamics simulations.17-21 In their reversible reference system † Present address: Department of Chemistry, Seoul National University, Seoul, Korea 151. X Abstract published in AdVance ACS Abstracts, April 1, 1996.
S0022-3654(95)02498-1 CCC: $12.00
propagator algorithm (r-RESPA), integrators for solving the Newton’s equations of motion in the multiple time scale problem are derived from the Trotter expansion of the classical Liouville propagator. The r-RESPA algorithm for classical simulations shows the stability of the integrators compared to the previous methods and has been applied to various multiple time scale problems with promising results.22-25 The RESPA-type algorithms have been applied to ab initio molecular dynamics (AIMD) simulations. Carter and coworkers have proposed nonreversible and reversible multiple time scale AIMD implementations.13,26,27 They used a small time step to propagate the electronic degrees of freedom (described by SCF or GVB coefficients), while a large time step is used in the integration of the classical equations of motion for the nuclear degrees of freedom. Tuckerman and Parrinello have applied multiple time scale methodology for integrating the Car-Parrinello equations.28 Their approach seeks a separation of time scales within the electronic dynamics and exploits the separation to speed up the propagation of the electronic degrees of freedom. For a simple model system with a few degrees of freedom, one can solve the Schro¨dinger equation exactly for the combined electron-nuclear dynamics. Single electron properties such as the absorption cross section and the charge-transfer rate can be obtained by time-dependent propagation of the total wavefunction. The multiple time scale problem naturally arises for the quantum wavepacket propagation of electron-nuclear dynamics. We use a split operator to perform unitary propagation with different time steps for light and heavy degrees of freedom. Splitting of the time evolution operator has been used extensively for the quantum wavepacket propagation.29,30 Usually, the kinetic energy operator and the potential energy operator are split in order to take advantage of the very efficient FFT method for the kinetic energy propagation.31 In the multiple time scale algorithm suggested here, the splitting will be done between the electronic and the nuclear Hamiltonians. In the present study, we will apply the multiple time scale wavepacket propagation to a 1D charge-transfer model. The characteristics of the algorithm will be examined and sample calculations will be done for the evaluations of the optical spectra by time correlation functions and the reaction rate of the system. © 1996 American Chemical Society
7868 J. Phys. Chem., Vol. 100, No. 19, 1996
Shin and Metiu
II. Method To illustrate the method, we consider a system with two degrees of freedom represented by the nuclear coordinate R and the electronic coordinate x. The Hamiltonian of the system is given by
H(x,R) ) K(x,R) + V(x,R)
(1a)
scheme, assuming that calling FFT routines is the timeconsuming part, is proportional to [γNN(Ne log Ne) + Ne(NN log NN)]. On the other hand, if the propagation is done by ordinary propagators to compute γ steps of length δt for all the degrees of freedom at each time step, the calculation time scales like γ[NN(Ne log Ne) + Ne(NN log NN)]. Therefore, the ratio of calculation time for the ordinary method to that of the MTS scheme will be
where
p2 ∂2 p2 ∂2 K(x,R) ) 2m ∂x2 2M ∂R2
r) (1b)
In the propagation method proposed by Feit and Fleck,29,30 the time evolution operator is factorized in the following way:
exp[iH(x,R)t]) {exp[iH(x,R)∆t]}P ) {exp[iV∆t/2] exp[iK∆t] exp[iV∆t/2]}P (2) where ∆t ) t/P. For the multiple time scale propagation, we first decompose the Hamiltonian into the nuclear part and the electronic Hamiltonian such that
H(x,R) ) HN(R) + He(x,R)
(3)
HN(R) ) KN + VN ) -
p2 ∂2 + VN(R) 2M ∂R2
(4a)
He(x,R) ) Ke + Ve ) -
p2 ∂2 + Ve(x,R) 2m ∂x2
(4b)
The short time split propagator is now given by
exp(iH∆t) ) exp(iHe∆t/2) exp(iHN∆t) exp(iHe∆t/2)
(5)
The time step ∆t is appropriate for the nuclear degrees of freedom but it is too long for the electronic ones. For this reason the time propagator for the electronic degrees of freedom is further decomposed by using a smaller time step δt:
exp(iHe∆t/2) ) [exp(iHeδt/2)]γ
(6)
where γ ) ∆t/δt is the ratio of the two time steps. The propagation scheme provided by eqs 5 and 6 can be summarized as follows: one propagates the electronic degrees of freedom (fast motion) for half of the long time step, then the slow nuclear degrees of freedom for a full, long time step, and then one propagates the fast degrees of freedom to the end of the long time step. During the last propagation the nuclear configurations are different from those used in the first one. In the actual implementation of the multiple time scale propagation introduced above, the individual nuclear and electronic propagators are evaluated by the usual split operator method. The reason for the time saving in the MTS algorithm is that we need the FFT call for the nuclear degrees of freedom only once for every long time step. Usually the number of grid points for the nuclear degrees of freedom (NN) is much bigger than that for the electronic degrees of freedom (Ne). The FFT method scales like N log N when N is the number of grid points. The number of time propagations of the electronic degrees of freedom for NN grid points of the nuclear coordinate is γNN for one large time step. The nuclear degrees of freedom needs to be propagated only once for Ne grid points of the electron. Thus the total calculations time for one large time step using the MTS
γ[log Ne + log NN] [γ log Ne + log NN]
(7)
For example, if we have Ne ) 32 and NN ) 256, then r has the value of 2.24 with γ ) 10 and 2.56 with γ ) 100. This is a modest saving, which may not be worthwhile in calculations requiring little computer time. It is, however, welcome in long calculations. We can consider two limiting cases. First, as (log NN/log Ne) f ∞, the ratio approaches the maximum value: r f γ. However, this limit cannot be obtained easily because of the logarithmic scale in the condition. One can obtain the limiting value of the ratio as γ f ∞. The above equation gives the value of r f 1 + log NN/log Ne, which suggest that one can get bigger savings for the cases of large differences in NN and Ne. For this reason, greater savings are expected for complex systems with more degrees of freedom. III. Numerical Results We apply the MTS propagation method outlined above to a 1D charge-transfer model we have used recently.7 The system consists of three positive ions and one electron. Two ions are fixed in position separated by length L, while one ion and the electron can move along the line joining the fixed ions. The variables describing the system are the electron position x and the position R of the moving ion. The Hamiltonian of the system is given by eqs 3 and 4. The nuclear potential VN(R) is the Coulomb interaction of the moving ion with the two fixed ions. For electronic potential, Ve(x,R), we use a pseudopotential appropriate for the interaction of an electron with cations as developed by Blake et al.32 We take the mass of a hydrogen atom for the ions. The other parameters are similar to those in the previous work.7 A. Propagation Error. As a first check of the algorithm, we examine the propagation error for a wavepacket in the bound system. For a relatively small separation between fixed ions (L ) 5 Å), the model system shows a bound potential for both the electron and the mobile nucleus. A grid of 32 points with ∆x ) 0.92 Å for the electron and that of 128 points with ∆R ) 0.06 Å for the nuclear coordinate are used. The initial wavefunction is a Gaussian in R on the ground electronic state. The position of the initial wavepacket is displaced from the equilibrium minimum position of the potential and the packet is not a stationary state. A time step of ∆t ) 0.5 au (0.012 fs) for both the electron and the nuclear propagation is used for the generation of the reference wavefunction, which is essentially a numerically exact result. The reference wavefunction is compared with that generated by MTS propagation. We calculate the overlap error, defined as 1-|〈Ψref(x,R;t)|ΨMTS (x,R;t)〉|, to study the accumulated error in the MTS method. In MTS propagation, the time step for the electron is fixed at δt ) 0.5 au, while the time step for the nucleus is ∆t ) γδt. Figure 1 displays the error as a function of time for different values of the ratio (γ) of the two time steps. A linear dependence of the accumulated error is found, which is expected for a stepwise integrator.33 The graph shows that MTS propagation gives fairly small propagation error for not too large γ (
