Moving Boundary Truncated Grid Method for Wave Packet Dynamics

Jan 17, 2018 - In addition to the boundary conditions ψij(t) = 0 for i = 1, i = N, j = 1, or j = M, ... As displayed in. Figure 1(b), if the extrapol...
0 downloads 0 Views 4MB Size
Subscriber access provided by UNIV OF OREGON

Article

Moving Boundary Truncated Grid Method for Wave Packet Dynamics Tsung-Yen Lee, and Chia-Chun Chou J. Phys. Chem. A, Just Accepted Manuscript • DOI: 10.1021/acs.jpca.7b11932 • Publication Date (Web): 17 Jan 2018 Downloaded from http://pubs.acs.org on January 18, 2018

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a free service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are accessible to all readers and citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

The Journal of Physical Chemistry A is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 31 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

Moving Boundary Truncated Grid Method for Wave Packet Dynamics Tsung-Yen Lee and Chia-Chun Chou∗ Department of Chemistry, National Tsing Hua University, Hsinchu, 30013, Taiwan (Dated: January 12, 2018)

Abstract The moving boundary truncated grid method is developed to significantly reduce the number of grid points required for wave packet propagation. The time-dependent Schr¨ odinger equation (TDSE) and the imaginary time Schr¨ odinger equation (ITSE) are integrated using an adaptive algorithm which economizes on the number of grid points. This method employs a variable number of grid points in the Eulerian frame (grid points fixed in space) and adaptively defines the boundaries of the truncated grid. The truncated grid method is first applied to the time integration of the TDSE for the photodissociation dynamics of NOCl and a three-dimensional quantum barrier scattering problem. The time-dependent truncated grid precisely captures the wave packet evolution for the photodissociation of NOCl and finely adjusts according to the process of the wave packet bifurcation into reflected and transmitted components for the barrier scattering problem. The truncated grid method is also applied to the time integration of the ITSE for the eigenstates of quantum systems. Compared to the full grid calculations, the truncated grid method requires fewer grid points to achieve high accuracy for the stationary state energies and wave functions for a two-dimensional double well potential and the Ar trimer. Therefore, the truncated grid method demonstrates a significant reduction in the number of grid points needed to perform accurate wave packet propagation governed by the TDSE or the ITSE.

ACS Paragon 1Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

I.

INTRODUCTION

Wave packet evolution has provided considerable insight into a wide range of dynamical processes in chemical physics, including reactive scattering, photodissociation, electronic nonadiabiatic dynamics, and the field of femtochemistry.1 In conventional computational techniques, the dynamics of a wave packet is obtained by solving the time-dependent Schr¨odinger equation (TDSE) using a large basis set or on a large space fixed grid. In order to avoid the traditional exponential scaling of computational effort with respect to the number of degrees of freedom, various trajectory approaches have been developed. In real and complex quantum trajectory methods, fewer grid points are needed for the time integration of the TDSE because these grid points move along quantum trajectories with velocities matching the flow velocity of the probability fluid.2–11 A synthetic quantum trajectory method based on the bipolar decomposition of the wave function has been developed for wave packet scattering problems.12–15 In contrast with the Eulerian and Lagrangian pictures, the arbitrary Lagrangian-Eulerian (ALE) picture that allows for adaptive grid point motion provides a different computational approach to the TDSE.16–18 The ALE method has been combined with the introduction of artificial viscosity in the quantum hydrodynamic equations of motion to obtain accurate time-dependent transmission probabilities for wave packet scattering.19–24 Moreover, several algorithms have been implemented for adjusting the number and distribution of points used in non-moving grids.25–27 For example, localized basis sets are added or eliminated based on fixed density parameters28,29 or time-dependent density parameters relative to the current maximum density.28 In addition to the trajectory or moving grid methods for the time integration of the TDSE, the quantum trajectory formulation has been developed to integrate the imaginary time Schr¨odinger equation (ITSE) for the stationary state energies and wave functions of quantum systems.30–32 The imaginary time quantum dynamics with the momentum-dependent quantum potential has been used to calculate the zero-point energies for chemical systems and model clusters consisting of up to 11 atoms.33,34 In the hydrodynamic formulation of quantum mechanics,2 quantum trajectories follow changes in the time-evolving probability density. This is an important feature of the quantum trajectory method for efficient description of large amplitude motion or reactive dynamics. However, the Lagrangian imaginary time quantum trajectories are governed by the inverted classical potential, and they fall ACS Paragon 2Plus Environment

Page 2 of 31

Page 3 of 31 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

off the barrier top and leave the region of high probability density. Thus, the imaginary time quantum trajectories diverge with time, and this leads to the inefficient sampling or trajectory representation for quantum mechanical bound state problems. In the current study, the moving boundary truncated grid method is developed to significantly reduce the number of grid points required for accurate wave packet propagation. This method employs a variable number of grid points in the Eulerian frame and adaptively defines the boundaries of the truncated grid. At each time step, the computational domain is truncated according to the values of the probability density and the spatial derivatives of the wave function at the grid points. The information of the wave packet within the truncated grid is employed to extrapolate approximate values of the wave function at additional grid points outside the grid boundaries. Based on the wave function in the truncated grid and the extrapolation results, the wave packet is then propagated using conventional time integrators. The process is repeated until the final solution is reached. The truncated grid method is first applied to the time integration of the TDSE for the photodissociation dynamics of NOCl and a three-dimensional quantum barrier scattering problem. The time-dependent truncated grid precisely captures the wave packet evolution for the photodissociation of NOCl, and the resulting energy-resolved absorption spectrum is in excellent agreement with the full grid calculation. In the barrier scattering problem, the truncated grid finely adjusts according to the process of the wave packet bifurcation into reflected and transmitted components, and excellent time-dependent transmission probabilities are obtained. On the other hand, the truncated grid method is also applied to the time integration of the ITSE for the eigenstates of a two-dimensional double well potential and the Ar trimer. Compared to the full grid calculations, the truncated grid method requires fewer grid points to achieve high accuracy for the calculations of energies and wave functions for the ground and first-excited states of the double well potential. Moreover, the ground state energy and wave function of the Ar trimer are obtained using an extremely small number of grid points in the three-dimensional computational domain. These computational results demonstrate that appropriate Eulerian grid points containing significant density are activated and deactivated during the evolution of the wave packet without any advance knowledge of the dynamics. Therefore, the truncated grid method proposed in this study leads to a significant reduction in the number of grid points required for the time integration of the TDSE and the ITSE. ACS Paragon 3Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 4 of 31

The outline of the remainder of this study is as follows. In Sec. II, the truncated grid method is developed. The algorithms for the truncation of the computational grid, the extrapolation of the wave function, and the propagation of the wave packet are describe. In Sec. III, applications of the truncated grid method are made to the time integration of the TDSE and the ITSE. Computational results and analysis are presented for the photodissociation dynamics of NOCl and a three-dimensional quantum barrier scattering problem. In addition, results based on the integration of the ITSE are presented for the eigenstates of a two-dimensional double well potential and the Ar trimer. Finally, we make some comments and conclude with a discussion about future research directions in Sec. IV.

II.

COMPUTATIONAL METHOD

