THE JOURNAL OF CHEMICAL PHYSICS 141, 234107 (2014)

Dynamics in the quantum/classical limit based on selective use of the quantum potential Sophya Garashchuk,a) David Dell’Angelo, and Vitaly A. Rassolov Department of Chemistry and Biochemistry, University of South Carolina, Columbia, South Carolina 29208, USA

(Received 15 August 2014; accepted 26 November 2014; published online 16 December 2014) A classical limit of quantum dynamics can be defined by compensation of the quantum potential in the time-dependent Schrödinger equation. The quantum potential is a non-local quantity, defined in the trajectory-based form of the Schrödinger equation, due to Madelung, de Broglie, and Bohm, which formally generates the quantum-mechanical features in dynamics. Selective inclusion of the quantum potential for the degrees of freedom deemed “quantum,” defines a hybrid quantum/classical dynamics, appropriate for molecular systems comprised of light and heavy nuclei. The wavefunction is associated with all of the nuclei, and the Ehrenfest, or mean-field, averaging of the force acting on the classical degrees of freedom, typical of the mixed quantum/classical methods, is avoided. The hybrid approach is used to examine evolution of light/heavy systems in the harmonic and double-well potentials, using conventional grid-based and approximate quantum-trajectory time propagation. The approximate quantum force is defined on spatial domains, which removes unphysical coupling of the wavefunction fragments corresponding to distinct classical channels or configurations. The quantum potential, associated with the quantum particle, generates forces acting on both quantum and classical particles to describe the backreaction. © 2014 AIP Publishing LLC. [http://dx.doi.org/10.1063/1.4903764] I. INTRODUCTION

In general, the cost of numerically solving the Schrödinger equation for coupled anharmonic systems scales exponentially with the system size, and the cost of storing the wavefunction alone may render exact quantum-mechanical study unfeasible for large systems.1, 2 However, for heavy particles, such as nuclei, the classical dynamics is often accurate and insightful,3 and the cost of solving the Newton’s equations of motion scales linearly with the dimensionality. (For realistic chemical systems, however, the overall cost is likely to be dominated by evaluation and scaling of the potential.) Classical dynamics methods are used to study systems of thousands of atoms,4 whereas a full-dimensional quantum study of a pentatomic reaction, a collision of hydrogen and methane,5, 6 is the largest reactive scattering study to date. We are interested in dynamics of molecular systems in the intermediate regime when the classical dynamics captures the general behavior, but the quantum-mechanical (QM) effects, such as tunneling and the zero-point energy (ZPE), are important in a few degrees of freedom (DOFs) describing, for example, protons at low energy or temperature. Development of a rigorous theoretical framework for evolution of interacting quantum and classical particles is a fundamental problem of theoretical chemistry and physics with long rich history (beyond our scope). For the purposes of this article, it is important to distinguish two separate aspects of the description of quantum-classical system: (i) approximation of the Hamiltonian, and (ii) approximation of the wavefunction. a) [email protected]

0021-9606/2014/141(23)/234107/12/$30.00

It is common to consider just the Hamiltonian approximations, i.e., equations of motion for the quantum and classical particles, with “wavefunction” for the classical degrees of freedom collapsing to a single point. For example, the “classic” solution is the Ehrenfest or mean-field approach, where evolution of classical particles is represented with a single trajectory evolving according to the forces averaged over a wavefunction describing the quantum particles.7 The main problem of the Ehrenfest-style quantum/classical approaches8 is the absence of the backreaction of the quantum DOFs (QDOFs) on the classical DOFs (CDOFs), manifested in a single outcome, i.e., “no branching” for the CDOFs. Representation of the classical particles with multiple trajectories in the CDOFs allows multiple outcomes in a statistical sense, but precludes interaction between the processes associated with different mean-field trajectories, which may be important, for example, for reactions with multiple paths. Other approaches of dealing with the backreaction problem are the multiconfiguration time-dependent mean field approach9, 10 and the surfacehopping.11 The former approach requires careful analysis of a system to determine a compact practical basis; the latter approach is developed for a few discrete quantum states among which trajectories hop according to heuristic prescriptions. We consider multidimensional dynamics involving a continuum of states in both, QDOFs and CDOFs, and leading to multiple outcomes in the CDOFs, exemplified by reactive processes occurring in the condensed phase. In this regime, semiclassical and quasiclassical trajectory approaches12, 13 – based on classical trajectory evolution from which the QM effects are assessed – are particularly attractive for practical reasons. The trajectory description of heavy particles is efficient while

141, 234107-1

© 2014 AIP Publishing LLC

234107-2

Garashchuk, Dell’Angelo, and Rassolov

