Multidimensional Quantum Trajectory Dynamics in Imaginary Time

J. Phys. Chem. C 2010, 114, 20595–20602

20595

Multidimensional Quantum Trajectory Dynamics in Imaginary Time with Approximate Quantum Potential† Sophya Garashchuk* and Tijo Vazhappilly Department of Chemistry and Biochemistry, UniVersity of South Carolina, Columbia, South Carolina 29208 ReceiVed: June 1, 2010; ReVised Manuscript ReceiVed: August 25, 2010

The quantum trajectory dynamics in imaginary time with the momentum-dependent quantum potential [J. Chem. Phys. 2010, 132, 014112], which can be used for calculations of low-lying energy eigenstates or for Boltzmann evolution, is extended to many dimensions. Evolution of a nodeless initial wave function represented in exponential form is described in terms of a trajectory ensemble. In the Lagrangian frame of reference the trajectories evolve according to the classical equations of motion under the influence of the inverted classical potential and the quantum potential defined by the gradient of the trajectory momenta. The initial wave function is sampled by the trajectories randomly, which makes the approach scalable to high dimensions. Another practical feature is approximate computation of the quantum potential and force from the global least squares fit of the trajectory momenta in a small basis. The importance sampling and the Eulerian frame-of-reference formulation, compensating for divergent dynamics in bound potentials, are introduced to reduce the number of trajectories in multidimensional applications. The zero-point energies are computed for model and chemical (H2O and SO2) systems. 1. Introduction The time-dependent Schro¨dinger equation (SE) is equivalent to the diffusion equation with the quantum-mechanical (QM) Hamiltonian under the transformation of the real-time variable t f -ıτ

∂ψ(x b, τ) ˆ ψ(x H b, τ) ) -p ∂τ

(1)

QM evolution in imaginary time for τ g 0 according to eq 1 has been used to determine low-lying energy levels and eigenfunctions1 including high-dimensional systems using Monte Carlo methods.2-4 Shifting the potential so that all eigen-energies are positive and using the eigenstate expansion of a wave function

ψ(x b, τ) )

∑ cje-E τψj(xb), j

ˆ ψj(x H b) ) Ejψj(x b)

(2)

j

it is easy to see that a wave function, governed by eq 1, decays to the lowest energy state

ˆ |ψ〉 〈ψ|H ) E0 τf∞ 〈ψ|ψ〉 lim

(3)

QM calculations of the zero-point energy (ZPE) and tunneling splitting in malonaldehyde5 and diffusion Monte Carlo study of CH5+6 are excellent illustrations of the imaginary-time propagation concept. The imaginary-time evolution with trajectories introduced by Miller,7 unfolding on the inverted classical potential, is also employed in studies of thermal †

Part of the “Mark A. Ratner Festschrift”.

reaction rates and other properties of complex systems in the Boltzmann (thermal) distribution.8,9 Due to exponential scaling of the conventional basis-type QM evolution techniques,10 approximate or semiclassical trajectorybased methods are of interest (for example see11). Several approximate approaches12-17 are based on the de Broglie-Bohm formulation of the time-dependent SE.18 The polar representation of a wave function in terms of real amplitude and phase, ψ(x b, t) ) A(x b, t) exp(ıS(x b, t)/p), leads to the quantum trajectory dynamics, whose signature is that it follows evolution of the highdensity regions of a wave function. The trajectory ensemble is governed by the classical and quantum potentials. A global approximation of the quantum potential

Ureal ) -

b, t) p2 ∇2A(x 2m A(x b, t)

(4)

makes this framework practical.17 Extensions of the quantum trajectory concepts to the diffusion eq 119,20 focused on the “independent” quantum-trajectory implementations - extraction of the ZPE from evolution of S(x b, τ) and its high-order derivatives for a stationary or nearly stationary trajectory. The imaginary-time evolution of the quantum trajectory ensemble has been given in one dimension in the Lagrangian frame of reference and used for eigenstate and thermal reaction rate calculations.21,22 Below we present the multidimensional formulation in the Lagrangian and Eulerian frames of reference suitable for high-dimensional calculations of the ZPE and corresponding wave functions using Monte Carlo sampling of initial wave functions. Section II describes the multidimensional Lagrangian quantum trajectory formulation in imaginary time and, essential for high-dimensional implementation, approximate determination of the quantum potential using the least squares fit to the trajectory momentum. Numerical examples, including ZPE calculation of triatomic molecules and multidimensional model systems, importance sampling and Eulerian formulation

10.1021/jp1050244  2010 American Chemical Society Published on Web 09/10/2010

20596

J. Phys. Chem. C, Vol. 114, No. 48, 2010

Garashchuk and Vazhappilly

that reduce the number of trajectories are presented in section III. Section IV concludes. II. Multidimensional Formulation A. Equations of Motion in the Lagrangian Frame-ofReference. Let us consider the diffusion eq 1 in N -dimensional space with the QM Hamiltonian defined as 2 ˆ ) - p ∇ · M-1 · ∇ + V H 2

U is a time-dependent constant: contrary to the real-time dynamics, the quantum force in this case is zero and the trajectory motion is purely classical. (The quantum potential U still contributes to the action functions and defines the ZPE.) Transformation of eq 7 into the Lagrangian frame of reference gives the evolution of the imaginary-time action function