Without loss of generality, we consider the evolution of a wave packet governed by the two-dimensional TDSE ∂ψ ~2 i~ =− ∂t 2m

(

∂ 2ψ ∂ 2ψ + 2 ∂x2 ∂y

) + V (x, y)ψ,

(1)

where ψ(x, y, t) is the wave function and V (x, y) is the potential energy surface. The initial condition for this equation is provided by the initial wave function ψ(x, y, 0). In general, the wave function asymptotically approaches zero as r → ∞, where r is the radial coordinate. In conventional numerical techniques for propagation of the wave packet, the wave function is set to zero at the boundaries of a computational grid as long as these boundaries are far enough from the interaction range. Suppose that the rectangular computational domain extends from x = a to x = b and from y = c to y = d with the grid spacings ∆x and ∆y, and it is discretized by N grid points along x direction and M grid points along y direction with x1 = a, xN = b, y1 = c, and yM = d. The wave function ψ(x, y, t) at the grid point (x, y) = (xi , yj ) is denoted by ψij (t), which approximates the solution to the TDSE. However, a large grid has to be employed to avoid spurious or unphysical reflections of the wave packet at the boundaries. In addition, grid points introduced at early times where the wave packet initiates may not be needed at later times because it has moved to a different location. These results lead to wasted effort in calculating solutions for grid points that are of little significance to the overall dynamics of the wave packet. ACS Paragon 4Plus Environment

Page 5 of 31 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

(b)

(a)

TB

PG

TBE

TG

Target

Target

FIG. 1. (Color online) (a) The truncated grid (TG), the boundary of the truncated grid (TB ), and the propagation grid (PG); (b) Target points are those points with zero probability density near the TBE points. The extrapolation process is applied across the boundary of the truncated grid to obtain the approximate wave function at the target points.

A.

Truncation of the computational grid

In order to concentrate on the significant portions of the time-evolving wave packet, we choose a cutoff value ρmin for the probability density. In addition to the boundary conditions ψij (t) = 0 for i = 1, i = N , j = 1, or j = M , if the probability density at a grid point is smaller than the cutoff value (ρij (t) = |ψij (t)|2 < ρmin ), then the wave function at this point is set to zero. However, the probability density also becomes very small near nodal regions. It is noted that the absolute value of the spatial derivative of the wave function generally takes on large values near these regions. Thus, another parameter dmin is introduced to avoid the removal of grid points near nodal regions. Hence, the wave function is set to zero at the grid points satisfying all the following conditions:

ρij (t) < ρmin , ∂ψ(xi , yj , t) < dmin , and ∂ψ(xi , yj , t) < dmin . ∂x ∂y

(2)

In this way, the truncated grid (TG) can capture the main features of the time-evolving wave packet, as shown in Fig. 1(a). ACS Paragon 5Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 6 of 31

FIG. 2. (Color online) The initial wave packet and contours of the S1 electronic potential energy surface as a function of the coordinates R and r at the equilibrium angle θeq = 127.4◦ for the photodissociation dynamics of NOCl. The boundary points of the truncated grid are shown as black dots.

B.

Extrapolation for the wave packet

As each time step, the information of the wave packet within the truncated grid is employed to extrapolate approximate values of the wave function at additional grid points outside the grid boundaries. These extrapolation results will be used for time propagation of the wave packet. As shown in Fig. 1(a), the set of the boundary points for the truncated grid is denoted by TB. However, not all the boundary points in the set TB are considered for the extrapolation process. As the truncation criteria for the computational grid presented in Sec. II A, two additional parameters ρext and dext are introduced. We focus on the set TBE containing those boundary points satisfying one of the following conditions: ρij (t) ≥ ρext , ∂ψ(xi , yj , t) ∂ψ(xi , yj , t) ≥ dext , or ≥ dext , ∂x ∂y

(3)

where (xi , yj ) belongs to the set TB. The set TBE does not include the grid points with i = 2, 3, i = N − 1, N − 2, j = 2, 3, or j = M − 1, M − 2. Thus, the extrapolation process is not performed at these outermost grid points. As shown in Fig. 1(b), the grid points with zero probability density near the TBE points are referred to as target points. The wave packet is extrapolated across the boundary of ACS Paragon 6Plus Environment

Page 7 of 31 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

the truncated grid to determine the wave function at the target points. The extrapolation method begins with invoking the exponential form of the wave function, ψ = exp(A), where A is a complex-valued function. The values of the wave function at two grid points are used to perform the linear extrapolation for A at the target point. Even when the wave function ψ is rapidly oscillating, A is frequently slowly varying with respect to the x and y coordinates. As displayed in Fig. 1(b), if the extrapolation can be performed along several directions, the approximate wave function value at this point is given by the average of these extrapolation results. In addition, in order to increase the stability of the extrapolation method, if the absolute value of the wave function at the target point is greater than the smallest absolute value of the wave function in the surrounding grid points with nonzero probability densities (say |ψmin |), then the extrapolated wave function at the target point is evaluated by ψtarget = ψmin e−ηh

(4)

where η is a reducing factor and h is the corresponding grid spacing. This modification can reduce errors arising from the extrapolation method and avoid sudden and dramatic changes in the wave function near the boundary of the truncated grid. In this study, η = 1 (in atomic units) will be used in the computational examples presented in Sec. III.

C.

Time propagation of the wave packet

The time integrator, such as the Euler method and the fourth-order Runge-Kutta method, is applied to the target points to obtain the temporary wave function ϕij at time t + ∆t. It is noted that the actual wave packet may expand or contract beyond one grid spacing within a time step ∆t. If the wave function ϕij at a target point satisfies one of the following conditions: |ϕij |2 ≥ ρmin , ∂ϕ(xi , yj , t) ∂ϕ(xi , yj , t) ≥ dmin , or ≥ dmin , ∂x ∂y

(5)

then we add the target point into the set TBE. The extrapolation method described in Sec. II B is applied to these new TBE points, based on their original extrapolation results. The process is repeated until no new TBE points are generated. Then, we include the final target points and their neighboring grid points into the truncated grid to form a propagation ACS Paragon 7Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

FIG. 3. (Color online) Snapshots of the amplitude of the wave packet at four times for the photodissociation of NOCl on the excited potential energy surface. The boundary points of the truncated grid are shown as black dots.

grid (PG), as shown in Fig. 1(b). Finally, the time integrator is applied to the propagation grid to obtain the wave function at the next time step, ψij (t + ∆t) = ϕij . We briefly summarize the computational procedure. At each instant of time t, the computational grid is truncated according to the criteria in Eq. (2). From the wave function in the truncated grid, the wave packet is extrapolated to obtain the wave function at the target points. The wave function at the target points is evolved one time step. If the temporary wave function at a target point satisfies one of the conditions in Eq. (5), we add the target point into the set TBE and redo the extrapolation of the wave packet. The extrapolation process is repeated until the wave packet no longer expands or contracts (there is no new TBE point). Then, the time integrator is applied to the propagation grid to obtain the final wave packet at time t + ∆t. Thus, the wave function have been updated, and then we use ACS Paragon 8Plus Environment

Page 8 of 31

Page 9 of 31

(b)

Truncated Grid Pts (%)

(a)

40 t (fs)

