Article pubs.acs.org/JPCA

Quantum Dynamics with Gaussian Bases Defined by the Quantum Trajectories Bing Gu and Sophya Garashchuk* Department of Chemistry & Biochemistry, University of South Carolina, Columbia, South Carolina 29208, United States ABSTRACT: Development of a general approach to construction of efficient high-dimensional bases is an outstanding challenge in quantum dynamics describing large amplitude motion of molecules and fragments. A number of approaches, proposed over the years, utilize Gaussian bases whose parameters are somehowusually by propagating classical trajectories or by solving coupled variational equationstailored to the shape of a wave function evolving in time. In this paper we define the time-dependent Gaussian bases through an ensemble of quantum or Bohmian trajectories, known to provide a very compact representation of a wave function due to conservation of the probability density associated with each trajectory. Though the exact numerical implementation of the quantum trajectory dynamics itself is, generally, impractical, the quantum trajectories can be obtained from the wave function expanded in a basis. The resulting trajectories are used to guide compact Gaussian bases, as illustrated on several model problems. general wave functions with the system size.12 The traditional exact approaches based on direct-product bases or discrete variable representation of nuclear wave functions are practical for systems of 4−5 atoms (about 12 DOFs) due to exponential scaling of the basis size with dimensionality of a system. Over the years many efforts went into the development of more efficient adaptive bases for exact time-evolution methods and of approximate and semiclassical dynamics methods. One of the most developed exact method balancing computational cost and accuracy is the multiconfiguration time-dependent Hartree method (MCTDH).13−15 The MCTDH, based on contraction of a general basis to single (or a few) particle functions, and its extension (the multilayer MCTDH16) have been very useful and successful in studies of bound high-dimensional molecular systems with PES expressed in a “natural” form related to the basis used. A closer to this work extension of the MCTDH utilizes the variationally determined time-dependent Gaussian bases and is termed the multiconfiguration Gaussian (vMCG) approach.17−19 The Gaussian basis functions are widely used in quantum dynamics, as they often offer an intuitive “classical” interpretation of the bases, but the use of Gaussian bases has been explored in the time-independent context as well. The idea goes back to Davis and Heller,20 who defined the basis functions in the classical phase space. One of the later developments is the construction of compact Gaussian bases by adapting the parameters of the functions to the external

1. INTRODUCTION Importance of the quantum-mechanical effects associated with the nuclei is gaining recognition in chemistry and physics, as researchers manipulate matter, light, and electric and magnetic fields at the atomistic level to develop advanced materials and molecular structures with desired properties. A few recent examples of systems with significant nuclear quantum effects (NQE) are the quantum sieving, which is the hydrogen/ deuterium isotope separation, based on the zero-point-energy difference,1−4 the dependence of optoelectronic properties of P3HT/PCBM heterojunction on the H/D substitution of the polymer P3HT,5 the proton conductance in low-dimensional boron nitrides and silicon-based structures.6 NQE are expected to be the most pronounced for light nuclei at low temperatures when a typical energy of a process is comparable to the separation between rovibrational energy levels for the nuclei. To describe the quantum process unfolding in time rigorously, one has to solve the timedependent Schrödinger equation (TDSE) for the nuclei evolving on a single or multiple electronic potential energy surfaces (PES): ∂ Ĥ ψ (x,t ) = iℏ ψ (x,t ) ∂t

2

Ĥ = K̂ + V̂ (x)

ℏΔ K̂ = − 2m (1)

K̂ and V̂ are the kinetic energy and potential energy operators, for simplicity written for a single PES and for a particle mass m in all degrees of freedom (DOFs); Δ is the Laplacian operator. Although the propagation schemes, such as the split-operator and Chebyshev expansion of the Hamiltonian operator, are fairly well-established,7−11 the numerical representation of a wave function remains a a bottleneck to the quantum studies of high-dimensional systems, due to the scaling properties of © XXXX American Chemical Society

Special Issue: Ronnie Kosloff Festschrift Received: October 13, 2015 Revised: December 23, 2015

A

DOI: 10.1021/acs.jpca.5b10029 J. Phys. Chem. A XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry A potential according to a semiclassical criterion.21,22 Interestingly, it has been shown that a Gaussian basis with strongly overlapping functions gives better accuracy for the eigenstates. Poirier has developed a method of constructing efficient bases using wavelets, obtained by canonical orthogonalization of the coherent state Gaussians.23 Together with the truncation scheme, the method is highly efficient for the rovibrational spectra calculations and it has been further improved by introducing the momentum-symmetrized Gaussian basis.24,25 Impressive results were obtained for model systems of up to 27 dimensions and several polyatomic molecules, such as CH3CN.26 A related approach is based on the von Neumann lattice of phase space Gaussians and incorporates periodic boundary condition generating bases with superior convergence properties and efficient in the classical limit.27 In the time-dependent domain the idea of using localized functions of Gaussian form, or Gaussian wavepackets (GWP), to represent nuclear wave functions evolving in time became popular thanks to Heller.28−31 On the one hand, a GWP exactly describes evolution of an initially Gaussian wave function in a time-dependent parabolic potential, such as the local harmonic approximation to a general potential around a trajectory position, which is a highly useful model of molecular vibrations. On the other hand, the GWP dynamics is easily connected to classical mechanics, which is adequate to describe nuclear motion in many situations and is widely used to simulate molecular systems: in a parabolic potential the center of a quantum GWP moves classically. Finally, Gaussian functions are convenient to work with in numerical implementation, because many integrals can be performed analytically. Though using a single Gaussian to approximate a wave function has an obvious limitation (a Gaussian is a localized wave function and thus cannot describe a process involving bifurcation of the probability density), the idea of using multiple Gaussians to represent a wave function is used in numerous exact and approximate methods. Some of the exact methods are the coupled coherent Gaussians (CCG)32,33 with the GWP moving along the classical or multiple Ehrenfest trajectories, the matching-pursuit/split-operator-Fourier-transform (MP-SOFT)34,35 with GWPs found through a reexpansion of a wave function, ab initio multiple spawning (AIMS)36 for nonadiabatic dynamics, based on classically evolving position and momentum of the Gaussian centers, and the vMCG MCTDH approach,17−19 where Gaussian parameters are determined variationally. The last has the advantage of rigorously conserving the wave function energy, but the resulting coupled equations are tedious to implement numerically and the Gaussian parameters may lose their physical meaning, compared to the GWPs centered on classical trajectories. The more recent Gaussian-based methods are the basis expansion leaping multiconfigurational Gaussians (BELMCG) of Frankcombe37 and the trajectory-guided timeindependent Gaussian basis of Saller and Habershon.38 In the former method (BEL-MCG), as in MP-SOFT, the Gaussian basis changes in time in a discrete manner when the reexpansion of the wave function in terms of new set of Gaussians is performed with the goal of keeping the basis size small. In the latter method, the time-independent Gaussian basis is constructed from the functions centered at the select positions and momenta of the classically evolved trajectories, sampling the Wigner transform of the initial wave function in phase space.