the basis representation of highly oscillatory (in the classical limit of ¯ → 0) wavefunctions is not. The initial conditions for the trajectories can be chosen randomly, circumventing exponential scaling typical of conventional direct-product basis methods of quantum dynamics.14 Newton’s equations of motion are easy to solve, especially in Cartesian space. We work with a trajectory formulation of the time-dependent Schrödinger equation (TDSE), due to Madelung, de Broglie, and Bohm,15–17 which is conceptually appealing. A wavefunction is described by an ensemble of quantum trajectories (QTs) whose dynamics exactly accounts for all QM effects. While generally impractical as an exact numerical approach, this QT formulation, also called Bohmian or hydrodynamic, provides a universal framework for a hybrid quantum/classical description of light “quantum” and heavy “classical” particles comprising a large molecular system. The QT formulation of the SE has been used by Burant and Tully to develop nonadiabatic dynamics,18 where the light particle is an electron whose wavefunction is represented in a timeindependent basis of electronic states. In the context of timedependent wavepackets (continuum of QM states), the QT formulation has been used in the mixed quantum-classical approaches8, 19 based on the quantum evolution of multiple, reduced dimensionality wavefunctions in the QDOFs, guided by the CDOF trajectories evolving without the quantum potential, independently of each other. A closely related approach, the energy-conserving mixed quantum/classical trajectory dynamics,20 is based on the wavefunction factorization into a fast “quantum” and slow “classical” components. As a result, a full-dimensional wavefunction is represented by multiple QT sub-ensembles, each one being guided by a mean-field trajectory in the CDOFs similar to methods of Refs. 8 and 19. The energy of the QT sub-ensembles is rigorously conserved and their independence of each other is convenient for implementation enabling, for example, studies of collisions of H with C37 H15 .21 However, inability of different sub-ensembles to exchange probability density and energy may be a serious limitation if the quantum effects are significant, for example, when the vibrational ZPE and energy flow affect the bond breaking.22–24 In this paper, we consider a classical limit of quantum dynamics defined by compensating for the quantum potential in the TDSE. Selective use of the quantum potential for the light/heavy particles yields a hybrid quantum/classical dynamics of a full-dimensional wavefunction without any restrictions on its form and without the mean-field approximation. Formulated in terms of QTs, this hybrid dynamics of a wavefunction is based on a single ensemble of trajectories moving in full space of QDOFs and CDOFs with the quantum forces derived only from the QDOFs and acting on both QDOFs and CDOFs. The quantum corrections are determined approximately in full dimensionality on timeindependent spatial domains determined by the positions of heavy particles. This makes the approach practical and, importantly, precludes unphysical interactions between the fragments of a total wavefunction corresponding to distinct classical channels or configurations. To illustrate the proposed method and to contrast it with other methods, let us consider an evolution of a Gaussian

J. Chem. Phys. 141, 234107 (2014)

wavepacket, ψ = exp ( − α(x − y)2 − β(y − y0 )2 − ıp0 (y − y0 )), in a potential V = k(x − y)2 /2. The light particle of mass m1 is described by the coordinate x, and the heavy particle of mass m2 is described by the coordinate y. This is a simple model of an “atom” propagating in free space, with an “electron” (x, m1 ) bound to a “nucleus” (y, m2 ) by a harmonic potential. From the onset, we will consider the motion in y coordinate to be fully classical, and the motion in x to be exact, i.e., fully quantum. If at a later time such a system were to encounter a potential barrier in the classical coordinate, the probability of passing through the barrier would be determined by the kinetic energy in the y-coordinate. Therefore, let us consider the kinetic energy in y DOF during the propagation. Regardless of the specific approximation on the classical motion, the evolution of the quantum part of the system is largely determined by the Gaussian “width” parameter α,  whose coherent value is αc = k/m1 /2. In the non-coherent regime, α = α c , there are oscillations of the wavepacket width, due to the time-dependence of α, accompanied by the translational motion of the wavepacket manifested in the time-dependence of y0 , with possible evolution of other parameters depending on approximations. In the coherent regime, the oscillation of the wavepacket may vanish. Now, let us compare the proposed hybrid quantum/classical approach with three other quantum-classical methods, most closely related to it. In the (i) Ehrenfest approach,7 β → ∞ and the system becomes a quantum one-dimensional harmonic oscillator with uniform classical motion in the y coordinate, regardless of α. In the (ii) closely related approaches of Prezhdo and Brooksby8 and Gindensperger et al.,19 the system is subdivided into the ensemble of trajectories in y, each associated with the quantum wavefunction in x. In the limit of fully classical motion in y, each of these trajectories behaves like an Ehrenfest trajectory that we have just considered. The (iii) mixed quantum trajectory dynamics with approximation on domains of Ref. 20 will describe some energy transfer between the quantum and classical DOFs for any value of α. Finally, the proposed approach will describe the energy transfer between the quantum and classical DOFs in the noncoherent α regime, and the translational motion of the otherwise unchanged wavepacket for the coherent α = α c . We believe this to be the most realistic description of evolution of a quantum/classical system, with the energy transfer from the quantum to classical DOFs describing a backreaction due to quantum evolution. This backreaction should vanish when the quantum system is in a stationary state, as in the coherent α = α c regime. The hybrid quantum/classical formalism is presented in Sec. II, and is applied to dynamics of light/heavy systems evolving in the harmonic and double-well potentials, using exact QM and approximate QT time propagation in Sec. III. Section IV concludes. II. THEORY A. The quantum trajectory formalism