20 0

20 15 10 5

4 r

3 2

4

8

6

0 0

R

10

20

30 t (fs)

40

50

FIG. 4. (Color online) (a) Time evolution of the boundary of the truncated grid for the photodissociation dynamics of NOCl; (b) The percentage of grid points in the computational domain used in the truncated grid as a function of time.

(a)

(b) −4.0 −4.2 log (Error)

σ

−4.4

10

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

−4.6 −4.8 −5.0

0

0.5

1

1.5 E (e.V.)

2

2.5

0

10

20

30 t (fs)

40

50

FIG. 5. (Color online) (a) The energy-resolved absorption spectrum for NOCl derived from the truncated grid method (blue circles) and the full grid calculation (black curve); (b) Log of the relative error between the wave functions obtained from the truncated grid method and the full grid calculation as a function of time.

this information as the initial condition to advance one more time step. ACS Paragon 9Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

III.

Page 10 of 31

COMPUTATIONAL RESULTS

In this section, the moving boundary truncated grid method will be applied to the time evolution of wave packets governed by the TDSE or the ITSE.

A. 1.

Time-dependent Schr¨ odinger equation Photodissociation dynamics of NOCl

The first example deals with the photodissociation dynamics of NOCl on an electronically excited potential energy surface.35–41 The dynamics of the system is described by the Jacobi coordinates R, r, and θ, where R is the distance (dissociative coordinate) from Cl to the center of mass of NO, r is the bond length (vibrational coordinate) of NO, and θ indicates the angle between r and R. The Hamiltonian acting on the wave function for the NOCl molecule takes the form

[ ~2 1 ∂2 1 ∂2 Hψ = − (Rψ) + (rψ) 2 md R ∂R2 mv r ∂r2 ( )] 1 ∂ ∂ψ + sin θ + V (R, r, θ)ψ. Iθ sin θ ∂θ ∂θ

(6)

The reduced masses md and mv and the moment of inertia Iθ for the NOCl molecule are given by 1 1 1 = + , md mN + mO mCl 1 1 1 = + , mv mN mO 1 1 1 = + , 2 Iθ md R mv r 2

(7) (8) (9)