In the area of semiclassical Gaussian-based dynamics the thawed and frozen Gaussians28,29 and the initial value representation method of Herman and Kluk39,40 are the most popular. One important conclusion emerging from the semiclassical work is, although a single thawed Gaussian, i.e., a Gaussian with the complex time-dependent width parameter, is the exact general solution of the quantum harmonic oscillator model, use of multiple thawed Gaussians in a basis representation of a wave function is numerically tedious and unstable. Thus, in our approach described below, we limit ourselves to the “frozen” Gaussian basis of functions with fixedin-time width. We note, at this point, that in the area of timeindependent rovibrational calculations, it has been demonstrated that adaptation of the Gaussian width parameters to the potential increases efficiency of a basis.21,41 Therefore, in the context of time-dependent bases, a middle ground between frozen and thawed Gaussian basis functions merits exploration. The remainder of this paper is organized as follows. In section 2 we describe the formalism, including the key features of the quantum trajectory dynamics and their incorporation into the frozen Gaussian bases. Section 3 gives details of the numerical implementation and results for standard in chemistry one- and two-dimensional model systems. Discussion and summary are presented in section 4.

2. FORMALISM In this section the basic features of time-dependent bases and of the QT dynamics are reviewed first. Then, they are combined to construct the QT-guided Gaussian basis and the working equations of the resulting dynamics are derived. 2.1. Schrö dinger Equation in the Time-Dependent Basis. To solve the TDSE in a time-dependent and, in general, nonorthogonal basis, one starts by representing a wave function in terms of the basis functions gk(x,t), where k enumerates the basis functions, k = 1, 2, ..., Nb. Nb is the number of the basis function. For simplicity we assume here that Nb is constant in time and that x is the vector of coordinates describing a particle of mass m in the f-dimensional Cartesian space. The wave function at time t, projected onto the subspace spanned by the basis functions is Nb

ψ (x,t ) =

∑ ck(t ) gk(x,t ) k=1

(2)

where {ck(t)} are the expansion coefficient. Let us assume that each basis function has no explicit time dependence but instead depends on the time-dependent parameters z = {z1, ..., zNp}, completely specifying the basis function, gk (x,t ) ≡ gk (x,zk1(t ),...,zkNp(t ))

(3)

Np is the number of parameters taken to be the same for all basis functions. Denoting ż ≡ dz/dt, the time derivative of a basis function is dgk dt

Np

=

∑ zk̇ μ μ=1

∂gk ∂zkμ

(4)

Substitution of eqs 2 and 4 into the TDSE (eq 1) and multiplication by ⟨gn| give the following matrix equation: iℏM B

dc = (H − iℏD)c dt

(5) DOI: 10.1021/acs.jpca.5b10029 J. Phys. Chem. A XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry A where Np

Mkn = ⟨gk |gn⟩

Dkn = ⟨gk |

∑ zṅ μ| μ=1

∂gn ∂znμ

⟩ (6)

and H is the Hamiltonian matrix, Hkn = Kkn + Vkn

p·p dS(x,t ) = − V (x) − U (x,t ) dt 2m

(12)

∇·p dρ(x,t ) =− ρ(x,t ) m dt

(13)

As seen from eq 12, all quantum effects in dynamics of ψ(x,t) are formally generated by the quantum potential U,

(7) f

ℏ2 ∂2 Kkn = − ⟨gk | ∑ 2 |gn⟩ 2m ∂xj j=1

U (x,t ) = −

Vkn = ⟨gk |V (x)|gn⟩

gk =

p dx = dt m

⎛ det A ⎜− 1 exp f ⎜ 2 π ⎝

⎞

∑ (xμ − qkμ(t )) Aμν(xν − qkν(t ))⎟⎟ μν

⎠

(15)

(9)

(10)

Substituting eq 9 into TDSE (eq 1), one obtains coupled equations determining the amplitude and phase, which after identification of the trajectory momentum with the gradient of the phase, p = ∇S

4

The integer f refers to the number of DOFs; the DOF indices, labeled by the Greek letters μ and ν, change from 1 to f. The basis function gk is centered at the position of the kth quantum trajectory, qk, their number defining the size of the basis Nb. For simplicity, assume the matrix Aμν is diagonal, and its nonzero elements are Aμμ = αμ. Note that the form of the Gaussian in eq 15 above, is simpler than the standard frozen GWP related to a general Gaussian, solving the quantum harmonic oscillator.48 The latter includes the momentum pk(t) of the trajectory specifying the GWP center. This momentum defines a linear in the x phase, ı S lin = pk (t )(x − q k(t )) + γk(t ) ℏ

The dynamics equations are conveniently formulated in terms of the probability density ρ(x,t), ρ(x,t ) = ψ *(x,t ) ψ (x,t ) = A2 (x,t )

(14)

entering the formulation on par with the external classical potential V. The potential U(x,t) is a nonlocal time-dependent function dependent on the evolving wave function, described by an ensemble of QTs. As follows from eq 13, the probability density within the volume element of each trajectory is conserved: ρ(x,t)δxt = ρ(x,0) δx0.47 Therefore, if the exact QTs could be computed, then a wave function will be well-represented by the same ensemble of QTs, which sampled subspace of nonzero density ρ(x,0) at initial time, throughout the entire propagation. An ensemble of QTs can be interpreted as an adapting grid optimally suited to describe ψ(x,t). Implementation of eq 12 and following from it Newton’s equations of motion defining a QT, requires computation of the quantum potential and its gradient. In general, U is singular when ψ takes zero values, which makes the QT formulation impractical as an exact method for quantum systems exhibiting interference. In this work, we use the concept of QTs to guide a Gaussian basis without computing the quantum potential itself: the QT momentum, p, is determined from the wave function represented in a basis, and its value is used to increment positions of the QTs according to eq 11. In other words, the guiding trajectories are reconstructed from ψ(x,t) expressed by eq 2 without solving Newton’s equations of motion. 2.3. Implementation with Real Frozen Gaussians. We construct an adaptive Gaussian basis out of time-dependent multidimensional Gaussian basis functions of the following form