We will define the classical limit of the time-dependent Schrödinger equation obtained from its trajectory or

234107-3

Garashchuk, Dell’Angelo, and Rassolov

J. Chem. Phys. 141, 234107 (2014)

hydrodynamics formulation due to Madelung, de Broglie, and Bohm.15–17 This formulation is restated here for further reference. The QT formalism is presented for a general multidimensional quantum system described in full Cartesian space, followed by an adaptation for a molecular system comprised of light and heavy nuclei. The positions of Np quantum particles are sequentially arranged into a column-vector z of dimensionality Nd = 3Np . Assuming for simplicity a diagonal form of the kinetic energy operator, the particle masses, {m1 , m2 . . . mN }, define a diagonal matrix G as p

G11 = G22 = G33 = m−1 1 ,

G44 = G55 = G66 = m−1 2 .... (1) The range of integration for each spatial coordinate is [ − ∞, ∞]. The gradient, ∇ z , is a column-vector of spatial derivatives, ⎛ ⎞ ⎞ ⎛ ∂/∂z1 z1 ⎜ ⎟ ⎟ ⎜ ⎜ ∂/∂z2 ⎟ ⎜ z2 ⎟ ⎜ ⎟ ⎟ ⎜ (2) z = ⎜ . ⎟ and ∇z = ⎜ . ⎟ . ⎜ .. ⎟ ⎜ .. ⎟ ⎝ ⎠ ⎠ ⎝ zN ∂/∂zN d

d

The QT formalism for arbitrary coordinate systems and kinetic energy operators has been given elsewhere.25 Wavefunction evolution is governed by the TDSE with the time-independent potential, V , which will be referred to as “classical,” ı¯

∂ ψ(z, t) = Hˆ ψ(z, t). ∂t

(3)

The operator Hˆ is the usual QM Hamiltonian,

of the classical potential V and the quantum potential U, defined as U =−

¯2  T ¯2 T ∇z G∇z A = − r Gr + ∇zT Gr , 2A 2

where r is the non-classical component of the momentum operator,26 r = A−1 ∇z A.

The wavefunction amplitude A can be reconstructed using the continuity equation for the probability density, A2 = |ψ|2 : the probability of finding a particle within the volume element dν t ,

p = ∇z S,

(6)

and switching to the Lagrangian frame-of-reference, ∂ d = + pT G∇z , (7) dt ∂t one obtains the following 2Nd equations of motion for the quantum trajectories: dzt = Gpt , dt

dνt = dz1t dz2t . . . dzN t ,

(4)

The superscript T denotes a transposition, so that the ∇ T is a row-vector. The trajectory formulation is based on the polar representation of a wavefunction,

ı S(z, t) (5) ψ(z, t) = A(z, t) exp ¯ in terms of two scalar functions – a non-negative amplitude A and a real phase S. Substituting Eq. (5) into the TDSE of Eqs. (3) and (4), identifying the trajectory momentum p with the gradient of the phase,

(8)

d pt = − ∇(V + U )|z=z . (9) t dt Subscript t labels quantities associated with trajectories. Each trajectory is defined by its position zt and momentum pt . The dynamics of the QTs is governed by the combined influence

(11)

The potential U is a non-local time-dependent function which couples evolution of multiple QTs representing a single wavefunction. Interaction of QTs mediated by the quantum potential is responsible for all QM effects in a system. Note that the QT formulation is formally equivalent to the TDSE and has a natural classical dynamics limit of U = 0. The limit of vanishingly small quantum potential (¯ → 0 or of the infinite particle mass), is appropriate for heavy particles such as nuclei and will be invoked in Sec. II B. The QTs representing a single wavefunction ψ(z, t) form an ensemble: all trajectories should be propagated concurrently if the quantum potential is evaluated on-the-fly. In principle, an ensemble of QTs contains complete information about a wavefunction. The wavefunction phase St at a trajectory position zt is given by 

t 1 T pτ Gpτ − ( V + U )|z=z dτ. (12) St = τ 2 0

2

¯ Hˆ = − ∇zT G∇z + V (z). 2

(10)

d

(13)

associated with each trajectory, i.e., a trajectory weight wt is constant in time,26 wt = (A(zt ))2 dνt ,

dwt /dt = 0.

(14)

Equation (14) allows calculation of expectation values of the ˆ over the position- and momentum-dependent operators, , QT ensemble without explicit knowledge of the wavefunction amplitude,

ˆ t= 

Ntraj

. . . |ψ(zt )| (zt , pt )dνt ≈ 2

  (k) (k)  zt , pt w (k) . k=1

(15) In Eq. (15), the multidimensional integration over the coordinate space is replaced with the summation over the ensemble of a finite number of trajectories, Ntraj , representing an initial wavefunction. The time-independent weight of the kth trajectory, w (k) , is defined by an initial wavefunction and a sampling scheme of the trajectory positions. Equation (15) is used to develop practical approximations to the quantum potential. Computational cost of quantum dynamics for general (anharmonic, coupled) problems scales exponentially with the system size due to inherent complexity of wavefunctions. Recasting TDSE in terms of trajectories does not change this complexity. Since the cost of classical dynamics scales linearly with the system size, the unfavorable scaling of exact