dS 1 p +V+U ) b p · M-1 · b dτ 2

(5)

(12)

The energy of a trajectory can be defined as M is a symmetric matrix of particle masses, possibly with off-diagonal elements and ∇ is the vector of spatial derivatives, ∇i ) ∂/∂xi, i ) 1, ..., N. A vector multiplying a matrix from the left is the vector-column transposed. Generalization to arbitrary coordinate systems is similar to the Bohmian formulation of the real-time SE.23 Substitution of the exponential wave function, unlike the real-time quantum trajectory formulation there is no amplitude prefactor

ψ(x b, τ) ) e-S(xb,τ)/p

(6)

ε)

(7)

Using τ as a time variable and defining the trajectory momentum as in real-time Bohmian dynamics

b p (x, t) ) ∇S(x, t)

(8)

in the Lagrangian frame of reference

d ∂ ) +b p · M-1 · ∇ dτ ∂τ

(9)

dx b )b p · M-1, dτ

dp b ) ∇(V + U) dτ

(10)

Compared to the real-time equations of motion of the quantum trajectories,24 in imaginary time dynamics unfolds on the “inverted” classical potential, -V, in the presence of the momentum-dependent quantum potential (MDQP)

U)

p p p ∇ · M-1 · ∇S ) ∇ · M-1 · b 2 2

(∫ ∇·M

dΩτ ) exp

τ

0

-1

)

·b p dt dΩ

(11)

∇U is the quantum force acting on the trajectories in addition to the classical force ∇V. The quantum potential U is a nonlocal quantity dependent on the ensemble of trajectories representing a wave function. U is responsible for all QM effects in the system and being proportional to pM-1 vanishes in the classical limit. The “classical” evolution can be defined by setting U ) 0. Moreover, as noted in ref 25 for a Gaussian wave function,

(14)

Therefore, in calculations of expectation values the explicit contribution of the quantum potential to the wave function compensates the change in the volume element

ˆ 〉τ ) 〈O

∫ O(xb)e-2S(xb,τ) dΩ ) ∫ O(xb(τ))e-2S(xb(τ)) ˜ dΩτ ) ∫ O(x b(τ))e-2S(xb(τ)) dΩ

(15)

b x(τ) is a trajectory position at time τ. S˜ is the “classical” action function, computed ignoring the quantum potential, along the quantum trajectories moving according to the combination of the classical and quantum forces

dS˜ 1 p +V ) b p · M-1 · b dτ 2

Equation 7 leads to the following equations of motion,

(13)

Once a square-integrable initial wave function is represented as an ensemble of trajectories, their evolution according to eqs 10 and 12 leads to several notable properties. (1) Let us define the volume element associate with the initial position of a trajectory, dΩ ) dx1dx2...dxN. In the Lagrangian formulation its time-dependence along this trajectory is

into eq 1 with the Hamiltonian of eq 5 and division by ψ(x b, τ) gives

∂S(x b, τ) 1 b, τ) + V + ) - ∇S(x b, τ) · M-1 · ∇S(x ∂τ 2 p b, τ) ∇ · M-1 · ∇S(x 2

ˆψ H ∂S 1 p +V+U ) )- b p · M-1 · b ψ ∂τ 2

(16)

In a discrete implementation contribution of the kth trajectory to an average quantity is proportional to a time-independent “weight” defined by the initial (possibly) nonuniform sampling (k) w(k) ) exp(-2S˜(k) 0 ) dΩ

(17)

Then the average quantity is simply

ˆ 〉τ ) 〈O

∑ O(bx (k)τ ) exp(-2S˜(k)τ )w(k)

(18)

k

Superscript k labels the trajectory-specific quantities. (2) The average momentum is zero at all times. Using integration by parts for the ith spatial component of the trajectory momentum b p one obtains

Multidimensional Quantum Trajectory Dynamics

〈pi〉 )

J. Phys. Chem. C, Vol. 114, No. 48, 2010 20597

b) -S e dΩ ∫ e-S ∂S(x ∂xi )-

1 2

∫ e-2S|x )-∞(dxi)-1 dΩ ) 0 xi)∞

(19)

i

The change of the wave function norm, Nτ )〈ψ|ψ〉τ, is proportional to the energy of the trajectory ensemble as follows from eqs 7, 13, and 18

dNτ ) -2 dτ

∫ ∂S(xdτb, τ) e-2S(xb,τ) dΩ ) -2 ∫ ε(x, τ)e-2S(xb,τ) dΩ ) -2〈ε〉τ

(20)

(3) Single-valuedness of the wave function given by eq 6 implies that trajectories do not cross. For an initially real positive wave function, as the wave function decays into the ground state of the system according to eqs 1 and 2 the energy of each trajectory approaches the ZPE, E0

ˆ ψ ) limε(x, τ) ) E0 limψ-1H

τf∞

τf∞

(21)