(8)

Equation 5 specifies evolution of the expansion coefficients and, thus, the evolving-in-time wave function in a basis for any time dependence of the basis function parameters. The choice of this time dependence is crucial for accuracy and conservation properties of the ensuing dynamics. The normalization of the wave function in a finite basis is conserved regardless of its time dependence.42 The energy is conserved for (i) any time-independent basis or for (ii) a time-dependent basis with variationally determined parameters (such as in vMCG), and also for (iii) a complete (on the space relevant to the problem) basis with arbitrary time dependence. For a finite time-dependent basis the energy conservation may serve as an indicator of its completeness as time progresses, i.e., of the dynamics accuracy. Though the variationally determined compact basis is conceptually appealing, in practice it is known that the equations defining its parameters are stiff and the numerical solutions are unstable.17,43 Thus, using positions and momenta of classical trajectories, often sampling the phase space of an initial wave function, or Ehrenfest-type trajectories, is a popular choice for defining a Gaussian basis. The resulting dynamics equations are much easier to implement, though there is a possibility of missing regions of space inaccessible to classical trajectories yet involved in quantum dynamics. In the remainder of this section we introduce the quantum trajectories (QT) and, capitalizing on their inherent property of following the evolution of the probability density, use them as guiding trajectories to define compact Gaussian bases (GB) referred to as QTGB. 2.2. Bohmian Mechanics. The Madelung−de Broglie− Bohm, also referred to as the hydrodynamic or the QT, formulation of the TDSE44−46 is based on the polar representation of a complex wave function in terms its amplitude A(x,t) and phase S(x,t), which are both real functions of x and t, ⎛i ⎞ ψ (x,t ) = A(x,t ) exp⎜ S(x,t )⎟ ⎝ℏ ⎠

ℏ2 ∇2 A(x,t ) 2m A(x,t )

The time-dependent function γk(t) is associated with the classical action function, and because the basis functions contribute to ψ(x,t) with coefficients ck(t), this part of the phase can be included in ck(t). The choice of the basis functions in eq 15 to be real, i.e., without pk(t) as an explicit parameter, means

(11)

yield the quantum Hamilton−Jacobi equation for S and the continuity equation for the probability density: C

DOI: 10.1021/acs.jpca.5b10029 J. Phys. Chem. A XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry A that there is no term dp/dt in eq 4 and no need for tedious and unstable computation of ∇U. The matrix elements of eqs 6 and 8 are expressed as follows. The arguments of functions are omitted for clarity. The overlap matrix elements and the timederivative terms are f

Mkn =

⎛ αμ ⎞ (qkμ − qnμ)2 ⎟ ⎠ 4

∏ exp⎜⎝− μ=1

Dkn = ⟨gk | − i

∂ |g ⟩ = ∂t n

f

∑i μ=1

αμ 2m

pnμ (qnμ − qkμ)Mkn

follow the description of several practical aspects of the QTGB implementation. The atomic units are used henceforth. 3.1. Practical Details. In a numerical implementation of the Gaussian bases one often has to deal with an ill-conditioned overlap matrix, which happens when two GWPs strongly overlap or if a basis becomes overdetermined.20 Formally, the QTs never cross, which should alleviate this problem: as two QTs approach each other, the exact quantum force becomes very large and “repels” these trajectories. In reality, however, this behavior makes the QT dynamics inherently unstable. The QTs reconstructed from the exact wave function according to eq 20 would exhibit the noncrossing behavior in the exact limit, but in practice, the Gaussian functions may come close enough to generate an ill-conditioned overlap matrix M, which has to be inverted to implement eq 5. To avoid this problem, following the work of others,34,37,38 ψ(x,t) is occasionally reexpanded in a new set of Gaussian basis functions. After a certain time interval T, we resample ψ(x,t) in terms of new QTs. Besides preventing the “collision” of two GWPs, the reexpansion procedure removes basis functions where the probability is negligible and adds more basis functions where the wave function becomes delocalized in space. The exact QTs track such changes in the wave function density through the probability density conservation. However, because our basis functions are of fixed width, they do not describe these changes within a basis function itself−these changes are captured by the expansion coefficients. Therefore, because our goal is a compact representation of ψ(x,t) in terms of small number of GWPs, not in terms of a complete QT ensemble, instantaneous adjustments of the number of basis functions and of their width outside the time-propagation loop are useful. In the examples below we have used a simple reexpansion algorithm: start with an outlying initial point q1 and compute the probability density, |ψ(q1,t)|2 = |∑kckgk(q1)|2. If this value is below a predefined threshold ϵ, |ψ(q1,t)|2 < ϵ, disregard the point. Otherwise, take this position as the center of a GWP that will be included in the new basis set. Move to the right by an interval Δq, and repeat the last step until we get to the region where the density is negligible. In high dimensionality, a new set {q′k} can be generated through the random importance sampling. The basis size Nb′ resulting from this procedure can be different from Nb of the previous propagation. The advantage is that there is a flexibility to change the number of basis functions depending on the wave function localization. The new expansions coefficients, {ck′}, are found by minimizing the error between the expansions in the new and old bases,

(16)

(17)

The matrix elements for the kinetic energy operator are f

Kkn =

∑ μ=1

αμ ⎛ αμ ⎞ ⎜1 − (q − qnμ)2 ⎟Mkn ⎠ 4m ⎝ 2 kμ

(18)

The potential energy matrix elements are evaluated approximating V with the second-order expansion near the Gaussian center qk or qn, known as the local harmonic approximation: ⎛ Vkn = ⎜⎜V (q̅) + ⎝

f

∑ μ=1

⎞ 1 ∂ 2V (q̅) ⎟ Mkn 4αμ ∂xμ 2 ⎟⎠