234107-4

Garashchuk, Dell’Angelo, and Rassolov

QT dynamics comes from the quantum potential U.2 Therefore, approximate and computationally efficient treatment of the quantum potential is required for a practical method. One of the cheapest approaches is to use a global fit to the nonclassical momentum defined by Eq. (11), thus assuming its identical functional form for all trajectories in the distribution. In its simplest form relevant to quantum molecular dynamics, the linear fitting function is used, yielding quadratic form of the quantum potential, and linear quantum force according to Eq. (10). The method, termed the Linearized Quantum Force (LQF),26 is outlined in the Appendix. One of the shortcomings of the LQF approximation is unphysical interaction of isolated wavefunction fragments through the common global fitting procedure. A simple fix is to perform the fits to the non-classical momentum on separate spatial domains.27 This approach has added cost due to commutators of the kinetic energy and the domain functions, which have to be included into the definition of the quantum potential. For a system combining classical and quantum DOFs, however, the domain formalism simplifies if the domains are defined only in the CDOFs, because the commutators of the kinetic energy and domain functions vanish. At the same time, the unphysical long-distance interactions disappear in systems for which strong delocalization in the QDOF is accompanied by significant changes in the CDOFs. This is the case of a typical proton transfer reaction in a biochemical system, with quantum description of the reactive proton, or in a Marcus-type electron transfer reaction in solution, with quantum treatment of the transferring electron.28 In this work, we develop a formalism and investigate its performance in the model systems. The combined quantum-classical TDSE is introduced in Sec. II B, and the CDOF domains are described in Sec. II C. Atomic units, specifically ¯ = 1, will be used in the remainder of the paper. B. The quantum/classical TDSE in full space

In the QT framework, a natural (if not mathematically rigorous) transition to classical mechanics is to set U ≡ 0, which provides a “quasiclassical” description of the system: a particle is still described by a wavefunction, represented with trajectories carrying amplitudes and phases, but the trajectories evolve under the influence of the classical potential independently of each other. Therefore, a hybrid quantum/classical dynamics can be defined by neglecting contributions to the quantum potential, associated with the DOFs describing heavy, nearly classical particles. A straightforward way is to set to zero all terms involving derivatives over the classical DOFs in Eq. (10). Below, instead of a single position vector z, we use x to denote the coordinates of Nq quantum dimensions, describing Nq /3 light particles, and y to label the remaining Nd − Nq CDOFs describing (Nd − Nq )/3 heavy particles. Subscripts x denote quantities associated with the QDOFs, and y with the CDOFs. The full quantum potential of Eq. (10) is divided into contributions from the quantum and classical DOFs, U = Ux + Uy 1 Ux = − (rx + ∇x )T Gx rx , 2

1 Uy = − (ry + ∇y )T Gy ry . 2 (16)

J. Chem. Phys. 141, 234107 (2014)

The matrices Gx and Gy are the diagonal matrices of inverse masses for the quantum and classical DOFs. To obtain the hybrid quantum/classical limit Uy is simply set to zero, while Ux is retained, Uy ≈ 0,

U ≈ Ux .

(17)

The same quantum/classical dynamics is obtained by subtracting the quantum potential associated with the classical DOFs from the full-dimensional TDSE, ∂ ψ ı ψ( x , y, t) = Hˆ ψ( x , y, t) + ∇yT Gy ∇y |ψ|. (18) ∂t 2|ψ| In the absence of the CDOF, Eq. (18) is reduced to the standard TDSE in the x-coordinate; in the absence of the QDOF Eq. (18) describes dynamics of a quasiclassical wavefunction.29 Numerical solution of Eq. (18) as well as computation of exact Bohmian trajectories is challenging, due to, generally unstable dynamics near the focal points (or when dispersion is small) in CDOFs as discussed in Sec. III, but it will be used for benchmarking of the approximate QT implementation at short times. The concept of the quantum potential leads the hybrid quantum/classical TDSE (18) which has no assumptions or approximations regarding the wavefunction, and in this respect is different from the Mixed Quantum Classical Bohmian (MQCB) dynamics of Gindensperger, Meier, and Beswick19 and of, essentially identical to it, Prezhdo and Brooksby.8 In MQCB dynamics, the quantum potential in the CDOFs is neglected as well as in Eq. (18), but the central approximation of MQCB dynamics is the mean-field coupling of the QDOFs and CDOFs, absent in the hybrid quantum/classical TDSE. In MQCB, this decoupling is achieved by the meanfield-type averaging of the classical force acting on the classical degrees of freedom over a fragment of the wavefunction defined in QDOFs. Multiple such pieces – each guided by a mean-field trajectory – are evolved in QDOFs only. The branching problem of the Ehrenfest dynamics is solved by, essentially, sampling of the probability density in the CDOFs with the Ehrenfest-type trajectories, but there is no force acting from the QDOFs on the CDOFs beyond the mean-field. As shown in Ref. 20 the MQCB approximation is formally traced to a factorization of a full-dimensional wavefunction into slow component, approximated with the mean-field trajectories in the CDOFs, and a fast component, propagated in time using exact or approximate methods of quantum dynamics. While having multiple, independent of each other, fragments of a wavefunction in QDOF (each guided by its own CDOF trajectory) is advantageous for numerical work and a quantum/quasiclassical wavefunction may be defined in fulldimensional space, the fact that these fragments cannot exchange energy or probability density, is an inherent limitation of MQCB dynamics. Equation (18) describes evolution of a wavefunction in full dimension (QDOF and CDOF) with no constraints or assumption about the wavefunction form. C. The quantum/classical dynamics in the trajectory formulation with domains in classical DOFs