The latter property is important for the ZPE calculations with small number of trajectories. B. Global Approximate Implementation. Our goal is a cheap implementation of the trajectory formalism which can give estimates of the quantum effects on nuclear dynamics in high-dimensional systems. Therefore, we compute the quantum potential and force in eqs 10 and 12 from the global fit to the trajectory momenta. Of course, more accurate semilocal, such as approximation on subspaces,26 or local, such as trajectorycentered least squares fit,27 implementations can also be used. During the numerical trajectory propagation, at each time step we approximate b p(x b, τ) performing the least squares fit28 with b, τ) as a weighting function and using eq 18 to compute ψ2(x the necessary average quantities. If the ith component of the b, τ) is approximated in a small basis bf(x b) of momentum pi(x b, τ) ≈ bf(x b) · b bi(τ), then the optimal values of the size Nb as pi(x bi)2〉 for i ) 1, ..., N, are the coefficients b bi, minimizing 〈(pi - bf · b

B ) S-1P

(22)

b1, b b2, ..., b bN) are the coefficients Columns of the matrix B ) (b of the Least Squares Fit to the components of the momentum b1, b b), p2(x b), ..., pN(x b)). The matrix P ) (P P2, ..., vector b p ) (p1(x b p with the basis functions, Pij ) 〈pi|fj〉. PN) contains averages of b S is a symmetric matrix of the time-dependent basis function overlaps, Sij ) 〈fi|fj〉. Placing the minimum of V at b x ) 0 the fitting basis bf can be taken as the Taylor basis

bf ) (1, x1, x2, ..., xN, x12, x1x2, ...)

The fitting procedure is practical in many dimensions provided the basis size is small, one-two functions per dimension. The linear basis, bf )(1, x1, x2, ..., xN), of the size Nb ) N + 1 is exact for Gaussian wave functions, and can give useful estimates for systems with small anharmonicities often encountered in nuclear dynamics or, at least, indicate how anharmonic the given system is. The quadratic basis bf(x b) is the smallest basis that leads to a nonzero quantum force affecting the dynamics. Its scaling is Nb ) N + 1 + N(N + 1)/2. More fitting functions can be introduced into selected degrees of freedom, known to be substantially anharmonic. As shown for several anharmonic one-dimensional systems for Nb e 6,21 the imaginary-time quantum trajectory approach yields more accurate ground and excited energy eigenvalues as the basis size increases. However, the global high-order polynomial fitting will be impractical in high dimensions, because of the basis size scaling and because of the convergence and stability of the trajectories: high-order polynomial representation of b p will lead to very large quantum forces on the fringes of the MDQP wave function. Semilocal lower order approximations, defined on spatial cells or domains, might be a better choice if very high accuracy is required. Here we consider only small, linear through cubic, global fitting bases that are cheap in high dimensions. C. Gaussian Wavepacket in a Well. To visualize the imaginary-time quantum trajectories describing an eigenstate let us consider evolution of a Gaussian wave function in the harmonic well V, characterized by the “coherent wavepacket” width parameter, Rc

Rc )

mw , 2

V ) V0 +

mw2x2 2

(24)

The time-dependent Gaussian wave function

ψ(x, τ) ) exp(-Rτ(x - qτ)2 - nτ)

(25)

is described by three parameters that are functions of time: qτ, position of the wavepacket center, Rτ, Gaussian “width” characterizing the wave function localization and nt describing the wave function norm

R0 cosh wτ + Rc sinh wτ R0 sinh wτ + Rc cosh wτ

(26)

q0 cosh wτ + (Rc /R0) sinh wτ

(27)

Rτ ) Rc

qτ )

For the coherent initial wavepacket, R0 ) Rc, eqs 26 and 27 are simplified

Rτ ) Rc

and

qτ ) q0 exp(-wτ)

(28)

(23) Evolution of the norm parameter is found from

This basis can be also used to represent wave function with nodes, such as excited states, as described in ref 21. The main cost (in addition to solving Newton’s equations of motion) is the solution of the system of linear equations SB ) P, performed once per time-step and yielding results for all dimensions and all trajectories.

dnτ pRτ mw2q2τ ) + dτ m 2 which for R0 ) Rc becomes

(29)

20598


Rcq02 pw nτ ) n0 + + V0 τ + (1 - e-2wτ) 2 2

(

)


(30)

The wave function decay is dominated by exp((-V0 - pw/ 2)τ), as the last term in eq 30 rapidly becomes a constant. Unlike the real-time representation of an eigenstate described in terms of stationary quantum trajectories, the imaginary-time quantum trajectories evolve in time and their divergence describes the decay of the wave function norm. For a coherent Gaussian, ∇p ) mw, the trajectories diverge exponentially with time and the volume element changes as dxτ ) exp (wτ) dx0 as shown on Figure 1a. The panel b gives an example of a more complicated trajectory dynamics corresponding to a double-peaked ground state of the double well potential, V ) (x2 - 1)2/4, for a particle of mass m ) 20. The trajectories in the central region converge (without crossing), whereas the outer trajectories exhibit divergent behavior similar to the harmonic well. Time-dependence of the Gaussian wavepacket parameters and normalized energy, E/N, converging to the ZPE values are shown on Figure 2a. As clear from eq 30 the wave function decay can be adjusted by changing the constant V0. Snapshots of the trajectory energies defined by eq 13 for several instances of time are given on Figure 2b. The energies ε(k) of all trajectories approach the ZPE, E0, at longer times, a feature behind using “independent” quantum trajectories in ZPE calculations.19,20 Overall, the Lagrangian quantum-trajectory formulation for real and imaginary-time SE is well-suited to calculations of the thermal reaction rates, where real-time Hamiltonian dynamics follows evolution with the QM Boltzmann operator.22 However, for the ZPE calculations (or low temperature Boltzmann evolution of bound degrees of freedom) divergent trajectory dynamics leads to problems with sampling of the wave function in many dimensions. If the initial distribution of the trajectories is uniform (or even Gaussian uniform28), most trajectories evolving in the inverted bound potential according to eqs 10 leave the area of the high wave function density; the trajectories also acquire large action functions, S, making their contribution to the expectation values negligible. This representation problem is mitigated by the importance sampling in section III.A.