(19)

where q̅ defines the midpoint of two Gaussians, q̅ = (qk + qn)/ 2. The QT position is incremented at each time step according to its momentum determined from ψ(x,t): p dq = dt m

⎛ ∇ψ ⎞ p = ∇ S = ℑ⎜ ⎟ ⎝ ψ ⎠

(20)

The local harmonic approximation to V is frequently made in practical computation of the potential energy matrix elements of the dynamics approaches utilizing basis functions localized in coordinate space. In the CCG method, for example, the potential matrix elements are computed within the linear local expansion of V(q)̅ (no second derivative of V). The MCTDH relies on the expansion of a PES in a product basis and there are dedicated programs to fit a PES in this form.49−51 The quadratic expansion of V (the local harmonic approximation) will be used in applications of section 3. The QTGB simulation begins with an expansion of an initial wave function in terms of the Gaussian functions of fixed width α whose centers are the initial positions of the QTs. The number of the QTs depends on the desired accuracy of the expansion. For high-dimensional systems the sampling of the initial QT positions will be random or quasi-random. The width parameter α affects the basis size (the number of QTs) and the accuracy of the quadratic approximation in evaluation of the potential energy matrix of eq 19. The value of α has to be chosen with some prior knowledge about the surface and with several trials. The initial momenta of QTs are defined by the initial wave function. Equations 5 and 16−20 are solved to evolve a wave function in time.

I = ∥∑ ck′gk′ − ψ (x,t )∥2

(21)

k

The resulting matrix equation on {c′k} is (22)

M′c′ = b′

where ′ = ⟨gk′|gn′⟩ Mkn

3. IMPLEMENTATION AND EXAMPLES The QTGB has been applied to several model systems, including the Morse oscillator representing the vibration of H2, scattering on the Eckart barrier, and the double-well systems in one and two dimensions. The results, compared to the exact quantum-mechanical results obtained with the split-operator propagator and the grid representation of wave functions,52,53

bk′ = ⟨gk′|ψ ⟩

(23)

The time evolution of the expansion coefficients is performed with the second-order time differencing scheme,52 c(t +Δt ) = c(t −Δt ) + 2Δt

dc dt

(24)

with dc/dt defined by eq 5. D

DOI: 10.1021/acs.jpca.5b10029 J. Phys. Chem. A XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry A To complete the propagation set up in the new basis, the QT momenta are defined by differentiating the phase of ψ(x,T) and evaluating it at x = qk′ according to eq 20. Because, the main contribution to ψ(x,t) at a QT position comes from the GWP centered at that position, the right-hand-side of eq 20 is nonsingular. Moreover, to define the trajectory momentum, we use a convoluted wave function, as outlined below in one dimension: ⎛ β ⎞1/2 ψβ̃ (x ,t ) = ⎜ ⎟ ⎝ 2π ⎠

Table 1. Propagation and Wavefunction Parameters for the One-Dimensional Systems for a Particle of Mass ma model

Nb

α

β

Morse oscillator double well Eckart barrier

10

64

8

0.2

9.16

8 15.4

2 1

0.001 0.001

0.5 1.0

16b 17c

Δt

α0

q0

p0

1.3

0

−2.5 −5.0

0 5.0

m 925 1 1

a

∞

∫−∞

e

−β /2(x − y)2

ψ (y ,t ) dy

All quantities are given in atomic units. The parameters of the initial wavefunction α0, q0, and p0 are defined in eq 28. Nb is the number of the basis functions given eq 15; α specifies their width, and β defines the convolution. bNb increases up to 23 for t = [0, 3.0] au. cNb increases up to 40 for t = [0, 1.6] au.

(25)

In the limit of infinite β, ψ̃ β → ψ. The convolution is done to ensure smoothness of the momentum as a function of position for a sparse set of QTs. Using the convoluted wave function ψ̃ β of eq 25 in eq 20, the momentum is defined as follows: ⎛ β ⎞1/4 ⎛ β p = − ⎜ ⎟ ℑ⎜ ⎝π ⎠ ⎝ ψ (x ,t )

∞

2

⎞

∫−∞ (x − y)e−β/2(x−y) ψ (y ,t ) dy⎟⎠ (26)

With ψ(y,t) written as a sum of Gaussian functions, the integrals above are analytical. When β → ∞, the convolution has no effect on the obtained momentum; if β → 0, then all the momenta are equal to the average value of p and deviate significantly from the QT values. We use large values of β to smooth out the momentum as a function of the GWP centers. This value depends on the width of the Gaussian functions and should be such that contributions in the convolution procedure come from several GWPs. Note that because {pk} are not part of the basis function definition, the convolution procedure or other modifications of the momenta do not affect the accuracy of the reexpansion. They will affect, however, the completeness of the basis at later times: having momenta closer to the QT values yields an accurate basis representation of the wave function for longer times. 3.2. Morse Oscillator. The first model is the Morse oscillator representing the vibration of a nonrotating H2 molecule, described by an anharmonic single-well potential, V (x) = D(1 − e−a(x − x0))2

Figure 1. Morse oscillator: the imaginary part of the autocorrelation function given by eq 29. The results of the QT-guided Gaussian basis (Nb = 10) and the exact quantum (QM) dynamics are shown as the solid line and circles, respectively. The result of propagation in the quadratic approximation to the potential (HO) is shown with the dashed line.

harmonic expansion of the potential near x = x0, also shown in Figure 1, illustrates the anharmonicity of the system. 3.3. Scattering on a Barrier. Scattering on a barrier, essential for reactive dynamics of open systems, unfolds over vast regions of coordinate space and is a challenging problem for the trajectory based methods due to QM tunneling at low energies. We test the QTGB performance on the Eckart barrier, a widely used one-dimensional model mimicking the transition state of the H + H2 reaction. The potential energy has the following form:

(27)

The parameter values are D = 0.176 Eh, a = 1.02 a0−1, x0 = 1.4 a0. The particle mass is the reduced mass of the diatomic, m = 925 au. The initial wave function is a displaced Gaussian wavepacket, ψ (x ,0) =

4

2α0 exp( −α0(x − q0)2 + ip0 (x − q0)) π

V (x) = D cosh−2(ax)