Our implementation of Eq. (18) in the trajectory framework is based on the LQF in which the nonclassical

234107-5

Garashchuk, Dell’Angelo, and Rassolov

J. Chem. Phys. 141, 234107 (2014)

momentum r is approximated by linear functions of coordinates (Appendix). The LQF is defined variationally on full space. Thus, the corresponding quadratic quantum potential and linear quantum force are of the mean-field character, yet there is no approximation made to the wavefunction itself. Representation of the wavefunction in polar form, Eq. (5), however, is inefficient for a wavefunction with multiple isolated fragments and will result in a very complicated r. In this case, the globally defined approximation to r is very inaccurate, similar to approximating a sum of two nonoverlapping Gaussian functions with a single Gaussian. To alleviate this problem it has been suggested to approximate r semi-globally, using domains or subspaces in the CDOFs,30 yielding an approximate U of a multiconfigurational meanfield character. For example, if a CDOF describes a reaction coordinate, based on the physics of the problem the domains can be associated with the reactant, product, and interaction regions. This setup prevents unphysical interaction of spatially separated wavefunction fragments and yields more accurate LQF. This idea has been implemented for the LQF dependent only on the QDOF coordinates for model two-dimensional systems of reactive and non-reactive scattering.31 The procedure has been shown to converge with respect to the number of domains, but large number of domains and additional stabilization forces were necessary to achieve acceptable energy conservation. The convergence/stability issues were attributed to large changes in the quantum force parameters from domain to domain. In this work, we approximate the QDOF components of the nonclassical momentum, rx , in a linear basis involving all coordinates. This improves convergence and enables use of larger domains and studies of high-dimensional systems (up to 14 DOFs). The LQF acting in the lth domain is determined by the trajectories belonging y ). The to this domain, defined by the function l ≡ l ( resulting quantum potential,  1 U≈ Ul ( x , y)l ( y ), Ul = − (r˜ x,l +∇x )Gx r˜ x,l (19) 2 l generates a quantum force acting on the QDOFs and on the CDOFs. Each Ul is determined variationally, and the total energy is formally conserved. For simplicity let us consider domains, associated with a single coordinate y (correlating, for example, with the separation of the reactive fragments). The domains (a total of L domains) are defined by a smooth representation of the step-function, such as the Fermi or logistic function. Labeling the domain boundaries with Yl , l = 0. . . L, all trajectories characterized by positions ( xk , yk ) contribute mostly to the lth domain, if Yl − 1 < yk < Yl . Contribution to the neighboring domains vanishes exponentially. The domain functions add up to 1, l l ≡ 1. The quantum potential of Eq. (19) exerts a force on the QDOFs and on the CDOFs  ∇x Ul l (y), (20) ∇x U = l

∇y U =

 l

∇y Ul l (y) +

 l

Ul ∇y l (y).

(21)

The last term in Eq. (21) is due to the domain boundary region and is required for conservation of total energy; this term vanishes if the quantum forces of adjacent domains are identical in the boundary region, which is the case for exact quantum force. Spatial domains defined in the CDOFs have been introduced in Ref. 31 but the approximate quantum potential had no explicit dependence on the CDOFs, thus the quantum force had zero component acting on the classical particle. In other words, there was no “backreaction” of QDOFs on CDOFs beyond the influence of coupling in the classical potential. Explicit dependence of the quantum potential on the CDOFs provides such backreaction, which increased accuracy of dynamics in the double-well system of Sec. III C. III. EXAMPLES AND DISCUSSION

First, we consider evolution of a Gaussian wavepacket in the quantum/classical limit of Eq. (18) for a model system of coupled harmonic oscillators introduced in Ref. 32 and its modification by coupling to additional “bath” DOFs. For a Gaussian evolving in a quadratic potential in this limit, the quantum potential exhibits singularities rendering conventional numerical solution of Eq. (18) or QT dynamics unstable. Therefore, the harmonic oscillator problem is examined by solving the time-evolution equations for the parameters of a Gaussian wavepacket described in Sec. III A. Then, we numerically examine dynamics in a double well coupled to the harmonic bath QDOFs, using exact and approximate wavefunction evolution with LQF on domains. A. Dynamics of a Gaussian wavepacket with quantum and classical DOFs

A general multidimensional quadratic potential is given by a real symmetric matrix V2 , V =

1 T z V2 z. 2

(22)

A Gaussian wavefunction is a solution to the TDSE for the potentials that are at most quadratic, such as V of Eq. (22). It can be specified in terms of a time-dependent complex symmetric matrix At , At = A + ıA , two real time-dependent  t and Pt defining a trajectory of the Gaussian cenvectors Q ter, a real overall phase St and a real function γ t normalizing ψ t)  t )T At (z − Q ψ(z, t) = exp(−(z − Q  t ) + ıSt + γt ). +ı PtT (z − Q

(23)

Initial value of γ 0 normalizes ψ(z, 0), γ0 =

1  ln (2/π )Nd det A |t=0 . 4

(24)

Solution of the hybrid quantum/classical Eq. (18) gives the following evolution equations for the parameters: t dQ = GPt , dt

d Pt t, = −V2 Q dt

(25)

Garashchuk, Dell’Angelo, and Rassolov

1.5

(28)

The matrices Gx and Gy are diagonal matrices of the inverse masses for the quantum and classical DOFs, respectively, as in Eq. (1). The total energy of the Gaussian wavefunction is

(b) 1 0.5 1.5 QM QC BO

(c) 1 0.5 0

2 Time [a. u.]

4

FIG. 1. Expectation values of the quantum potential associated with the light particle, initially in the ground vibrational state (a11 = α c ) for the masses of heavy particle equal to m2 = {10, 40, 160} a.u. shown in panels (a), (b), and (c), respectively. The results are obtained with the full quantum (QM, dash) and hybrid quantum/classical (QC, solid) dynamics. The straight line (BO, dot-dash) shows Ux  in the Born-Oppenheimer limit of the infinitely heavy classical particle.

(30)

 t , Pt ) describIn Eq. (30), Ecen is the energy of a trajectory (Q ing the Gaussian center; the next two terms are the kinetic and potential energy due to the wavefunction delocalization; the last term is the expectation value of the quantum potential Ux , defined for the selected DOFs. The hybrid quantum/classical evolution of a Gaussian wavefunction is useful for testing of other propagation methods, and provides a simple model of the quantum DOF interacting with the classical DOF. Note, that for a Gaussian wavefunction the trajectory  t , Pt ) of its center evolves according to the classical equa(Q tions of motion (the quantum force acting on this trajectory is zero), and the expectation values of the position and momentum operators averaged over the wavefunction of Eq. (23) are  t and Pt , respectively. Thus, we analyze the expecequal to Q tation value of the quantum potential Ux , as a measure of the “quantumness” of the system, which for a Gaussian wavefunction is directly related to the dispersion. In the remainder of the paper, the coordinates x and y will be used, respectively, for the quantum and classical DOFs of the system; a vector of coordinates z will represent the quantum bath DOFs. B. Dynamics in a coupled quadratic potential

Let us consider the hybrid quantum/classical evolution according to Eq. (18) for the two-dimensional harmonic oscillator model of Kohen, Stillinger, and Tully.32 The potential, k(x − y)2 Ky 2 + (31) 2 2 describes a light particle bound to a heavy particle moving in a harmonic potential. The force constants are k = 5Eh /a20 and K = 15 Eh /a20 , and the particle masses are m1 = 1 and m2 = {10, 40, 160} a.u. The initial wavefunction is defined as direct product of Gaussian functions,  (32) ψ(x, y, 0) = exp − a11 x 2 − a22 y 2 + γ0 . V (x, y) =

(a)

1 0.5

(27)

Using vectors x and y for positions of the quantum and classical DOFs, as in Sec. II B, the matrices Gqm and Gcl are defined in full dimension Nd     0 0 Gx 0 Gqm = . (29) , Gcl = 0 Gy 0 0

1  E = Ecen + Tr A−1  A GA 2 1  −1 1 + Tr A V2 + Tr(A Gqm ). 8 2



1.5

(26)



dS 1  t ) − Tr(A Gqm ), = PtT GPt − V (Q dt 2 dA 1 ı t = 2At GAt − V2 − 2A Gcl A , dt 2 dγt = Tr(A G). dt

J. Chem. Phys. 141, 234107 (2014)



234107-6

The width parameter for the heavy particle is taken a22 = 44.44 a−2 0 , and for the light particle it is equal to the value of the QDOF ground state, a11 = α c ,  (33) αc = m1 k/2 = 1.118 a0−2 .  t , Pt ) is purely clasSince dynamics of the Gaussian center (Q sical, the initial wavefunction ψ(x, y, 0) can be centered at the origin of coordinates, because the matrix At , describing the wavepacket localization, contains the main features of quantum dynamics. As the classical particle becomes heavier, dynamics of the light particle develops the vibrationally adiabatic character, i.e., contributions from its excited vibrational states become negligible. In other words, as m2 → ∞ the light particle has enough time to adjust to the changing position of the heavy particle, making the Born-Oppenheimer approximation exact. This trend is illustrated in Fig. 1, showing the expectation value of the quantum potential associated with the light particle, Ux , as a function of time for the three values of m2 . For the ground state of a harmonic system, Ux  is equal to one half of the ZPE. Therefore, deviations of Ux  from a constant value correlate with vibrational nonadiabaticity of the light particle, as the wavepacket undergoes “breathing” motion due to interaction with the moving heavy particle. The difference between the fully quantum and the hybrid quantum/classical dynamics decreases with the increase of m2 . Because this two-dimensional model is based on the quadratic classical potential, the hybrid quantum/classical wavefunction exhibits high degree of localization in the CDOF, as trajectories go through nearly focal points at certain moments of time. The wavefunction dispersion σ ,  (34) σ = y 2  − y2 , takes very small values at these instances of time, as shown on a logarithmic scale for σ in Fig. 2(a) (black solid line). At these moments of time the width matrix At exhibits singular behavior: some elements become very large which

234107-7

Garashchuk, Dell’Angelo, and Rassolov

J. Chem. Phys. 141, 234107 (2014) TABLE I. Parameters of the initial wavefunction for the double-well model.  , The wavefunction, ψ(x, y, z, 0), is given by Eq. (23) with the vectors z, Q 0  and P0 replaced with the vectors (x, y, z), (Qx , Qy , Qz ), and (Px , Py , Pz ), respectively. Parameters for all bath DOFs are the ones listed for z.

σ

0.1

Nb=0 Nb=2 Nb=10 Nb=20

0.01

(a)

[Ea]

0

1

2 3 Time [a. u.]

DOF Quantum, x Classical, y Bath, z

4

P [Eh /a0 ]

Mass [a.u.]

a11 = 1.0 a22 = 6.75 a33 = 1.58

Qx = −1.4 Qy = −1.4 Qz = 0.0

Px = 0.0 Py = 10.0 Pz = 7.8

m1 = 1 m2 = 10 mb = 10

C. The double-well system

1. The system and implementation

0.6

0.5 (b) 1

 [a ] Q 0

5

Nb=0 Nb=2 Nb=10 Nb=20

0

A [a−2 0 ]

2 3 Time [a.u.]

4

5

FIG. 2. Dynamics of the quantum/classical system coupled to the quantum bath DOFs. (a) Dispersion in the CDOF y and (b) the quantum potential of the light particle are shown for Nb = {0, 2, 10, 20} bath DOFs. The initial wavefunction is a squeezed Gaussian in x, a11 = 2α c ; the CDOF mass is m2 = 10 a.u.

corresponds to the dispersion σ being small. Then, the solution of Eq. (27) is numerically unstable necessitating very small time steps. (For the grid-based quantum propagation methods, the numerical problems are even more severe in this situation, because of the non-local operator of the kinetic energy.) The singularity issue is alleviated for systems of higher dimensionality but, in general, any numerical solution of the quantum/classical TDSE (18) is prone to difficulties. To illustrate the effect of additional DOFs, we introduced Nb “bath” harmonic oscillators coupled to the CDOF y. The bath DOFs are labeled with coordinates zi , i = [1, Nb ] and treated as quantum DOFs. The bath oscillators are initialized in the ground state. (These DOFs are introduced to explore high-dimensional dynamics and are not meant to describe the temperature of an equilibrated system.) The bath Hamiltonian, added to that of the system, is  Nb  2  (pˆ z,i )2 z k b i + cb zi y . + (35) Hˆ b = 2mb 2 i=1 The parameter values are mb = 10 a.u., kb = 1 a−2 0 , and . The effect of the quantum bath DOFs on the discb = 1 a−2 0 persion in the CDOF and on the quantum potential Ux , associated with the light particle, is illustrated on Fig. 2. The results are obtained for m2 = 10 a.u. and for a squeezed initial wavefunction in the QDOF, a11 = 2α c . The focal points in the CDOF disappear (no vanishingly small values of σ ) as the bath size increases Nb = {0, 2, 10, 20}. The effect of the bath size on the quantum potential of the light particle (Fig. 2(b)) is not as pronounced, since the bath DOFs are not directly coupled to it, but Ux  becomes a smoother function of time for larger Nb .

The second model system is relevant to reactive processes in condensed phase; here, we focus on the energy flow between quantum and classical DOFs, not on the QM tunneling. The potential consists of a double well in the CDOF to which a harmonic QDOF is linearly coupled. The functional form of the potential is V (x, y) = y 2 (ay 2 − b) +

b2 c(x − y)2 + , 2 4a

(36)

−2 and the parameter values are: a = 1 Eh a−4 0 , c = 4 Eh a0 , and −2 b = 4 Eh a0 . The particle masses are set to m1 = 1 a.u. in the QDOF x, and to m2 = 10 a.u. in the CDOF y, as in the model of Sec. III B. We examine the reactant-well, or survival, probability of a wavefunction, initially localized in the reactant (left) well, evolving in the hybrid quantum/classical dynamics framework according to Eq. (18). The initial wavefunction is a product of one-dimensional Gaussians, centered at the reactant minimum of the potential. The parameters are listed in Table I. The QDOF wavefunction describes the ground state in the absence of couplings. The momentum for the CDOF y is chosen so that the corresponding kinetic energy exceeds the barrier height by 25%. Thus, the heavy particle is treated classically, but is fully coupled to the light quantum particle. The energy exchange between the two particles may influence reactivity even for the reaction path along the CDOF. The effect of energy flow between the quantum and classical particles is further explored by adding the harmonic bath QDOFs coupled to the CDOF. The initial wavefunction for these bath DOFs are taken as the ground states of the respective DOFs, assuming zero coupling; the initial momenta are specified in Sec. III C 3. The initial positions of the QTs representing these multidimensional Gaussian wavefunctions, are chosen randomly according to the normal deviates.33 To assess the accuracy of the hybrid quantum/classical QT-based dynamics for this two-dimensional system, we performed the wavefunction evolution according to Eq. (18) using the split-operator/Fast Fourier Transform method14, 34 on a two-dimensional equidistant spatial grid. This exact numerical implementation yielded unstable dynamics because of the singularities in the quantum potential term, which is numerically ill-defined when a wavefunction exhibits high degree of localization in the CDOF as discussed in Sec. III B. The most accurate procedure of incorporating the quantum potential term in Eq. (18) that we came up with, was (i) to find the second derivative of |ψ| with respect to the CDOF in

Garashchuk, Dell’Angelo, and Rassolov

J. Chem. Phys. 141, 234107 (2014)

coordinate space by the fifth order finite difference, (ii) to add a small positive constant δ to |ψ| in the denominator and multiply by ψ, and (iii) to add the resulting function, Uc , as a non-homogeneous term into the kinetic energy step (37)

δ = 10−8 is a representative value of this parameter. (Other ways that were explored and found less stable include differentiation of |ψ| in the momentum space and determination of the phase of ψ to analytically cancel small amplitudes, |ψ|, in the numerator and denominator.) With the exact implementation according to Eq. (37) the reactant-well probabilities were converged for short times, t < 4 a.u. These results were used to assess the accuracy of the QT implementation with the approximate quantum correction defined on spatial domains in the CDOF, described in Sec. II C. The QT dynamics was then extended to include the harmonic bath as was done in the harmonic oscillator model. To implement the LQF on domains in the CDOF, the spatial domains are defined by a smoothed step-function in y, 1 (tanh(κ(y − Yl + dY )) − tanh(κ(y − Yl ))). (38) 2 The domain borders are specified by a set of values {Yl }, l = 0, 1, . . . , L, evenly spaced by dY ≡ Yl − Yl − 1 . At time t the main contribution to the lth domain comes from trajectories with current positions (xt , yt ) such that yt ∈ [Yl − 1 , Yl ]. Parameter κ controls the steepness of the domain walls; the domain function approaches the step function as κ → ∞. Smaller values of κ result in a more delocalized domains, which (at the cost of the LQF accuracy) improves convergence with respect to the number of trajectories. In the calculations here, the value of κ is chosen so that the overlap between the nextto-nearest neighbor domains is negligible, and the adjacent domain functions add up to 1 in the region of their overlap. Using l =

η ≡ κdY,

η  1,

(39)

the overlap of the adjacent domains and next-nearest neighbor domains is  ∞ −1

∞ 1/2 , (40) l l+1 dy l l dy ≈ η −1 −∞ −∞



∞

−∞

l l+2 dy

−1

∞

−∞

l l dy

≈

ηe−2η . η−1

A B C D F

L

κ [a−1 0 ]

Ntraj [103 ]

Nb

cb [a−2 0 ]

10 10 [1, 15] 10 10

25 [15, 75] 50 50 50

[5, 80] 80 80 40 40

0 0 0 {4, 8, 12} 4

... ... ... 0.1 {0.1, 0.2, 0.4}

interval y = [ − 2, 2] a0 . The outer borders are placed at infinity, Y0 = −∞ and YL = ∞. The values of κ are considered in the interval κ = [15, 75] a−1 0 . The remaining parameters of calculations are summarized in Table II. First, we performed test calculations to investigate the energy conservation in the course of dynamics with 10 domains, L = 10, and convergence with respect to the number of trajectories, Ntraj . The width of the step-function for 10 domains is dY = 0.4 a0 . Parameters relevant to this calculation are given in line A of Table II. The total energy was conserved within 0.02%, while the probability (at the second maximum t = 10 a.u.) is converged within 0.4%. Fig. 3(a) shows P(t), including the sampling error, for the number of trajectories, Ntraj = {5, 10, 20, 40, 80} × 103 . The average absolute deviation between the largest and smallest trajectory ensembles is below 1% (the maximum deviation is 2%), while the average statistical error is 0.5% and 0.1%, respectively. We do not achieve an ideal uniform convergence, because at any time 0.5

The convergence studies are performed for the bare twodimensional double-well. The survival or reactant-well probability is defined according to  w (k) , (42) P (t) = yt(k)