Figure 1. Exact imaginary-time quantum trajectories describing evolution of a Gaussian wave function (a) in the quadratic well and (b) in the double well potentials. The trajectory positions as functions of time τ are shown with solid lines. The dashes outline the ground state wave functions; their amplitudes are shown in arbitrary units along the horizontal axes.

III. Multidimensional Implementation and Examples A. Importance Sampling. In the ZPE calculations it is reasonable to start with a Gaussian wave function

ψ(x, 0) )

( ) 2R0 π

1/4

exp(-R0(x - q0)2)

(31)

for every dimension. Ideally, its parameters, R0 and q0, describe the ground state of the quadratic approximation to the potential, whose ZPE can be used to estimate the shift V0 (section II.C), reducing the decay of the wave function. The propagation time necessary to obtain ZPE will generally depend on how close the initial wave function is to the ground state of the full potential. To make the trajectory representation of a wave function more efficient by compensating their divergent behavior, one can use importance sampling. As illustrated on Figure 1 at long times only the central trajectories will have nonnegligible contributions to the wave function. Given that the least squares fit procedure of section II.B is rigorous for finite sums as well as for the integrals we can define Gaussian importance sampling to initialize more trajectories near q0. These trajectories will not move far from the center of the well during

Figure 2. Gaussian wavepacket evolution in the quadratic potential. (a) The norm N, normalized energy, the wavepacket center and dispersion are shown as functions of time. (b) The energy of trajectories as a function of position are shown for τ ) 0, 1, 2, and 3; ε(x) approaches the ZPE value at long τ.

propagation. The parameter for the importance sampling is estimated from the time-dependence of the Gaussian width parameter given by eq 28. Introducing the importance sampling function

g(x) ) (2γ/π)1/4 exp(-γ(x - q0)2),

γ ) R0 exp(R0β) (32)

g2(x) will define the uniform Gaussian sampling of ψ2(x, 0) given by eq 31. Then, the initial values of S are S(x, 0) ) (R0 - γ)(x


Figure 3. Normalized energy of two Morse oscillators computed from the MDQP dynamics of 5000 trajectories with uniform Gaussian sampling (β ) 0) and using importance sampling defined by eq 32 for the parameter values equal to the final propagation time, β ) 0.2, and twice that value, β ) 0.4.

Figure 4. Ground state energy of water obtained with the imaginary time propagation of quantum trajectories. The results are obtained from propagation of 300 000 quantum trajectories with MDQP using linear, Nb ) 4, and quadratic, Nb ) 10, bases. Calculation using importance sampling with 104 trajectories for MDQP determined in the quadratic basis is also shown. β ) 400 is used for the importance sampling described in text.

- q0)2, and the initial momenta are defined by the full action function as p ) 2R0(x - q0). β is the time parameter equal to (or on the order of) the final propagation time. For a harmonic well these γ and β will give a uniform Gaussian distribution of trajectories according to ψ2(x, τ) at time τ ) β, thus efficiently describing the ground eigenstate. For an illustration, let us consider two identical uncoupled Morse oscillators of mass m ) 1 (the coupled case is considered in section III.C)

V(x) ) D(exp(-zx) - 1)2,

D ) 160,

z ) 1.0435 (33)

The potential mimics nonrotating H2 bond with the ZPE of 9.197 units per dimension. The harmonic oscillator energy for the quadratic expansions of this potential is 9.333 units. The initial positions of trajectories are sampled according to the quasi-random Gaussian distribution.28 Figure 3 shows the timedependence of normalized energy as a function of time obtained with 5000 trajectories for several values of β. Nonzero β values do not fully represent the wave function at τ ) 0 (thus, the difference in the short-time energy values), yet they capture anharmonicity of the Morse potential and, obviously, improve convergence of ZPE with the number of trajectories. Calculations were performed using quadratic fitting of the trajectory momenta, bf )(1, x, y, x2, y2). The energies and standard