The parameter values, after rescaling the particle mass to m = 1 au, are D = 16 Eh and a = 1 a0−1. The initial wave function is a Gaussian wavepacket given by eq 28, whose parameters are listed in Table 1. The barrier is steep and narrow, which makes it particularly challenging for the frozen Gaussian basis at low energies, because the transmitted component of a wave function becomes diffuse whereas the interference of the incoming and reflected components results in localized features of the probability density in the barrier region. Thus, within the restrictions of the frozen Gaussians and the local harmonic approximation to the potential, we have to work with highly localized basis functions. The basis width parameter α equals 15.4 a0−2, whereas the wave function width parameter is α0 = 1.0 a0−2. To capture the main features of the dynamics, we had to reexpand the wave function more frequently than in other models, i.e., every 50 time steps (Δt = 0.001 au), and to

(28)

whose parameters along with the propagation parameters are listed in Table 1. The propagation is accomplished using Nb = 10 functions in the QTGB with the width parameter α = 64 a0−2. To assess the accuracy of the wave function over time, the autocorrelation function C(t) is computed, using the property fulfilled for the real initial ψ0, ψ(x,−t/2) = ψ†(x,t/2), C(t ) = ⟨ψ0|e−iHt |ψ0⟩ =

∞

∫−∞ ψ 2(x ,t /2) dx

(30)

(29)

The imaginary part of C(t) obtained from the QTGB and the exact QM dynamics are shown in Figure 1. The agreement is very good and we did not encounter situations of the overlap matrix M becoming ill-conditioned. Thus, no basis reexpansions have been made. 0(C(t )) for the same ψ0 evolving in the E

DOI: 10.1021/acs.jpca.5b10029 J. Phys. Chem. A XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry A increase the number of basis functions from Nb = 17 to Nb = 40 to reach the final time of t = 1.6 au. We also had to use a wider convolution (β = 1.0 a0−2) compared to that used for the bound problems. The accuracy of reexpansion was I < 5 × 10−5, where I is the residual error between the wave functions represented in the new and old bases, given by eq 21. The time evolution was performed using relatively small number of basis functions. Figure 2 shows the wave function at t = 1.2 au and

parameters are listed in Table 1, is a Gaussian (eq 28) localized in the left well. The initial energy for this wave function is E0 = 0.9318 Eh, which is about 68% of the barrier height. The feature of this model is the wide barrier width (the distance between the two minima is around 4.7 a0), which causes difficulty for semiclassical methods.34 The QTGB dynamics is accomplished with Nb = 16 basis functions of the width α = 16 a0−2. The propagation details are given in Table 1 and the wave function reexpansion has been performed every 1500 time steps. Figure 3 shows the wave

Figure 2. Wave function scattering on the Eckart barrier obtained with QTGB (red dashed line) and exact QM (black solid line) dynamics. The blue dot-dashed line outlines the rescaled potential.

Figure 3. Probability density, |ψ(x,t)|2, for the one-dimensional double well obtained with the QTGB and exact quantum (QM) dynamics at t = 3.0 au shown with the solid line and circles, respectively. The dashed line shows the initial density, and the dot-dashed line shows the shape of the potential given by eq 31.

1.6 au The main features of dynamics are reproduced, but there are unphysical oscillations due to sparsity of localized basis functions in the transmitted region, which will deteriorate the solution of the TDSE at long times. This obvious deficiency of the frozen Gaussian basis has been encountered before, for example, in scattering on a smoothed square barrier simulated with the BEL-MCG.37 To increase the accuracy, one can (i) increase the basis size: for example, 200 GWPs are used in ref 37. A more appealing option is (ii) to relax the frozen Gaussian restriction without complicating the equations of motion of the basis function centers. The latter may be accomplished by having GWP width parameter depend on the function’s position, i.e., α ≡ α(x), and adapting α to the wave function shape during the reexpansion, while keeping α frozen in time between the reexpansions to maintain the simplicity of equations for the trajectory motion. 3.4. One-Dimensional Double Well. The double-well potential often represents a reaction proceeding in condensed phase; for the light particles such as a proton or an electron, the quantum-mechanical tunneling may be significant at low energies. In one-dimension the classical trajectories with energies below the barrier top will not cross the barrier, while some of the QTs, following the probability density flow, gain energy from the QT ensemble and cross the barrier even if their initial energy was too low for that. The total energy of a QT ensemble is conserved, but it can be redistributed between the QTs via the quantum potential. Thus, in this system, the QT-guided Gaussian basis can be initially localized in the reactant well, yet describe the wave function density in the product well at some later time. For illustration we take the double-well potential simulating an electron transfer:34 1 4 1 2 V (x ) = x − x 16η 2 (31)

function at t = 3 au in the process of tunneling, which is in good agreement with the exact quantum results. The GWP centers {qk} or the QT positions, plotted in Figure 4a as functions of time up to t = 6.0 au, illustrate the adaptation of the basis and the effect of reexpansions. Note that around t = 5.4 au one of the QTs goes around a node in ψ indicative of interference. At the node, located near x = 1.5 a0 as marked in Figure 4a, the wave function changes sign, which generally leads to singularities in the quantum potential and a breakdown of the QT dynamics. In the QTGB dynamics, however, the interference pattern is reproduced through superposition of the basis functions. The propagation error increases with time but, overall, the dynamics is fairly stable. The absolute value of the autocorrelation function for t = [0, 24] au, shown in Figure 4b, was obtained using eq 29 from dynamics performed up to t = 12 au. Longer (up to t = 24 au) time evolution of comparable quality required twice shorter time steps with reexpansions performed every 1000 time steps. The residual error of the wave function reexpansion given by eq 21 is less than 5 × 10−5. The representative dependence of the wave function energy on time is shown in the inset of Figure 4b. Apart from the numerical accuracy the propagation the energy conservation indicates the completeness of a time-dependent basis; the energy, shown with respect to the bottom of the well, is conserved within 0.7% of its initial value. 3.5. Two-Dimensional Double Well. Finally, we examine the dynamics of a two-dimensional system with the potential consisting of a double well linearly coupled to a harmonic oscillator. The potential is taken from ref 54,

where η = 1.3544 Eh−1 a04. The barrier of height Vb = 1.3544 Eh is centered at x = 0. The initial wave function, whose