where mN , mO , and mCl are the masses of atoms N, O, and Cl, respectively. The coordinates and the Hamiltonian are the same for both the S0 and S1 potential energy surfaces. The NOCl system is reduced to two coordinates by freezing the angular variable at the equilibrium value θeq = 127.4◦ , and the S0 and S1 surfaces for NOCl have the same equilibrium value.37 Thus, we obtain the two-dimensional TDSE ] [ ∂ψ ~2 1 ∂2 1 ∂2 i~ =− (Rψ) + (rψ) ∂t 2 md R ∂R2 mv r ∂r2 + V (R, r, θeq )ψ. ACS Paragon10 Plus Environment

(10)

Page 11 of 31 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

This equation can be further simplified by introducing the dimensionless quantities38,39 2 2 mv req mv req ϕ = Rrψ, U = V, t0 = , ~2 ~ ( )1/2 R md r t , y= , τ= , x= mv req req t0

where md = 29661 and mv = 13710. Atomic units (~ = 1) are used throughout this paper, except where indicated otherwise. The constant req = 2.155 is the value of the vibrational coordinate at the molecular equilibrium on the S0 surface. Hence, the TDSE in Eq. (10) becomes the dimensionless version given by [ ] ∂ϕ 1 ∂ 2ϕ ∂2ϕ i =− + + U (x, y)ϕ, ∂τ 2 ∂x2 ∂y 2

(11)

where x and y range from zero to infinity. In the time-dependent study of photodissociation, the initial ground-state vibrational wave function on the S0 surface is vertically excited to the S1 surface within the FranckCondon approximation. The excitation from the S0 to S1 surface is assumed to occur instantaneously via a laser field. Thus, the initial wave packet on the S1 surface has the same shape as that on the S0 surface. Then, the dynamics of photodissociation is described by the time evolution of this initial bound state on the excited electronic surface. Because the S0 surface has a deep well, the initial wave packet determined by the ground state wave function of the S0 surface under the quadratic approximation is given by37–39 ϕ0 (x, y) =

( ω ω )1/4 { [ 1 2 exp −(1/2) a(x − xeq )2 2 π ]} +2b(x − xeq )(y − yeq ) + c(y − yeq )2 ,

(12)

where xeq = 2.944, yeq = 1, ω1 = 527.1, ω2 = 498.6, a = 506.1, b = 12.53, and c = 519.7. The S1 surface for NOCl is taken as39 U (x, y) =a2 (y − y1eq )2 + c00 (1 − Q(x)) [1 + c01 Q(x) + c10 (y − y1eq ) + c20 (y − y1eq )2 + c11 Q(x)(y − y1eq )

] + c21 Q(x)(y − y1eq )2 , ACS Paragon11 Plus Environment

(13)

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 12 of 31

where Q(x) = 1 − exp[−α(x − x1eq )]. The associated parameters are given by x1eq = 2.945, y1eq = 0.991, α = 2.198, a2 = 0.202 × 106 , c00 = 2450, c10 = −3.129, c01 = 0.217, c20 = −19.69, c11 = 3.260, c21 = −3.508. This simplified S1 potential and the parameters are identical to those used in Ref. 39. The initial wave packet and contours of the S1 surface are displayed in Fig. 2. The computational domain extends from x = 1.5 to x = 7.0 and y = 0.6 and y = 2.1 with the grid spacings ∆x = ∆y = 0.01, and the spatial derivatives in the two-dimensional TDSE in Eq. (11) were calculated using the second-order central finite difference. Time propagation was performed using the fourth-order Runge-Kutta method from t = 0 (0 fs) to t = 2400 (58 fs) with the time step ∆t = 0.05. The parameters for the grid truncation and the extrapolation of the wave function are given by ρmin = 0.001, dmin = 0.01, ρext = 0.01, and dext = 0.1, respectively. The boundary points of the truncated grid for the initial wave packet are shown in Fig. 2, and the truncated grid captures the significant portion of the wave function. Figure 3 displays the time evolution of the wave packet with the boundary of the truncated grid. As time evolves, the wave packet gradually spreads, and a large portion of the packet rapidly travels along the potential valley with a downward slope into the NO+Cl dissociation channel. In addition, Fig. 3 demonstrates that the time-dependent adaptive truncated grid follows the main features of the evolving wave packet. Fig. 4(a) presents the time evolution of the boundary of the truncated grid. As time advances, the truncated grid continues to expand and constantly changes its shape to capture the intricacies of the wave packet. As shown in Fig. 4(b), the truncated grid method significantly reduces the number of grid points for the photodissociation dynamics of NOCl compared to the full grid calculation. Even at the final time t = 58 fs, only 21% of the grid points in the computational domain were used. The energy-dependent absorption spectrum can be obtained by the Fourier transform of the time-dependent autocorrelation function C(t) given by [∫ ] 2 iωt −t2 /τw σ(ω) = ω Re C(t) e e dt ,

(14)

where C(t) = ⟨ϕ(0)|ϕ(t)⟩ is the overlap integral between the evolving wave packet at time ACS Paragon12 Plus Environment

Page 13 of 31 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

t and the initial wave packet at t = 0. The quantity τw = 67 fs characterizes the finite resolution of the spectrum because of the limited propagation time.42 The choice of the window function corresponds to a convolution of the spectrum with a Gaussian. Figure 5(a) presents the spectrum for NOCl obtained from the truncated grid method and the full grid calculation, and these two results are in excellent agreement. Moreover, in order to assess the accuracy of the truncated grid calculation, we consider the relative error defined by Error =

∥ ϕtruncated − ϕf ull ∥ , ∥ ϕf ull ∥

(15)

where ∥ · ∥ is the L2 norm. Figure 5(b) presents the log of the time-dependent error. At t = 0, the discrepancy between the truncated and full grid calculations arises from the truncation of the initial wave packet. As time proceeds, the error gradually increases; however, its value is less than 0.00015 for all times. Therefore, the truncated grid method can significantly reduce the total number of grid points relative to the full grid calculation, while the solution accuracy is retained for the photodissociation dynamics of NOCl.

2.

Three-dimensional barrier scattering

In the second application of the truncated grid method, the time-dependent study of a three-dimensional quantum barrier scattering problem will be presented. The anharmonic potential energy surface is modeled by an Eckart barrier along the translational coordinate with two Morse vibrational modes in the directions perpendicular to the reaction path,24 V (x, y, z) = V0 sech2 (αx) + De [1 − exp(−βy)]2 + De [1 − exp(−βz)]2 − 2De ,

(16)

where the values of the parameters are given by V0 = 10, α = 0.4, De = 40, and β = 0.25. The Eckart barrier maximum is located at x = 0, and these Morse oscillators are centered at y = 0 and z = 0, respectively. The initial wave packet is given by ( ψ0 (x, y, z) =

2βx π

)1/4 e−βx (x−xc )

2 +ik

0x

ψM (y)ψM (z),

ACS Paragon13 Plus Environment

(17)

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

FIG. 6. (Color online) Snapshots of the time evolution of the wave packet scattering from a threedimensional barrier at six times with the boundary of the truncated grid (blue dots). The isosurface of the probability density of the wave packet is plotted in each figure for P = 0.01 (red surfaces). The two isosurfaces of the potential energy surface in each figure have the values V = −60 (yellow surfaces) and V = −63 (green surfaces).

ACS Paragon14 Plus Environment

Page 14 of 31

Page 15 of 31

(a)

(b) 1 25

E=1.50V

0

0.8

Truncated Grid Pts (%)

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

E=1.25V0 E=1.00V0

0.6

P E=0.75V

0.4

0

E=0.50V0

0.2

E=0.25V0

0 0

0.5

1

1.5

2

2.5

3

3.5

20 15 10 5 0 0

0.5

t

1

1.5

2

2.5

3

3.5

t

FIG. 7. (Color online) (a) Time-dependent transmission probabilities obtained from the full grid calculations (blue curves) and the truncated grid method (black dots) for the three-dimensional barrier scattering; (b) The percentage of grid points in the computational domain used in the truncated grid as a function of time for E = 0.25V0 (purple curve), E = 0.50V0 (blue curve), E = 0.75V0 (green curve), E = 1.00V0 (yellow curve), E = 1.25V0 (orange curve), and E = 1.50V0 (red curve).

where βx = 1.5, xc = −5, and ψM is the ground state of the Morse potential. The Morse ground state wave function is expressed as ψM (r) = exp(C(r)) with24 ) ( 1 [ln(2λ) − βr] − λ exp(−βr), (18) C(r) = ln N + λ − 2 √ √ where λ = 2mDe /β~ and N = β(2λ − 1)/Γ(2λ). The initial wave packet in Eq. (17) is centered at (x, y, z) = (xc , 0, 0), and it has a momentum ~kx toward the product region, √ where kx = 2mE/~, m = 1, and E is the initial translational energy. The computational domain extends from x = −15 to x = 25, from y = −3 to y = 5, and from z = −3 to z = 5 with the grid spacings ∆x = ∆y = ∆z = 0.1. Analogously, the fourth-order Runge-Kutta method was utilized for propagating the wave packet from t = 0 to t = 3.5 with the time step ∆t = 0.001. The parameters for the truncated grid method used in the three-dimensional calculation are as follows: ρmin = 1 × 10−7 , dmin = 1 × 10−3 , ρext = 1 × 10−6 , and dext = 5 × 10−3 . Figure 6 shows the time evolution of the probability density of the wave packet for the initial energy E = V0 with the boundary of the truncated grid. In this figure, the isosurface is shown for the probability density at P = 0.01, and two isosurfaces for the potential energy surface are also presented for V = −60 and V = −63. ACS Paragon15 Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 16 of 31

2

1

y 0

−1

−2 −2

−1

0 x

1

2

FIG. 8. (Color online) Contours of the two-dimensional double well potential in Eq. (24).

At t = 0, the initial wave packet is heading toward the barrier. For t = 1, the wave packet encounters the barrier region. As time proceeds, the wave packet starts to interact with the potential surface and some amplitudes have transmitted to the product side of the barrier, as shown in this figure for t = 1.5 and t = 2. For t = 3, the total wave packet starts to split into the reflected and transmitted parts. Finally, the transmitted and reflected wave packets are well bifurcated at t = 3.5. As shown in this figure, the truncated grid correctly identifies the necessary grid points to capture the main features of the evolving wave packet. Figure 7(a) presents the time-dependent transmission probability computed by integrating the density of the wave packet over the product region ∫









P (t) =



ρ(x, y, z, t)dxdydz. −∞

−∞

(19)

0

As shown in this figure, these time-dependent probabilities for six different initial translational energies obtained from the truncated grid method are in excellent agreement with the full grid calculations. The relative errors in the final probabilities at t = 3.5 between the truncated and full grid calculations are all less than 0.045%. In addition, Fig. 7(b) indicates that the truncated grid method significantly reduces the number of grid points for the three-dimensional barrier scattering problem compared to the full grid calculations. At the final time t = 3.5, the wave packet spreads in space and bifurcates into reflected and transmitted components; however, only 20%∼25% of the grid points in the computational domain were required for accurate wave packet propagations. ACS Paragon16 Plus Environment

Page 17 of 31 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

FIG. 9. (Color online) Two-dimensional double well potential: (a) The initial wave functions for the ground and first-excited states for the time integration of the ITSE with the boundary points of the truncated grid (black dots); (b) The final wave functions for the ground and first-excited states with the boundary points of the truncated grid (black dots).

B.

Imaginary-time Schr¨ odinger equation

The imaginary time propagation method has been developed to obtain the ground-state energy and wave function.43–49 Performing a transformation from real time to imaginary time by introducing the new variable τ = it, we obtain the imaginary time version of the TDSE ~

∂ ~2 2 ψ = −Hψ = ∇ ψ − V ψ. ∂τ 2m

(20)

The formal solution of this equation is written as ψ(τ ) = e−Hτ /~ ψ(0), ACS Paragon17 Plus Environment

(21)

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

(a)

(b) 2

2

1.5

1.5

τ 1

τ 1

0.5

0.5

0

0

4

4 2

2 y

Page 18 of 31

0 −2 −4

y −4

−2

0

x

4

2

0 −2 −4 −4

−2

0

2

4

x

FIG. 10. (Color online) Two-dimensional double well potential: Time evolution of the boundary of the truncated grid for the two-dimensional double well potential: (a) the ground state; (b) the first-excited state.

where ψ(0) is the initial wave packet. We then expand the initial state into a series expansion in terms of the eigenfunctions of H

ψ(0) =



cn ϕn ,

(22)

n

where ϕn are the eigensolutions of the TDSE, Hϕn = En ϕn , and the expansion coefficients are given by cn = ⟨ϕn |ψ(0)⟩. Substituting this expansion into the formal solution of the ITSE in Eq. (21) yields

ψ(τ ) = e−Hτ /~

∑ n

cn ϕn =



cn e−En τ /~ ϕn .

(23)

n

It follows from this equation that each eigenfunction decays exponentially to zero at a rate proportional to its eigenvalue for large τ . Thus, the ground state relaxes most slowly and persists longest. As time proceeds, the wave function ψ(τ ) converges to the ground state wave function ϕ0 regardless of the choice of the initial wave function, as long as there is a numerically significant overlap between the initial state ψ(0) and the ground state ϕ0 . Hence, propagating an arbitrary initial wave function in imaginary time projects onto the ground state. ACS Paragon18 Plus Environment

Page 19 of 31

40 Truncated Grid Pts (%)

35

30

25 0

1

2

3

τ

4

5

FIG. 11. (Color online) Two-dimensional double well potential: The percentage of grid points in the computational domain used in the truncated grid as a function of time for the ground (red curve) and first-excited (blue curve) states.

(b)

(a) 1

−3.4

0.5 E

log10(Error)

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

0

−0.5

−3.6

−1 −1.5 0

−3.5

1

2

τ

3

4

5

0

1

2

τ

3

4

5

FIG. 12. (Color online) Two-dimensional double well potential: (a) Time-dependent energy expectation values for the ground and first-excited states obtained from the truncated grid method (red and blue circles) and the full grid calculations (black curves); (b) Log of the relative errors between the wave functions obtained from the truncated grid method and the full grid calculation for the ground (red curve) and first-excited (blue curve) states as a function of time.

1.

Two-dimensional double well potential

In order to illustrate the application of the truncated grid method to the time integration of the ITSE in Eq. (20), we consider a two-dimensional double well potential given by49 V (x, y) = 3x4 − 8x2 + y 2 − xy. ACS Paragon19 Plus Environment

(24)

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 20 of 31

The potential energy is a sum of a double well potential along the reaction coordinate x and a harmonic oscillator along the bath coordinate y with a coupling term linear in both coordinates. This potential has been used to describe the intramolecular hydrogen transfer process,50,51 and its contours are displayed in Fig. 8. The computational domain extends from x = −5 to x = 5 and from y = −5 to y = 5 with the grid spacings ∆x = ∆y = 0.1, and the ITSE in Eq. (20) was integrated using the first-order Euler method from τ = 0 to τ = 5 with the time step ∆τ = 0.0001. The mass used in this example is m = 0.5. The initial states employed to determine the ground and first-excited states of the system are separately given by ( )1/2 2 2 2 e−(x +y ) , ψ0 (x, y) = π ( )1/2 8 2 2 ψ0 (x, y) = xe−(x +y ) . π

(25) (26)

The parameters for the truncated grid method are as follows: ρmin = 1 × 10−7 , dmin = 1 × 10−3 , ρext = 1 × 10−6 , and dext = 5 × 10−3 . The initial wave functions and the corresponding boundary points of the truncated grid are shown in Fig. 9(a). For the ground state, the initial wave packet gradually evolves from a Gaussian function to the ground state wave function with two peaks, as displayed in Fig. 9(b). On the other hand, this figure also indicates that the corresponding initial wave function converges to the first-excited state of the system. Figure 10 presents the time evolution of the boundary points of the truncated grid for the ground and first-excited states of the double well potential. Because the initial wave functions in Eqs. (25) and (26) rapidly converge to the eigenstates of the system, the boundaries of the truncated grid do not show significant changes after a short propagation time. As shown in Fig. 11, because the initial wave packets continue to expand before approximately τ = 1, the numbers of grid points required for the wave packet propagation gradually increase. However, after the convergence of these wave functions to the eigenstates, the total numbers of grid points in the truncated grid remain almost constant. Figure 12(a) shows the time-dependent energy expectation values for the ground and first-excited states of the double well potential. As time progresses, these energies obtained from the truncated grid method quickly converge to stable plateaus with E0 = −1.3413597 and E1 = −0.7661102, respectively. These results are in excellent agreement with the full ACS Paragon20 Plus Environment

Page 21 of 31 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

τ=0

τ = 2000

5

5

4.5

4.5

R3 4

R3 4

3.5

3.5

3

3 5

2.5

4 3 R 2.5

3

3.5 R

1

4

4.5

5

2.5 2

2.5

3

τ = 10000

3.5 R 4 1

4.5

5

5 4 3 R

2

τ = 35000

5

5

4.5

4.5

R3 4

R3 4

3.5

3.5

3

3 5

2.5 2.5

3

3.5 R 4 1

3 R 4.5

5

5

2.5

4

4 3

2

2.5

3

3.5 R

1

4

4.5

5

R2

FIG. 13. (Color online) Snapshots of the time evolution of the wave packet for the Ar trimer at four times with the boundary of the truncated grid (blue dots). The isosurface of the probability density of the wave packet is plotted in each figure for P = 0.01 (red surfaces). The coordinates R1 , R2 , and R3 are given in angstrom.

grid calculations E0 = −1.3413648 and E1 = −0.7661157. Moreover, Fig. 12(b) shows the log of the time-dependent relative errors defined by Eq. (15). At the final time τ = 5, the relative errors between the truncated and full grid calculations for the ground and firstexcited states are less than 0.0004. Hence, the application of the truncated grid method to the ITSE can yield eigenstates and energy eigenvalues with high accuracy.

2.

Ar trimer

The final example will be devoted to the calculation of the ground state for the Ar trimer. The Hamiltonian for the Ar trimer with zero total angular momentum can be expressed in ACS Paragon21 Plus Environment

The Journal of Physical Chemistry

(a)

(b) 16

-220

Truncated Grid Pts (%)

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 22 of 31

-230

E -240

14 12 10 8 6

-250 0

10000

τ

20000

30000

0

10000

20000

30000

τ

FIG. 14. (Color online) Ar trimer: (a) Time-dependent energy expectation value (given in cm−1 ) for the ground state obtained from the truncated grid method (blue circles) and the full grid calculation (black curve); (b) The percentage of grid points in the computational domain used in the truncated grid as a function of time for the ground state.

terms of the atom-atom pair coordinates R1 , R2 , and R3 as52,53 [ 3 { ∑ ~2 1 ∂ 2 ∂ Ri H= − 2 m R ∂R ∂Ri i i i=1

] } Rj2 + Rk2 − Ri2 ∂ 2 + + V (Ri ) , 2Rj Rk ∂Rj ∂Rk

(27)

where i ̸= j ̸= k and m = 72915. In these coordinates, the volume element is given by dV = R1 R2 R3 dR1 dR2 dR3 , and the normalization condition becomes ∫ ∫ ∫ |Ψ(R1 , R2 , R3 )|2 R1 R2 R3 dR1 dR2 dR3 = 1,

(28)

(29)

where Ψ(R1 , R2 , R3 ) is one of the eigenfunctions of the Hamiltonian in Eq. (27). Thus, the probability density takes the form P (R1 , R2 , R3 ) = |Ψ(R1 , R2 , R3 )|2 R1 R2 R3 .

(30)

For the Ar trimer, the potentials in Eq. (27) are given by simple Morse functions52,53 [ ] V (Ri ) = D e−2α(Ri −Re ) − 2e−α(Ri −Re ) , ACS Paragon22 Plus Environment

(31)

Page 23 of 31 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

−1 where D = 99 cm−1 , α = 1.72 ˚ A , and Re = 3.75 ˚ A. In this example, the initial wave

function is taken to be a Gaussian function Ψ(R1 , R2 , R3 ) = N e−[(R1 −R0 )

2 +(R −R )2 +(R −R )2 2 0 3 0

],

(32)

where R0 = 3.78 ˚ A and N is the normalization constant. The computational domain extends from Ri = 4 to Ri = 11 with the grid spacing ∆Ri = 0.1 for i = 1, 2, 3. The initial wave packet in Eq. (32) was then propagated with the first-order Euler scheme from τ = 0 to τ = 35000 using the time step ∆τ = 0.5. The relevant parameters for the truncated grid method are given by ρmin = 0.0001, dmin = 0.001, ρext = 0.001, and dext = 0.01. Figure 13 displays the time evolution of the wave packet at four times. As shown in this figure, the wave packet quickly contracts before τ = 10000, and then gradually approaches the ground state of the Ar trimer. As time evolves, the truncated grid adjusts accordingly such that the main features of the wave packet can be captured. In addition, we calculated the average triangle sides ⟨Ri ⟩ for the ground state of the Ar trimer. The expectation values ⟨R1 ⟩ = ⟨R2 ⟩ = ⟨R3 ⟩ = 3.83028 ˚ A obtained from the truncated grid method are in excellent agreement with the full grid results ⟨R1 ⟩ = ⟨R2 ⟩ = ⟨R3 ⟩ = 3.83051 ˚ A. Figure 14(a) shows the time-dependent energy expectation value for the ground state of the Ar trimer. As time progresses, the energy obtained from the truncated grid method rapidly converges to a stable plateau with the ground state energy E0 = −252.502 cm−1 , which is in excellent agreement with the full grid result E0 = −252.511 cm−1 . Moreover, the number of grid points in the computational domain used in the truncated grid as a function of time is presented in Fig. 14(b). Because the wave packet rapidly contracts before τ = 10000, the total number of grid points in the truncated grid considerably decreases during this time interval. When the wave packet gradually approaches the ground state of the Ar trimer, the number of truncated grid points tends to a constant. After the convergence of the wave packet to the ground state at τ = 35000, only approximately 5% of the grid points in the computational domain were used in the three-dimensional system.

IV.

DISCUSSION AND CONCLUSIONS

The moving boundary truncated grid method was developed to significantly reduce the number of grid points required for wave packet propagation. In order to concentrate on the ACS Paragon23 Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

significant portions of the wave packet, the computational domain is truncated according to the values of the probability density and the spatial derivatives of the wave function at the grid points. Through the exponential form of the wave function, the information of the wave packet within the truncated grid is employed to extrapolate approximate values of the wave function at additional grid points outside the grid boundaries. Based on the wave function in the truncated grid and the extrapolation results, the wave packet can be propagated using conventional time integrators. Then, the method was applied to the time integration of the TDSE for the photodissociation dynamics of NOCl and a three-dimensional quantum barrier scattering problem and to the time integration of the ITSE for the eigenstates of a two-dimensional double well potential and the Ar trimer. Computational results indicate that the truncated grid method provides an adaptive algorithm which economizes on the number of grid points in fixed grids required for the time integration of the TDSE or the ITSE. In the truncated grid method, we need to choose several parameter values for the truncation of the computational grid and the extrapolation of the wave packet. The choice for the value ρmin depends mainly on the desired accuracy of solutions and the total number of grid points in the truncated grid. A small value for ρmin is required for obtaining solutions with high accuracy. The grid truncation also involves the parameter dmin for the spatial derivatives of the wave function. First, we can choose a large value for dmin ; thus, the boundary of the truncated grid in this case is primarily determined by ρmin . The corresponding boundary points essentially have ρ > ρmin and d < dmin . Then, we gradually reduce the value of dmin such that the truncated grid is still approximately identical to that determined by ρmin ; however, the boundary points mostly become ρ < ρmin and d > dmin . Through this process, the boundary points of the truncated grid originally dominated by ρmin gradually develop into those mainly dominated by dmin . In this way, we obtain the appropriate values for the parameters ρmin and dmin . Regarding the choice of the parameters for the extrapolation of the wave packet, experience so far indicates that setting ρext = 2ρmin ∼10ρmin and dext = 2dmin ∼10dmin can provide stable extrapolation results. The parameters ρmin and dmin are used to determine whether a grid point is included in the truncated grid, while the other two parameters ρext and dext are employed to determine whether the extrapolation of the wave packet is performed. As a result, the values of ρext and dext are slightly larger than those for ρmin and dmin . ACS Paragon24 Plus Environment

Page 24 of 31

Page 25 of 31 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

One of the most popular methods for quantum dynamics is the split operator method based on fast Fourier transforms (FFTs).54–57 There exist some issues concerning the application of the moving boundary truncated grid method to the split operator method. In FFT algorithms, the vector length in position space must match that in momentum space. Although we may reduce the number of points in momentum space to match the number of points in the truncated grid, this process decreases the stability of the numerical scheme and the accuracy of solutions. On the other hand, we may expand the truncated grid with zeroes to maintain the appropriate number of momentum points required for a stable propagation, but there is no efficiency gained by placing zeroes around the truncated grid. For the computational examples presented in this study for the integration of the TDSE and the ITSE, although the number of grid points is significantly reduced, the moving boundary truncated grid method needs slightly more computational time than the full grid calculation. At each time step, we employ the information of the wave packet within the truncated grid to perform the extrapolation of wave function and to determine the boundary of the truncated grid. This process leads to a 10%–150% increase in computational time, depending on the system and its dimensionality. It is expected that the computational time used in the truncated grid method greatly decreases with increasing dimensionality. Thus, the truncated grid method may provide an alternative way to avoid the traditional exponential scaling of computational effort with respect to the number of degrees of freedom. As the computational results in the current study demonstrate, the moving boundary truncated grid method not only significantly reduces the number of grid points for wave packet propagation but also correctly identifies the necessary grid points to capture the main features of the evolving wave packet. Although the total number of grid points is reduced relative to full grid calculations, the solution accuracy can be retained. In particular, since the boundaries adaptively move in response to the dynamics, no advance knowledge of the wave packet evolution is necessary. In addition, no advanced mathematical techniques or complicated computational algorithms are involved in the truncated grid method. Therefore, this method provides a simple and easy-to-implement alternative to propagation of a wave packet in adaptive Eulerian grids. In order to improve the moving boundary truncated grid method, straightforward methods for determining appropriate values of parameters should be developed, and the number of parameters may be reduced. It is also expected that this method can be readily applied to the evolution of classical and quantum phase space ACS Paragon25 Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

distributions.58–62 Furthermore, grid points are stationary in the current truncated grid method. Thus, combination of the truncated grid method and the ALE method associated with adaptive moving grids would be of interest. Additional analysis and applications will be reported in our future studies.

ACKNOWLEDGMENTS

We gratefully acknowledge Ministry of Science and Technology, Taiwan (Grant No. MOST 106-2113-M-007-018) for its financial support of this research.



[email protected]

1

Zhang, J. Z. H. Theory and Application of Quantum Molecular Dynamics; World Scientific: Singapore, 1999.

2

Wyatt, R. E. Quantum Dynamics with Trajectories: Introduction to Quantum Hydrodynamics; Springer: New York, 2005.

3

Lopreore, C. L.; Wyatt, R. E. Quantum Wave Packet Dynamics with Trajectories. Phys. Rev. Lett. 1999, 82, 5190–5193.

4

Wyatt, R. E. Quantum Wave Packet Dynamics with Trajectories: Application to Reactive Scattering. J. Chem. Phys. 1999, 111, 4406–4413.

5

Wyatt, R. E. Quantum Wavepacket Dynamics with Trajectories: Wavefunction Synthesis Along Quantum Paths. Chem. Phys. Lett. 1999, 313, 189–197.

6

Goldfarb, Y.; Degani, I.; Tannor, D. J. Bohmian Mechanics with Complex Action: A New Trajectory-Based Formulation of Quantum Mechanics. J. Chem. Phys. 2006, 125, 231103.

7

Goldfarb, Y.; Schiff, J.; Tannor, D. J. Unified Derivation of Bohmian Mechanics and the Incorporation of Interference Effects. J. Phys. Chem. A 2007, 111, 10416–10421.

8

Goldfarb, Y.; Tannor, D. J. Interference in Bohmian Mechanics with Complex Action. J. Chem. Phys. 2007, 127, 161101.

9

Rowland, B. A.; Wyatt, R. E. Analysis of Barrier Scattering with Real and Complex Quantum Trajectories. J. Phys. Chem. A 2007, 111, 10234–10250.

10

Wyatt, R. E.; Rowland, B. A. Quantum Trajectories in Complex Phase Space: Multidimensional ACS Paragon26 Plus Environment

Page 26 of 31

Page 27 of 31 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

Barrier Transmission. J. Chem. Phys. 2007, 127, 044103. 11

Rowland, B. A.; Wyatt, R. E. Complex Trajectories sans Isochrones: Quantum Barrier Scattering with Rectilinear Constant Velocity Trajectories. J. Chem. Phys. 2007, 127, 164104.

12

Poirier, B. Reconciling Semiclassical and Bohmian Mechanics. V. Wavepacket Dynamics. J. Chem. Phys. 2008, 128, 164115.

13

Poirier, B. Reconciling Semiclassical and Bohmian Mechanics. VI. Multidimensional Dynamics. J. Chem. Phys. 2008, 129, 084103.

14

Park, K.; Poirier, B.; Parlant, G. Quantum Trajectory Calculations for Bipolar Wavepacket Dynamics in One Dimension. J. Chem. Phys. 2008, 129, 194112.

15

Park, K.; Poirier, B. Quantum Trajectory Calculations for Bipolar Wavepacket Dynamics in One Dimension: Synthetic Single-Wavepacket Propagation. J. Theor. Comput. Chem. 2010, 9, 711–734.

16

Trahan, C. J.; Wyatt, R. E. An Arbitrary Lagrangian-Eulerian Approach to Solving the Quantum Hydrodynamic Equations of Motion: Equidistribution with “Smart” Springs. J. Chem. Phys. 2003, 118, 4784–4790.

17

Hughes, K. H.; Wyatt, R. E. Wavepacket Dynamics on Dynamically Adapting Grids: Application of the Equidistribution Principle. Chem. Phys. Lett. 2002, 366, 336–342.

18

Hughes, K. H.; Wyatt, R. E. Wavepacket Dynamics on Arbitrary Lagrangian–Eulerian Grids: Application to an Eckart Barrier. Phys. Chem. Chem. Phys. 2003, 5, 3905–3910.

19

Kendrick, B. K. A New Method for Solving the Quantum Hydrodynamic Equations of Motion. J. Chem. Phys. 2003, 119, 5805–5817.

20

Pauler, D. K.; Kendrick, B. K. A New Method for Solving the Quantum Hydrodynamic Equations of Motion: Application to Two-Dimensional Reactive Scattering. J. Chem. Phys. 2004, 120, 603–611.

21

Kendrick, B. K. Quantum Hydrodynamics: Application to N-Dimensional Reactive Scattering. J. Chem. Phys. 2004, 121, 2471–2482.

22

Derrickson, S. W.; Bittner, E. R.; Kendrick, B. K. Quantum Hydrodynamics: Capturing a Reactive Scattering Resonance. J. Chem. Phys. 2005, 123, 054107.

23

Kendrick, B. K. An Iterative Finite Difference Method for Solving the Quantum Hydrodynamic Equations of Motion. J. Mol. Struct. Theochem 2010, 943, 158–167.

24

Kendrick, B. K. Time-Dependent Wave Packet Propagation Using Quantum Hydrodynamics. ACS Paragon27 Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Theor. Chem. Acc. 2012, 131, 1075. 25

Pettey, L. R.; Wyatt, R. E. Wave Packet Dynamics with Adaptive Grids: The Moving Boundary Truncation Method. Chem. Phys. Lett. 2006, 424, 443 – 448.

26

Pettey, L. R.; Wyatt, R. E. Quantum Wave Packet Dynamics on Multidimensional Adaptive Grids: Applications of the Moving Boundary Truncation Method. Int. J. Quantum Chem. 2007, 107, 1566–1573.

27

Pettey, L. R.; Wyatt, R. E. Application of the Moving Boundary Truncation Method to Reactive Scattering: H + H2 , O + H2 , O + HD. J. Phys. Chem. A 2008, 112, 13335–13342.

28

McCormack, D. A. Dynamical Pruning of Static Localized Basis Sets in Time-Dependent Quantum Dynamics. J. Chem. Phys. 2006, 124, 204101.

29

Hartke, B. Propagation with Distributed Gaussians as a Sparse, Adaptive Basis for HigherDimensional Quantum Dynamics. Phys. Chem. Chem. Phys 2006, 8, 3627–3635.

30

Garashchuk, S. Quantum Trajectory Dynamics in Imaginary Time with the MomentumDependent Quantum Potential. J. Chem. Phys. 2010, 132, 014112.

31

Garashchuk, S.; Mazzuca, J.; Vazhappilly, T. Efficient Quantum Trajectory Representation of Wavefunctions Evolving in Imaginary Time. J. Chem. Phys. 2011, 135, 034104.

32

Garashchuk, S.; Dixit, V.; Gu, B.; Mazzuca, J. The Schr¨ odinger Equation with Friction from the Quantum Trajectory Perspective. J. Chem. Phys. 2013, 138, 054107.

33

Garashchuk, S.; Vazhappilly, T. Multidimensional Quantum Trajectory Dynamics in Imaginary Time with Approximate Quantum Potential. J. Phys. Chem. C 2010, 114, 20595–20602.

34

Garashchuk, S. Calculation of the Zero-Point Energy from Imaginary-Time Quantum Trajectory Dynamics in Cartesian Coordinates. Theor. Chem. Acc. 2012, 131, 1083.

35

Schinke, R.; Nonella, M.; Suter, H. U.; Huber, J. R. Photodissociation of ClNO in the S1 State: A Quantum-Mechanical Ab Initio Study. J. Chem. Phys. 1990, 93, 1098–1106.

36

Untch, A.; Weide, K.; Schinke, R. The Direct Photodissociation of ClNO(S1 ): An Exact ThreeDimensional Wave Packet Analysis. J. Chem. Phys. 1991, 95, 6496–6507.

37

Manthe, U.; Meyer, H.-D.; Cederbaum, L. S. Wave-Packet Dynamics within the Multiconfiguration Hartree Framework: General Aspects and Application to NOCl. J. Chem. Phys. 1992, 97, 3199–3213.

38

Dey, B. K.; Askar, A.; Rabitz, H. Multidimensional Wave Packet Dynamics within the Fluid Dynamical Formulation of the Schr¨ odinger Equation. J. Chem. Phys. 1998, 109, 8770–8782. ACS Paragon28 Plus Environment

Page 28 of 31

Page 29 of 31 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

39

Mayor, F. S.; Askar, A.; Rabitz, H. A. Quantum Fluid Dynamics in the Lagrangian Representation and Applications to Photodissociation Problems. J. Chem. Phys. 1999, 111, 2423–2435.

40

Chou, C.-C. Complex quantum Hamilton-Jacobi Equation with Bohmian Trajectories: Application to the Photodissociation Dynamics of NOCl. J. Chem. Phys. 2014, 140, 104307.

41

Chou, C.-C. Schr¨ odinger-Langevin Equation with Quantum Trajectories for Photodissociation Dynamics. Ann. Phys. (N.Y.) 2017, 377, 22–37.

42

Manthe, U.; Meyer, H.-D.; Cederbaum, L. S. Multiconfigurational Time-Dependent Hartree Study of Complex Dynamics: Photodissociation of NO2 . J. Chem. Phys. 1992, 97, 9062–9071.

43

Kosloff, R.; Tal-Ezer, H. A Direct Relaxation Method for Calculating Eigenfunctions and Eigenvalues of the Schr¨ odinger Equation on a Grid. Chem. Phys. Lett. 1986, 127, 223–230.

44

Roy, A. K.; Gupta, N.; Deb, B. M. Time-Dependent Quantum-Mechanical Calculation of Ground and Excited States of Anharmonic and Double-Well Oscillators. Phys. Rev. A 2001, 65, 012109.

45

Lehtovaara, L.; Toivanen, J.; Eloranta, J. Solution of Time-Independent Schr¨ odinger Equation by the Imaginary Time Propagation Method. J. Comput. Phys. 2007, 221, 148–157.

46

Auer, J.; Krotscheck, E.; Chin, S. A. A Fourth-Order Real-Space Algorithm for Solving Local Schr¨ odinger Equations. J. Chem. Phys 2001, 115, 6841–6846.

47

Aichinger, M.; Krotscheck, E. A Fast Configuration Space Method for Solving Local Kohn–Sham Equations. Comput. Mater. Sci. 2005, 34, 188–212.

48

Chin, S. A.; Janecek, S.; Krotscheck, E. Any Order Imaginary Time Propagation Method for Solving the Schr¨ odinger Equation. Chem. Phys. Lett. 2009, 470, 342–346.

49

Chou, C.-C.; Kouri, D. J. Multidimensional Supersymmetric Quantum Mechanics: A Scalar Hamiltonian Approach to Excited States by the Imaginary Time Propagation Method. J. Phys. Chem. A 2013, 117, 3449–3457.

50

Makri, N.; Miller, W. H. Time-Dependent Self-Consistent Field (TDSCF) Approximation for a Reaction Coordinate Coupled to a Harmonic Bath: Single and Multiple Configuration Treatments. J. Chem. Phys. 1987, 87, 5781–5787.

51

Rom, N.; Moiseyev, N.; Lefebvre, R. Tunneling Rates in a Two-Dimensional Symmetric DoubleWell Potential Surface by the Exterior Scaling Procedure. J. Chem. Phys. 1991, 95, 3562–3569.

52

Gonz´alez-Lezana, T.; Rubayo-Soneira, J.; Miret-Art´es, S.; Gianturco, F. A.; Delgado-Barrio, G.; Villarreal, P. Comparative Configurational Study for He, Ne, and Ar trimers. J. Chem. Phys. ACS Paragon29 Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

1999, 110, 9000–9010. 53

Garashchuk, S.; Light, J. C. Quasirandom Distributed Gaussian Bases for Bound Problems. J. Chem. Phys. 2001, 114, 3929–3939.

54

Feit, M. D.; Fleck Jr., J. A.; Steiger, A. Solution of the Schr¨ odinger Equation by a Spectral Method. J. Comput. Phys. 1982, 47, 412–433.

55

Kosloff, D.; Kosloff, R. A Fourier Method Solution for the Time Dependent Schr¨ odinger Equation as a Tool in Molecular Dynamics. J. Comput. Phys. 1983, 52, 35–53.

56

Kosloff, R.; Kosloff, D. A Fourier Method Solution for the Time Dependent Schr¨ odinger Equation: A Study of the Reaction H+ + H2 , D+ + HD, and D+ + H2 . J. Chem. Phys. 1983, 79, 1823–1833.

57

Kosloff, R. Time-Dependent Quantum-Mechanical Methods for Molecular Dynamics. J. Phys. Chem. 1988, 92, 2087–2100.

58

Donoso, A.; Martens, C. C. Quantum Tunneling Using Entangled Classical Trajectories. Phys. Rev. Lett. 2001, 87, 223202.

59

Donoso, A.; Martens, C. C. Classical Trajectory-Based Approaches to Solving the Quantum Liouville Equation. Int. J. Quantum Chem. 2002, 90, 1348–1360.

60

Trahan, C. J.; Wyatt, R. E. Evolution of Classical and Quantum Phase-Space Distributions: A New Trajectory Approach for Phase Space Hydrodynamics. J. Chem. Phys. 2003, 119, 7017– 7029.

61

Trahan, C. J.; Wyatt, R. E. Classical and Quantum Phase Space Evolution: Fixed-Lattice and Trajectory Solutions. Chem. Phys. Lett. 2004, 385, 280–285.

62

Hughes, K. H.; Wyatt, R. E. Trajectory Approach to Dissipative Quantum Phase Space Dynamics: Application to Barrier scattering. J. Chem. Phys. 2004, 120, 4089–4097.

ACS Paragon30 Plus Environment

Page 30 of 31

Page 31 of 31 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

Graphical TOC Entry

ACS Paragon31 Plus Environment