J. Phys. Chem. C, Vol. 114, No. 48, 2010 20599 deviations are summarized in Table 1. Since the quadratic basis does not reproduce b p exactly, we do not get the ideal plateau energy value at long times. The data in the table were averaged over the period of τ ) [0.18, 0.2] units. B. Triatomic Molecules. Next, the imaginary-time quantum trajectory approach with approximate MDQP is applied to triatomic molecules, nonrotating H2O and SO2. The quartic force field potential in dimensionless normal-mode coordinates is employed for both systems.29 In the normal-mode coordinates M is a unit matrix. The initial three-dimensional wave function is a Gaussian function, which is a direct product of onedimensional functions of eq 31. The initial trajectory positions are sampled according to a uniform Gaussian distribution using Box-Muller algorithm.28 With this sampling the initial value of the action function S(x b, 0) ) S˜(x b, 0) ) 0, while the initial b A time step of dτ ) 0.2 au momentum is b p(x b, 0) ) -ψ-1∇ψ. was used for all calculations. The initial parameters of the wave function were R0 ) 0.6a0-2 and q0 ) 0.0 in all three dimensions. R0 differed from the eigenstate value of 0.5a0-2 to test the convergence properties of the formulation. To estimate ZPE we started with the linear fit to b p(x b, τ). In three dimensions the fitting basis consists of four functions, Nb ) 4, and defines constant approximate U and zero quantum force. As a test, 50 000 trajectories, propagated up to τ ) 200 au on the quadratic potential (Harmonic Approximation, HA) of H2O, gave the exact QM value, 4712 cm-1. The linear fitting is exact in this case. The QM value of the ZPE for this molecule is 4652 cm-1. Further, the ZPE was determined for the full quartic potential, which required 300 000 trajectories for convergence. The linear basis calculation gave the ZPE estimate of 4563.2 cm-1 which is 89 cm-1 lower than the exact value. The absolute difference with the exact QM is worse than the HA result, but the change is in the right direction (HA overestimates the ZPE). Thus, a more accurate fitting of b p, producing nonzero quantum force, was necessary. Inclusion of the quadratic functions into the fitting basis (Nb ) 10) significantly improved the ZPE accuracy yielding 4642.7 cm-1. The ZPE obtained with Nb ) 10 is only 9 cm-1 lower than the accurate ZPE: the quadratic basis MDQP gives 115% of the change in ZPE due to the anharmonicity of the potential. For the SO2 molecule the MDQP approach exactly reproduced the HA value of 1537 cm-1 on the harmonic surface with 30 000 trajectories propagated to τ ) 300 au. For the full quartic potential the ZPE obtained with 104 MDQP trajectories using the quadratic basis, E0 ) 1531.3 cm-1 is very close to the exact QM value of 1530 cm-1. The number of quantum trajectories needed to converge results was smaller than for H2O, because SO2 is much more harmonic than H2O.29 In both triatomic examples we needed many thousands of trajectories to obtain ZPEs, even in the harmonic approximations to the potential, because of the divergent dynamics discussed in section II.C. Remarkably, Liu and Makri obtained accurate ZPE value for water by evolving a few “independent” trajectories and the wave function derivatives through the fourth and sixth order in imaginary time (though multiple trajectories would be needed to obtain the wave function); their method is based on the amplitude and “phase” wave function representation yielding analogs of Bohm’s real time equations,19 but the amplitude/phase re-expansion keeps the trajectories nearly stationary. To make the MDQP practical, the number of trajectories has to be reduced. We used the importance sampling of section III.A for the H2O molecule. In eq 32 parameter γ was set to γ ) R0exp(ωβ/2), where ω is the fundamental frequency of each vibrational mode and β is a time parameter

20600



TABLE 1: Zero-Point Energy, E0, and Its Standard Deviation, σ, for a System of Two Morse Oscillators Obtained with Importance Samplinga Ntraj

β