V (x ,y) = y 2 (ay 2 − b) + F

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

(32)

DOI: 10.1021/acs.jpca.5b10029 J. Phys. Chem. A XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry A Table 2. Two-Dimensional Double-Well Modela Wave Function and Propagation Parameters α0 0.5

q0

p0

m

−1.4 0.0 1 Basis Parameters and Frequencies

Δt 0.001

QT-guided Gaussian basis exact QM Nb α νg νe

10 × 10 16 0.4827 0.7110

12 × 12 16 0.4822 0.7180

16 × 16 32 0.4830 0.7209

0.4829 0.7163

a

The initial wavefunction and propagation parameters for three QTGB calculations, and the eigen-frequencies of the ground and first excited states are given in atomic units. The parameters α0, q0, and p0 are defined in eq 28. Nb and α are the number and width of the basis functions given eq 15.

perspective, the exact QM dynamics was performed on a 512 × 512 equidistant spatial grid. The autocorrelation functions are computed up to t = 6.0 au and transformed into the energy domain using harmonic inversion55,56 to increase the resolution of the eigen-energies. The positions of the Gaussian centers at t = 3.0 au, plotted in Figure 5, illustrate the adaptation of the QTGB functions that were started as the square “grid” of positions. The frequencies of the symmetric ground and first excited states obtained from the QTGB dynamics are listed in Table 2 as νn ≡ En/(2π). The simulation is performed with several sets of parameters listed in Table 2. The frequencies, which are very sensitive to the quality of the correlation function, are in good agreement with the QM results and show a small dependence on the chosen bases. The basis reexpansion was not needed to obtain the short-time dynamics from which the eigen-frequencies were extracted.

Figure 4. One-dimensional double well. (a) Positions of the basis function centers, i.e., QTs, as a function of time. The wave function reexpansion has been performed at t = {1.5, 3.0, 4.5} au. One of the QTs goes around the wave function node (ψ(x,t) = 0) at t = 5.4 au. (b) Absolute value of the autocorrelation function obtained according eq 29 from the dynamics up to t = 12 au. The total energy with respect to the bottom of the well as a function of time is shown in the inset.

4. CONCLUSIONS We have described an approach of constructing compact timedependent Gaussian bases (GB) guided by the quantum or Bohmian-like trajectories (QT). Such QTGB are expected to be particularly efficient when nuclear wave functions are represented in high dimensionality, because the probability within the volume element of each QT is conserved. This means that the basis functions follow the flow of the probability density and adapt to the dynamics of a specific wave function. In principle, if exact QTs could be computed on-the-fly, a single ensemble of QTs, accurately representing a wave function at the beginning of propagation on a discrete random “grid” in coordinate space, would describe the same wave function at later times without loss of accuracy. The scaling properties of QTGB with the system size should be superior to those of the conventional time-independent bases if a wave function delocalizes in coordinate space in the course of dynamics. Although the QT dynamics itself is, generally, unstable and its numerical implementation is tedious, the QT framework is useful when it comes to dynamics close to the classical regime, typical of the nuclei in large molecular systems. This is the target regime of the QTGB approach. The difficulties in computing accurate QTs on-the-fly stem from large gradients of the quantum potential, acting on the QTs in addition to physical interactions, which generate all quantum effects in otherwise classical dynamics. In the QTGB approach, however, we do not compute the quantum potential or force, but

with the parameter values set to a = 1, b = 4, and c = 4. The particle mass is m = 1 au. The contours of the potential energy are shown in Figure 5: the two minima are located at (−√2,

Figure 5. Two-dimensional double-well potential. Positions (x, y) of the basis function centers defined by the QTs at t = 3.0 au are superimposed on the contour plot of the potential.

−√2) and (√2, √2); the barrier height is Vb = 4 Eh. The initial wave function is a direct product of two Gaussian functions, each given by eq 28, placed in the left well: Ψ(x ,y ,0) = ψ (x ,0) ψ (y ,0)

(33)

The parameters of Ψ(x,y,0) and of the propagation are listed in Table 2. The time propagation is performed with three bases of sizes 10 × 10, 12 × 12, and 16 × 16 Gaussian functions. For G

DOI: 10.1021/acs.jpca.5b10029 J. Phys. Chem. A XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry A

(3) Noguchi, D.; Tanaka, H.; Kondo, A.; Kajiro, H.; Noguchi, H.; Ohba, T.; Kanoh, H.; Kaneko, K. Quantum Sieving Effect of ThreeDimensional Cu-Based Organic Framework for H-2 and D-2. J. Am. Chem. Soc. 2008, 130, 6367−6372. (4) Sumida, K.; Stueck, D.; Mino, L.; Chai, J.-D.; Bloch, E. D.; Zavorotynska, O.; Murray, L. J.; Dinca, M.; Chavan, S.; Bordiga, S.; et al. Impact of Metal and Anion Substitutions on the Hydrogen Storage Properties of M-BTT Metal-Organic Frameworks. J. Am. Chem. Soc. 2013, 135, 1083−1091. (5) Shao, M.; Keum, J.; Chen, J.; He, Y.; Chen, W.; Browning, J. F.; Jakowski, J.; Sumpter, B. G.; Ivanov, I. N.; Ma, Y.-Z.; et al. The Isotopic Effects of Deuteration on Optoelectronic Properties of Conducting Polymers. Nat. Commun. 2014, 5, 3180. (6) Hu, S.; Lozada-Hidalgo, M.; Wang, F. C.; Mishchenko, A.; Schedin, F.; Nair, R. R.; Hill, E. W.; Boukhvalov, D. W.; Katsnelson, M. I.; Dryfe, R. A. W.; et al. Proton Transport Through One-AtomThick Crystals. Nature 2014, 516, 227. (7) Kosloff, D.; Kosloff, R. a Fouirier Method Solution for the TimeDependent Schrodinger Equation As a Tool in Molecular Dynamics. J. Comput. Phys. 1983, 52, 35−53. (8) Leforestier, C.; Bisseling, R.; Cerjan, C.; Feit, M.; Friesner, R.; Guldberg, A.; Hammerich, A.; Jolicard, G.; Karrlein, W.; Meyer, H.; et al. a Comparison of Different Popagation Schemes for the TimeDependent Schrodinger Equation. J. Comput. Phys. 1991, 94, 59−80. (9) Kosloff, R. The Fourier Method. Numerical Gird Methods and Their Applications to Schrodinger Equation; NATO Advanced Research Workshop on Grid Methods in Atomic and Molecular Quantum Calculations, Corte, France, Sep 27−Oct 03, 1992; NATO: Brussels, 1993; pp 175−194. (10) Kosloff, R. Propagation Methods Fo Quantum Molecular Dynamics. Annu. Rev. Phys. Chem. 1994, 45, 145−178. (11) Billeter, S.; VanGunsteren, W. a Comparison of Different Numerical Propagation Schemes for Solving the Time-Dependent Schrodinger Equation in the Position Representation in One Dimension for Mixed Quantum- and Molecular Dynamics Simulations. Mol. Simul. 1995, 15, 301−322. (12) Rassolov, V. A.; Garashchuk, S. Computational Complexity in Quantum Chemistry. Chem. Phys. Lett. 2008, 464, 262−264. (13) Meyer, H. D.; Manthe, U.; Cederbaum, L. S. the MultiConfigurational Time-Dependent Hartree Approach. Chem. Phys. Lett. 1990, 165, 73−78. (14) Burghardt, I.; Meyer, H.-D.; Cederbaum, L. S. Approaches to the Approximate Treatment of Complex Molecular Systems by the Multiconfiguration Time-Dependent Hartree Method. J. Chem. Phys. 1999, 111, 2927−2939. (15) Meyer, H. D.; Worth, G. A. Quantum Molecular Dynamics: Propagating Wavepackets and Density Operators Using the Multiconfiguration Time-Dependent Hartree Method. Theor. Chem. Acc. 2003, 109, 251−267. (16) Wang, H. B.; Thoss, M. Multilayer Formulation of the Multiconfiguration Time-Dependent Hartree Theory. J. Chem. Phys. 2003, 119, 1289−1299. (17) Burghardt, I.; Nest, M.; Worth, G. A. Multiconfigurational System-Bath Dynamics Using Gaussian Wave Packets: Energy Relaxation and Decoherence Induced by a Finite-Dimensional Bath. J. Chem. Phys. 2003, 119, 5364−5378. (18) Roemer, S.; Ruckenbauer, M.; Burghardt, I. Gaussian-Based Multiconfiguration Time-Dependent Hartree: A Two-Layer Approach. I. Theory. J. Chem. Phys. 2013, 138, 064106. (19) Roemer, S.; Burghardt, I. Towards a Variational Formulation of Mixed Quantum-Classical Molecular Dynamics. Mol. Phys. 2013, 111, 3618−3624. (20) Davis, M.; Heller, E. Semi-Classical Gaussian-Basis Set Method for Molecular Vibrational Wave-Functions. J. Chem. Phys. 1979, 71, 3383−3395. (21) Garashchuk, S.; Light, J. C. Quasirandom Distributed Gaussian Bases for Bound Problems. J. Chem. Phys. 2001, 114, 3929−3939.

increment the QTs in time according to the momentum reconstructed from the wave function expanded in a basis according to eq 20. This (i) generates a feedback loop between the QTs serving as the Gaussian centers and the timedependent wave function without solving numerically challenging equations, often encountered in the variational approaches, and (ii) allows for sparse sampling of the wave function with QTs. Moreover, because we use real Gaussian basis functions, the accuracy of the wave function expansion does not depend on the choice of the trajectory momenta. Thus, it is possible to use modified QT momenta, for example, defined by the convolution procedure described above, to improve stability of propagation, provided that occasional reexpansions are performed to maintain completeness of the basis in time. The completeness can be estimated (and reexpansions can be performed as needed) by examining the wave function energy, because the energy is conserved in the limit of a complete timedependent basis. In this work we have illustrated the concept of QTGB on low-dimensional model problems. In future high-dimensional tests and applications several improvements will be made. Predefined criteria and algorithms for expansion and reexpansion of wave functions in QTGB are essential for practical treatment of delocalized wave functions, as has been demonstrated by the scattering application. A well-defined procedure of setting the trajectory momenta using convolution of a wave function or some other approach is needed for stable propagation. A more conceptual development will be adaptation of the basis function shape to dynamics, such as use of correlated Gaussians whose width matrix A may depend on the coordinates and have nonzero off-diagonal elements. In the proof-of-concept applications presented above, the basis functions adapt their positions but their shape remains “frozen” throughout the propagation and is the same for all functions. Making their widths commensurate with the distance between the QTs will take advantage of the probability conservation property of the QT dynamics yielding more compact QTguided Gaussian bases and will reduce their degeneracy due to strongly overlapping basis functions. Development of the QTGB dynamics and high-dimensional applications will be further pursued in future work.

■

AUTHOR INFORMATION

Corresponding Author

*S. Garashchuk. Phone: +1 803-777-8900. Fax: +1 803-7779625. E-mail: [email protected]. Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS This material is based upon work partially supported by the National Science Foundation under Grant No. CHE-1056188.

■

REFERENCES

(1) Oh, H.; Savchenko, I.; Mavrandonakis, A.; Heine, T.; Hirscher, M. Highly Effective Hydrogen Isotope Separation in Nanoporous Metal-Organic Frameworks with Open Metal Sites: Direct Measurement and Theoretical Analysis. ACS Nano 2014, 8, 761−770. (2) FitzGerald, S. A.; Pierce, C. J.; Rowsell, J. L. C.; Bloch, E. D.; Mason, J. A. Highly Selective Quantum Sieving of D-2 from H-2 by a Metal-Organic Framework As Determined by Gas Manometry and Infrared Spectroscopy. J. Am. Chem. Soc. 2013, 135, 9458−9464. H

DOI: 10.1021/acs.jpca.5b10029 J. Phys. Chem. A XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry A