E0 ( σ

β

E0 ( σ

β

E0 ( σ

1000 5000 25000

0.0 0.0 0.0

18.596 ( 0.0745 18.403 ( 0.0136 18.419 ( 0.0039

0.2 0.2 0.2

18.496 ( 0.0009 18.392 ( 0.0007 18.393 ( 0.0004

0.4 0.4 0.4

18.396 ( 0.0059 18.390 ( 0.0005 18.391 ( 0.0001

a E0 is averaged over τ ) [0.18, 0.2]. β is the convergence parameter of eq 32. The trajectory momentum is fitted with quadratic polynomials.

TABLE 2: Convergence of Zero-Point Energy, E0, of the H2O Molecule with Respect to the Number of Basis Functions, Nb, Used to Fit the Trajectory Momentum in Three Dimensionsa method -1

HA

Nb ) 4 Nb ) 10 Nb ) 16 Nb ) 20

E0 [cm ] 4712.0 4563.2

4646.3

4648.5

4650.2

QM 4652.0

a

HA is the ZPE of the harmonic approximation to the potential. Nb ) 16 contains the cubic functions involving only the most anharmonic coordinate.

on the order of the final propagation time. The calculations were repeated for the full water potential with MDQP in the quadratic basis using importance sampling. The MDQP approach has full quantum mechanics as its limit of a complete momentum-fitting basis bf. As shown earlier in one dimension21 larger bases, generally, give higher-accuracy ZPE. For H2O the polynomial fit to the trajectory momentum was extended from the quadratic to cubic basis by adding functions in the most anharmonic mode, x -mode, first. Addition of all the cubic basis functions containing x to the momentum fit (Nb ) 16) gave ZPE of 4648.5 cm-1, about 2 cm-1 higher than the quadratic fit. The z mode is also appreciably anharmonic and strongly coupled to the x mode. Addition of the remaining cubic terms (Nb ) 20) gave the ZPE of 4650.2 cm-1. There is an improvement in ZPE of approximately 4 cm-1 (out of 6 cm-1) compared to the quadratic basis calculation. Convergence of the ZPE for water with the basis size is summarized in Table 2. The exact QM limit of complete fitting basis was not pursued as impractical in high dimensions. For H2O the importance sampling reduced the number of the MDQP trajectories by a factor of 30 to 10000, which is still large for a three-dimensional system. This is explained by the long propagation time needed for the initial wave function to decay to the ground state because this system is rather anharmonic. For large τ the overall divergent dynamics of the Lagrangian trajectories results in large number of trajectories (even with importance sampling) as most of them leave the are of the high wave function density. A solution to this problem is to work in the Eulerian frame-of-reference described below. C. Dynamics in the Eulerian Frame of Reference. The importance sampling, compensating for the divergent dynamics of trajectories in the harmonic approximation to a bound potential, improves convergence of the ZPE with the number of trajectories in the Lagrangian frame-of-reference. This approach is convenient for the imaginary-time evolution of scattering systems and seamlessly connects to the real-time dynamics of the quantum trajectories, a feature important for the thermal reaction rate calculations.22 In real time the quantum trajectories track the flow of the wave function density, thus, defining the optimal “moving grid” for the time-dependent wave function. However, in the imaginary-time calculations of ZPE, moving trajectories is not an optimal wave function representation, because the target of the calculation, i.e., the ground state wave function is a stationary in space object. In this case the

Eulerian frame of stationary “trajectories” is a sensible approach. This approach (and a combination of the Lagrangian/Eulerian grids) has been used for the real-time quantum trajectory dynamics30 and in imaginary time (“zero-velocity trajectories”).20 Liu and Makri19 used wave function re-expansion to achieve nearly stationary trajectories in the imaginary-time ZPE calculations. In this section we combine the Eulerian frame of reference for the imaginary-time wave function dynamics with the global approximation to the quantum potential of section II.B. A real positive wave function given by eq 6, solving the diffusion eq 1 is represented by an ensemble of random points in b x. The evolution of S(x b, τ) is given by eq 7, which, defining the momentum as before, b p ) ∇S(x b, τ), is

∂S p 1 p + V + ∇ · M-1 · b p )- b p · M-1 · b ∂τ 2 2

(34)

Taking the gradient of eq 34 the time-dependence of b p in the Eulerian frame of reference

∂b p ) ∇(V + U0) ∂τ

(35)

Function U0 includes the quantum potential and the kinetic energy-like term appearing due to the choice of the frameof-reference

U0 )

p 1 p - b p ∇ · M-1 · b p · M-1 · b 2 2

(36)

For practical reasons U0 will be determined from the global fitting of the momentum described in section II.B. For consistency with the definition of b p ) ∇S(x b, τ), the approximate U0 will be also used in eq 34 rewritten as

∂S ) V + U0 ∂τ

(37)

Equations 37 and 35 with U0 of eq 36 determined from the Least Squares Fit to b p give the time-dependent wave function (6). The accuracy of the solution is determined by the accuracy of the momentum fitting, i.e., by the size of the basis bf(x b) as in the Lagrangian frame-of-reference formulation. The approach, tested on a system of 40 linearly coupled harmonic oscillators, reproduced analytical ZPE values with very sparse wave function sampling, only 50 points. Good convergence of the ZPE with respect to the number of trajectories is attributed to the fact that energy of each sampling point b x becomes equal to the ZPE. Therefore, in the normalized energy the norm of the wave function cancels, regardless of how well the wave function itself is described.


J. Phys. Chem. C, Vol. 114, No. 48, 2010 20601

TABLE 3: Zero-Point Energy for the Two-Dimensional Morse Oscillator System with Coupling Computed in the Eulerian Formulationa Coupling

QM

Ntraj ) 50

V12

β

E0

E0

Efit 0

0 0.25 0.5 0 0 0

0 0 0 0.5 0.75 1.0

18.395 20.400 22.085 17.784 16.807 13.285

18.312 20.331 22.025 17.702 16.723 13.212

18.313 20.332 22.027 17.703 16.725 13.267

Ntraj ) 5000 E0

Efit 0

∆

18.347 20.358 22.047 17.736 16.758 13.159

18.349 20.360 22.049 17.739 16.760 13.176

0.006 0.004 0.004 0.002 0.002 0.002

a The approximate MDQP is obtained from the quadratic fit to b p. The imaginary-time propagation was performed up to τ ) 0.5. Efit 0 values are obtained by analyzing the wavefunction norm, 〈ψ|ψ〉τ, over the time interval τ ) [0.4, 0.5] as described in text. ∆ is the difference between E0 obtained with 5000 and 500 trajectories (the same for E0 and Efit 0 ). For the bottom row the propagation time was τ ) 2.0, the time interval, from which Efit 0 was obtained, τ ) [1.6, 2.0], and ∆ is given in comparison to the calculation with 25 000 trajectories. The harmonic ZPE without coupling is 18.667 units.

As a more thorough test the Eulerian evolution with momentum fitting to has been applied to the two-dimensional Morse oscillator model of Section III A: the convergence with respect to the number of points was dramatically improved. The action function computed at 50 points, sampling a two-dimensional Gaussian given by eq 31 in each dimensions, with the quadratic fit of b p gave the ZPE estimate of E0 ) 18.300 units. This value should be compared to the exact QM value of 18.394 and to the harmonic oscillator value of 18.667. As seen from the first row of data in Table 3, the Eulerian ZPE estimate converges with 100 times fewer trajectories but is less accurate than the Lagrangian ZPE estimate (see Table 1) for the same initial wave function, R0 ) 9.0a0-2 and q0 ) 0.0. We believe there is a physical and a formal reason for the difference in the Lagrangian and Eulerian estimates using the same fitting basis. The physical reason is that fixed sampling points of the Eulerian formulation do not “explore” the potential as do the Lagrangian trajectories and, thus, less information on anharmonicity of the potential influences the ZPE value. Starting with the more delocalized wavepacket, R0 ) 4.5a0-2, gave the Eulerian ZPE of 18.37 units, which is closer to the exact QM result. The formal reason is that in the Eulerian formulation more quantities are approximated. U0 contains not just the pM-1 -dependent quantum p/2, potential, but also the kinetic energy-like term, b p · pM-1 · b which does not vanish in the classical limit of infinitely heavy particle. Perhaps, a more balanced use of exact and fitted momenta in the approximate Eulerian evolution eqs 37 and 35 can be developed. We further tested the Eulerian MDQP approach in the presence of coupling considering both kinetic and potential energy coupling. For the former the cross-terms describing bondcoupling were included into the kinetic energy

M-1 )

(

1 -β -β 1

)

(38)

For the latter we add the coupling term, V12D(x - y)2, to the Morse potentials in x and y. The results from the dynamics with 50 and 5000 trajectories, obtained now with the addition of xy function into the fitting basis, are given in Table 3. The discrepancy with the 500 trajectory calculation is also shown. The effect of coupling on the ZPE was a shift by 3-28% of

the uncoupled QM value. In all cases the quadratic approximation gave the ZPE with accuracy better than 0.5% with 50 trajectories and better than 0.3% with 5000 trajectories (with the exception of β ) 1 case, for which the accuracy is 0.9%). Since the average energy of a wavepacket does not reach a perfect plateau at long τ, we also examined ZPEs using eq 20 by fitting ln〈ψ|ψ〉τ to a straight line for the last 20% of the propagation time. These values, given in the table, are in close agreement with the final average energies of the wavepackets. The presence of the bond-coupling coupling, nonzero β, did not affect the accuracy of the method. The ZPE estimates for nonzero potential coupling V12 also were of about the same accuracy for moderate values of V12 considered. These results indicate that the MDQP approach at the quadratic level of approximation performed equally well with and without coupling raising expectations that the approach will work in Cartesian coordinates. IV. Conclusions We have described the multidimensional formulation of the Schro¨dinger equation in imaginary time based on the exponential representation of a positive real wave function, ψ(x b, τ) ) exp(-S(x b, τ)/p). In the Lagrangian frame of reference, defining the trajectory momentum as b p ) ∇S(x b, τ) as in classical mechanics, the wave function evolution is representable in terms of a trajectory ensemble evolving on the inverted classical potential, -V, to which the momentum dependent quantum p/2, is added. This particular potential (MDQP), U ) p∇ · M-1 · b form of the quantum potential, responsible for all QM effects in the system, is the main difference from the real-time quantum trajectory formulation of SE.18 Analytical definition of the MDQP and associated with it force from the global least squares fit to b p using small Taylor basis enables cheap estimation of QM effects, such as the zeropoint energy, in high dimensions. The idea of global approximation to MDQP is similar to the approximate quantum potential (defined by the wave function amplitude in eq 4 in real-time quantum trajectory dynamics31), but there are important differences. In real time the fitting procedure was defined in terms of the moments of the trajectory distribution over the entire space (because of integration by parts), whereas for the imaginarytime case the fitting procedure is rigorous for any discrete set of trajectories. This feature enabled us to use importance sampling and, in principle, allows semilocal or even local fitting of b p which would increase the accuracy of the method. Another difference between the imaginary and real-time dynamics is that the interference effects, responsible for unstable behavior of the quantum trajectories in real time, are quenched in imaginary time. Smooth imaginary-time trajectories imply relatively smooth functions b p allowing us to obtain ZPE estimates using quadratic fitting basis, its size scaling quadratically with dimensionality. As shown for water, addition of the cubic basis functions (possibly into just the strongly anharmonic modes) increases the accuracy of the approximate calculation. The globally defined MDQP also distinguishes this work from the “independent” quantum trajectory approaches to the imaginarytime dynamics of refs 19 and 20 We also observed that for multidimensional systems in imaginary time the overall divergent trajectory dynamics in bound potentials (trajectories are rolling off potential barriers, -V) leads to undersampling of a wave function at long times. Consequently, many thousands of trajectories uniformly distributed at the start of evolution had to be propagated to reproduce the final wave function already in two-three dimen-

20602