(47) Garashchuk, S.; Rassolov, V. A. Energy conserving approximations to the quantum potential: Dynamics with linearized quantum force. J. Chem. Phys. 2004, 120, 1181−1190. (48) Tannor, D. J. Introduction to Quantum Mechanics: A TimeDependent Perspective; University Science Books: Sausolito, CA, 2006. (49) Jackle, A.; Meyer, H. Product Representation of Potential Energy Surfaces. J. Chem. Phys. 1996, 104, 7974−7984. (50) Jackle, A.; Meyer, H. Product Representation of Potential Energy Surfaces. II. J. Chem. Phys. 1998, 109, 3772−3779. (51) Otto, F. Multi-Layer Potfit: An Accurate Potential Representation for Efficient High-Dimensional Quantum Dynamics. J. Chem. Phys. 2014, 140, 014106. (52) Kosloff, R. Time-Dependent Quantum-Mechanical Methods for Molecular Dynamics. J. Phys. Chem. 1988, 92, 2087−2100. (53) Feit, M. D.; Fleck, J. A., Jr.; Steiger, A. Solution of the Schrödinger Equation by a Spectral Method. J. Comput. Phys. 1982, 47, 412−433. (54) Garashchuk, S.; Dell’Angelo, D.; Rassolov, V. A. Classical Limit of Quantum Nuclear Dynamics Based on Selective Use of the Quantum Potential. J. Chem. Phys. 2014, 141, 234107. (55) Mandelshtam, V. A.; Taylor, H. S. Harmonic Inversion of Time Signals and Its Applications. J. Chem. Phys. 1997, 107, 6756−6769. (56) Mandelshtam, V. A.; Taylor, H. S. Harmonic inversion of time signals and its applications (vol 107, pg 6756, 1997). J. Chem. Phys. 1998, 109, 4128−4128.

(22) Poirier, B.; Light, J. Efficient Distributed Gaussian Basis for Rovibrational Spectroscopy Calculations. J. Chem. Phys. 2000, 113, 211−217. (23) Poirier, B. Using Wavelets to Extend Quantum Dynamics Calculations to Ten or More Degrees of Freedom. J. Theor. Comput. Chem. 2003, 2, 65−72. (24) Poirier, B.; Salam, A. Quantum Dynamics Calculations Using Symmetrized, Orthogonal Weyl-Heisenberg Wavelets with a Phase Space Truncation Scheme. II. Construction and Optimization. J. Chem. Phys. 2004, 121, 1690−1703. (25) Halverson, T.; Poirier, B. Accurate Quantum Dynamics Calculations Using Symmetrized Gaussians on a Doubly Dense Von Neumann Lattice. J. Chem. Phys. 2012, 137, 224101. (26) Halversona, T.; Poirier, B. Large Scale Exact Quantum Dynamics Calculations: Ten Thousand Quantum States of Acetonitrile. Chem. Phys. Lett. 2015, 624, 37−42. (27) Shimshovitz, A.; Tannor, D. J. Phase-Space Approach to Solving the Time-Independent Schrodinger Equation. Phys. Rev. Lett. 2012, 109, 070402. (28) Heller, E. J. Time-Dependent Approach to Semiclassical Dynamics. J. Chem. Phys. 1975, 62, 1544. (29) Heller, E. J. Frozen Gaussians: A Very Simple Semiclassical Approximation. J. Chem. Phys. 1981, 75, 2923. (30) Heller, E. J. Cellular dynamics - A new semiclassical approach to time-dependent quantum mechanics. J. Chem. Phys. 1991, 94, 2723− 2729. (31) Huber, D.; Heller, E. J. Generalized Gaussian Wave Packet Dynamics. J. Chem. Phys. 1987, 87, 5302−5311. (32) Shalashilin, D. V.; Child, M. S. Time Dependent Quantum Propagation in Phase Space. J. Chem. Phys. 2000, 113, 10028−10036. (33) Shalashilin, D. V.; Child, M. S. Real Time Quantum Propagation on a Monte Carlo Trajectory Guided Grids of Coupled Coherent States: 26D Simulation of Pyrazine Absorption Spectrum. J. Chem. Phys. 2004, 121, 3563−3568. (34) Wu, Y. H.; Batista, V. S. Matching-Pursuit for Simulations of Quantum Processes. J. Chem. Phys. 2003, 118, 6720−6724. (35) Wu, Y. H.; Batista, V. S. Matching-Pursuit Split-Operator Fourier-Transform Simulations of Excited-State Intramolecular Proton Transfer in 2-(2(’)-Hydroxyphenyl)-Oxazole. J. Chem. Phys. 2006, 124, 224305. (36) Ben-Nun, M.; Quenneville, J.; Martinez, T. J. Ab Initio Multiple Spawning: Photochemistry from First Principles Quantum Molecular Dynamics. J. Phys. Chem. A 2000, 104, 5161−5175. (37) Koch, W.; Frankcombe, T. J. Basis Expansion Leaping: A New Method to Solve the Time-Dependent Schrodinger Equation for Molecular Quantum Dynamics. Phys. Rev. Lett. 2013, 110, 263202. (38) Saller, M. A. C.; Habershon, S. Basis Set Generation for Quantum Dynamics Simulations Using Simple Trajectory-Based Methods. J. Chem. Theory Comput. 2015, 11, 8−16. (39) Herman, M. F.; Kluk, E. a Semiclassical Justification for the Use of Non-Spreading Wavepackets in Dynamics Calculations. Chem. Phys. 1984, 91, 27. (40) Kay, K. G. Integral expressions for the semiclassical timedependent propagator. J. Chem. Phys. 1994, 100, 4377−4392. (41) Hamilton, I. P.; Light, J. C. On Distributed Gaussian Bases for Simple Model Multidimensional Vibrational Problems. J. Chem. Phys. 1986, 84, 306. (42) Habershon, S. J. Chem. Phys. 2012, 136, 014109. (43) Shalashilin, D. V.; Burghardt, I. Gaussian-Based Techniques for Quantum Propagation from the Time-Dependent Variational Principle: Formulation in Terms of Trajectories of Coupled Classical and Quantum Variables. J. Chem. Phys. 2008, 129, 084104. (44) Madelung, E. Quantum Theory in Hydrodynamic Form. Eur. Phys. J. A 1927, 40, 322−326. (45) de Broglie, L. An Introduction to the Study Ot Wave Mechanics; E. P. Dutton and Company, Inc.: New York, 1930. (46) Bohm, D. A. Suggested Interpretation of the Quantum Theory in Terms of ”hidden” Variables, I and II. Phys. Rev. 1952, 85, 166−193. I

DOI: 10.1021/acs.jpca.5b10029 J. Phys. Chem. A XXXX, XXX, XXX−XXX