sions. To reduce the number of trajectories the importance sampling, guided by the harmonic ZPE estimate and final propagation time (see eq 32), was introduced using the properties of the MDQP fitting procedure. The importance sampling considerably improved convergence of the ZPE with respect to the number of trajectories. In the H2O calculation the number of trajectories was reduced by a factor of 30, though still 104 trajectories were needed. In the ZPE calculations, a promising strategy of reducing the number of trajectories is to use the Euler frame-of-reference and not propagate trajectories alltogether, as has been done in zero velocity complex action approach.20 The evolution eqs 37 and 35 have been solved in the Eulerian frame of reference using, once again, the global least squares fit to b p in the quadratic basis. The approximation to the elements quadratic in b p in eq 36 was used to maintain the correct relation between the wave function norm and energy given by eq 20. Evolution with just 50 points gave the ZPE estimates of accuracy better than 0.5% for a model two-dimensional system of Morse oscillators with kinetic energy coupling and with the bilinear potential coupling. (The same number of points reproduced exact ZPE for 40 bilinearly coupled harmonic oscillators.) This is explained by the fact that for the wave function given by eq 6 energies of all trajectories become equal to the ZPE; further reduction in the number of trajectories could be achieved with more accurate fitting of the momentum, deserving further investigation. Although the Lagrangian formulation is well suited for scattering imaginary-time dynamics (or delocalized ground state calculations), encountered in calculations of thermal rates,22 the ZPE calculations Eulerian formulation has merits. An eigenstate does not move, thus, it can be well represented by stationary points b x, if reasonable estimates for the location and localization of the ground state are available. In addition, the classical potential and forces have to be computed only once (!), rather than at each time step, which would be a big computational saving in direct dynamics calculations. Both Lagrangian and Eulerian formulations are compatible with the calculation of low-lying excited energy eigenstates as have been described previously21 for the former. Ongoing work in the group is focused on the ZPE calculations for larger molecular systems and multidimensional calculations of thermal reaction rates.

Garashchuk and Vazhappilly Acknowledgment. This work was supported by the University of South Carolina. References and Notes (1) Janecek, S.; Krotscheck, E. Comput. Phys. Commun. 2008, 178, 835. (2) Blume, D.; Lewerenz, M.; Niyaz, P.; Whaley, K. B. Phys. ReV. E 1997, 55, 3664. (3) Ceperley, D. M. Mitas, L. Monte Carlo methods in quantum chemistry; AdVances in Chemical Physics; Wiley: New York, 1996; Vol. 93. (4) Lester, W. A., Jr.; Mitas, L.; Hammond, B. Chem. Phys. Lett. 2009, 478, 1. (5) Viel, A.; Coutinho-Neto, M. D.; Manthe, U. J. Chem. Phys. 2007, 126, 024308. (6) Hinkle, C. E.; McCoy, A. B. J. Phys. Chem. A 2008, 112, 2058. (7) Miller, W. H. J. Chem. Phys. 1971, 55, 3146. (8) Makri, N.; Miller, W. H. J. Chem. Phys. 2002, 116, 9207. (9) Miller, W. H.; Schwartz, S. D.; Tromp, J. W. J. Chem. Phys. 1983, 79, 4889. (10) Light, J. C.; Carrington, T., Jr. AdV. Chem. Phys. 2000, 114, 263. (11) Miller, W. H. J. Phys. Chem. A 2001, 105, 2942. (12) Garashchuk, S.; Rassolov, V. A. Chem. Phys. Lett. 2002, 364, 562. (13) Donoso, A.; Yeng, Y. J.; Martens, C. C. J. Chem. Phys. 2003, 119, 5010. (14) Maddox, J. B.; Bittner, E. R. J. Chem. Phys. 2003, 119, 6465. (15) Liu, J.; Makri, N. J. Phys. Chem. A 2004, 108, 5408. (16) Trahan, C. J.; Hughes, K.; Wyatt, R. E. J. Chem. Phys. 2003, 118, 9911. (17) Garashchuk, S.; Rassolov, V. A. J. Chem. Phys. 2004, 120, 1181. (18) Bohm, D. Phys. ReV. 1952, 85, 166. (19) Liu, J.; Makri, N. Mol. Phys. 2005, 103, 1083. (20) Goldfarb, Y.; Degani, I.; Tannor, D. J. Chem. Phys. 2007, 338, 106. (21) Garashchuk, S. J. Chem. Phys. 2010, 132, 014112. (22) Garashchuk, S. Chem. Phys. Lett. 2010, 491, 96. (23) Rassolov, V. A.; Garashchuk, S.; Schatz, G. C. J. Phys. Chem. A 2006, 110, 5530. (24) Wyatt, R. E., Quantum Dynamics with Trajectories: Introduction to Quantum Hydrodynamics; Springer-Verlag: Berlin, 2005). (25) Goldfarb, Y.; Degani, I.; Tannor, D. J. J. Chem. Phys. 2006, 125, 231103. (26) Rassolov, V. A.; Garashchuk, S. J. Chem. Phys. 2004, 120, 6815. (27) Lopreore, C. L.; Wyatt, R. E. Phys. ReV. Lett. 1999, 82, 5190. (28) Press, W. H.; Flannery, B. P.; Teukolsky, S. A.; Vetterling, W. T. Numerical Recipes: The Art of Scientific Computing, 2nd ed.; Cambridge University Press: Cambridge, U.K., 1992. (29) Isaacson, A. D.; Truhlar, D. G.; Scanlon, K.; Overend, J. J. Chem. Phys. 1981, 75, 3017. (30) Trahan, C. J.; Wyatt, R. E. J. Chem. Phys. 2003, 118, 4784. (31) Garashchuk, S.; Rassolov, V. A. Chem. Phys. Lett. 2003, 376, 358.

JP1050244

Multidimensional Quantum Trajectory Dynamics in Imaginary Time

Recommend Documents