Modified Gaussian Wave Packet Method for Calculating Initial State

Subscriber access provided by UNIV OF NEW ENGLAND ARMIDALE

Dynamics

A modified Gaussian Wave Packet Method for Calculating Initial State Wavefunctions in Photodissociation Shanyu Han, Daiqian Xie, and Hua Guo J. Chem. Theory Comput., Just Accepted Manuscript • DOI: 10.1021/acs.jctc.8b00771 • Publication Date (Web): 20 Sep 2018 Downloaded from http://pubs.acs.org on September 25, 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 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 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.

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 39 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

Journal of Chemical Theory and Computation

ct-2018-007714, revised 9/18/2018

An Modified Gaussian Wave Packet Method for Calculating Initial State Wavefunctions in Photodissociation

Shanyu Han,1,2 Daiqian Xie,1,* and Hua Guo2,* 1

Institute of Theoretical and Computational Chemistry, Key Laboratory of

Mesoscopic Chemistry, School of Chemistry and Chemical Engineering, Nanjing University, Nanjing 210093, China 2

Department of Chemistry and Chemical Biology, University of New Mexico, Albuquerque, New Mexico 87131, USA

* Corresponding authors: [email protected] and [email protected]

1

ACS Paragon Plus Environment

Journal of Chemical Theory and Computation 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

Abstract A modified Gaussian wave packet relaxation method is proposed to calculate the ground state wavefunction using an expansion of frozen Gaussian wave packets. This new procedure consists of two steps. In the first step, a multi-dimensional Gaussian product placed at the ground state equilibrium geometry is propagated in imaginary time. The relaxation optimizes the widths of the one-dimensional Gaussians. In the second step, additional Gaussian wave packets with the same widths are placed near the equilibrium geometry, and the corresponding expansion coefficients are optimized using the same relaxation method. This new algorithm is tested in photodissociation of NOCl and NH3, and the results show good agreement with the exact results in the energy, wavefunction, and absorption spectrum. In particular, the highly structured absorption spectrum of NH3 is reproduced, underscoring the accuracy of both the initial wave packet and excited state propagation.

2


Page 2 of 39

Page 3 of 39 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 Despite the tremendous success of exact grid-based wave packet methods in

characterizing quantum dynamics of small molecular systems,1-6 it is difficult to extend them to larger systems because of their exponential scaling with dimensionality.7 Developing methods capable of treating high-dimensional systems is a main goal of computational chemistry and physics. Gaussian wave packet (GWP) based methods are capable of providing such an efficient tool for treating quantum dynamics, albeit approximately, with many degrees of freedom (DOFs). The idea of using Gaussian functions as a basis set to solve the nuclear time-dependent Schrödinger equation can be traced back to the pioneering work of Heller.8-9 He demonstrated that a superposition of simple frozen Gaussians, the widths of which are fixed, is an useful and general approach for time-dependent quantum dynamics.9 Many GWP-based methods using frozen Gaussians have since been developed. To name a few, Martínez and coworkers developed the algorithm of multiple spawning, which allows new Gaussians to be added during non-adiabatic transitions.10-11 This algorithm has been successfully applied to direct non-adiabatic dynamics, known as Ab Initio Multiple Spawning (AIMS).12-13 The Coupled Coherent States (CCS) method developed by Shalashilin and co-workers14 was extended to multistate non-adiabatic dynamics, known as Multi-Configurational Ehrenfest (MCE).15 The Ab Initio Multiple Cloning (AIMC) algorithm16 is a combination of MCE and AIMS. Burghardt17 introduced a Gaussian (G) wave packet method within the Multi-Configuration Time-Dependent Hartree (MCTDH) framework.18-19 This 3



Page 4 of 39

so-called G-MCTDH method has been successfully applied to system-bath problems, in which the bath modes are approximated by Gaussians while the system is characterized by the more accurate MCTDH.20-21 As a simplification of the G-MCTDH method, one can use an all-Gaussian basis, which is known as variational multi-configuration Gaussian (vMCG) method.22 This method has been used for vibronic coupled systems treating non-adiabatic dynamics with models21 as well as direct dynamics.22-25 These GWP methods can be divided into two categories, depending on the equations of motion (EOMs) for the Gaussian center and momentum. The methods that follow classical trajectories are called classical Gaussians. The force for the Gaussian comes from either a single Born-Oppenheimer state or an Ehrenfest average of multiple states. On the other hand, the G-MCTDH/vMCG method uses a more optimal equation derived by applying the Dirac-Frenkel variational principle.17, 20 The variational GWP approach has been proved to be more accurate than classical Gaussian methods because it introduces quantum correlations.21, 26 For instance, it has been demonstrated that the G-MCTDH/vMCG propagation can reproduce more detailed autocorrelation functions within the first few recurrences of the autocorrelation function.21 Here, we are primarily interested in photochemical processes, in which the initial wave packet on an excited state is taken within the Condon approximation as a vibrational

wavefunction

on

the

ground

electronic

state.27

The

current

G-MCTDH/vMCG implementation relies on a very approximate choice of the initial 4


Page 5 of 39 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


wave packet, namely a single GWP in multiple dimensions. Such a choice is reasonable for problems formulated in normal coordinates, but may lead to significant errors when other coordinates are used. Under such circumstances, the simulation results are not quantitatively accurate because of the errors in the initial wave packet. In this publication, we discuss a scheme that allows a significant improvement in the definition of the initial wave packet within the G-MCTDH/vMCG framework. We demonstrate that the photo-absorption spectrum can be reproduced with reasonable accuracy for both direct and indirect dissociation processes, thanks to a better representation of the initial wave packet. II.

Theory Our improvement is based on the premise that a better approximation of the

initial wave packet can be achieved by expanding it in terms of multiple, rather than a single, GWPs. Such a strategy has been attempted before in AIMS12 and matching pursuit,28 but how to determine the expansion coefficients within the framework of a GWP method remains an open question. The widths of the Gaussians were often not optimally determined in the G-MCTDH/vMCG and most other GWPs methods, although attempts have been made. Martínez and coworkers, for example, proposed a way to optimize Gaussian widths in Cartesian coordinates by maximizing the overlap with the normal-mode eigenstate,29 which is itself a Gaussian function for the ground state of an harmonic oscillator with a width of

2 / 2 if mass-weighted coordinate is

used. Worth and coworkers also proposed an alternative method based on minimizing the energy within the MCTDH framework, which is more general, especially for 5



anharmonic eigenstates, which do not have the Gaussian form.30 Their approach is based on relaxation, namely propagating the wavefunction in imaginary time.31 This method has been implemented within the MCTDH framework,32-34 and its G-MCTDH/vMCG version is implemented in Quantics.35 However, the current implementation of the relaxation method in Quantics is not sufficiently robust and general. In the variable mean-field integration scheme.36 both the expansion coefficients and Gaussian parameters, including the widths, centers, and momenta, are relaxed simultaneously in every step.21 In this work, a modified relaxation method within the G-MCTDH/vMCG framework is proposed. Our approach is quite straightforward: a single multi-dimensional Gaussian is first placed at the potential minimum of the ground electronic state potential energy surface (PES) and the width in each dimension is optimized by relaxation. Then, additional GWPs with the same widths as determined from the first step are added to the expansion, and their expansion coefficients are optimized again using relaxation. This method is motivated by the consideration that the optimization of the nonlinear parameter (widths) is better performed separately from that of the linear parameters (expansion coefficients). It is tested in a triatomic system NOCl and a tetratomic system NH3. The calculated ZPEs and ground state wavefunctions are found to converge reasonably well to the exact ones. The wavefunctions were then used as the initial wave packets to propagate in real time on the excited state PES. The absorption spectrum is thus obtained via Fourier transform of the real time autocorrelation function.37 The NOCl results show a marked improvement over the previous one 6


Page 6 of 39

Page 7 of 39 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


without relaxation.24 More remarkably, the exact absorption spectrum of NH3, which has multiple peaks due to the indirect nature of the dissociation, is reproduced very well. a. G-MCTDH/vMCG relaxation Since the basic theory of G-MCTDH/vMCG and corresponding working equations have been extensively discussed in Refs. 17, 21, 26, we only present the relevant equations for relaxation. Following the notions in a recent review,26 a time-dependent wavefunction of N DOFs can in general be expanded as a linear combination of N-dimensional GWPs ( g j (x, t ) ): n

Ψ ( x, t ) = ∑ A j (t ) g j ( x, t ) ,

(1)

j =1

in which Aj (t ) are the corresponding time-dependent coefficients and x denotes the N-dimensional coordinates to be defined. The GWPs, g j (x, t ) , are defined in matrix notion as follows:

g j (x, t ) = exp ( xT ⋅ ς j ( t ) ⋅ x + ξ j ( t ) ⋅ x + η j ( t ) ) , where the superscript

T

denotes the transpose.

(2)

Λ j = {ς j , ξ j ,η j }

are the

time-dependent Gaussian parameters. ς j is the N × N width matrix that couples all DOFs, ξ j is a N-dimensional linear parameter vector, and η j is a scalar. A GWP with a ς j matrix allowing for coupling between different DOFs through off-diagonal terms is called a thawed Gaussian. The width matrix ς j of a Gaussian can be chosen to be diagonal, which is called a separable Gaussian or frozen Gaussian depending on whether it is time-dependent or time-independent. Assuming no coupling between

7



Page 8 of 39

different DOFs in the width matrix, g j (x, t ) can be given in a product of one-dimensional GWPs as:

  g j (x, t ) = ∏ g j ( xκ , t ) = exp  ∑ (ς jκ xκ 2 + ξ jκ xκ ) + η j   κ  κ

(3)

with κ running over the system DOFs represented by x. To the best of our knowledge, the equations of motion (EOMs) for the G-MCTDH/vMCG relaxation method has not been explicitly presented in any publication though already coded in Quantics.35 Here, we give the explicit form for propagating the wave packet in imaginary time. The EOMs for the Aj (t ) coefficients and Gaussian parameters are written as follows:17 −1   A& j (t ) = −  ∑  S jl  ( H lm − iτ lm ) -E1  Am (t ) ,  lm 

(4a)

& = − [C]−1 Y , Λ

(4b)

where S is the overlap matrix for two GWPs, the matrix elements of which are given by:

S jl = g j g l

(5)

and H is the Hamiltonian matrix with matrix elements defined as H jl = g j Hˆ g l .

Finally,

(6)

τ jl is the restriction matrix

τ jl = g j g& l = ∑ g j α

∂g l & λlα , ∂λlα

(7)

where λlα ∈Λ denotes the Gaussian parameters of gl with index

α and λ&lα

denotes its time derivative. For a separable Gaussian with N DOFs, the index

8


α

Page 9 of 39 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


denotes all time-dependent parameters including one scalar parameter ηl and N linear parameters

ξlκ , and N widths. Thus, λlα is explicitly given as

λl ,α =1 ≡ ηl , (8a)

λl ,α >1 ≡ ξl ,κ =α −1 , α = 1 ~ ( N + 1) , (8b)

λl ,α > N +1 ≡ ς l ,κ =α − N −1 , α = 1 ~ (2 N + 1) .

(8c)

For frozen Gaussian, the widths are fixed. Finally,

∑∑ A ⋅ H E= ∑∑ A ⋅ S j

j

* j

jl

⋅ Al

* j

jl

⋅ Al

l

is the

l

expectation value of the Hamiltonian. For Eq. (4b), the matrices C and Y are given by:

(

C jα ,l β = ρ jl G (jlαβ ) − G (α 0)S -1G (0β )  Y jα =

∑ρ

jl

l

(H

(α 0 ) jl

−  G (α 0 ) S -1 H 

jl

jl

),

(9a)

),

(9b)

where ρ jl = Aj Al are the density matrix elements and the definition of other matrices *

are given below:

G(jlαβ ) = G(jlα 0) = H (jlα 0) =

∂g j ∂gl , ∂λjα ∂λlβ ∂g j ∂λ jα

(10a)

gl ,

(10b)

∂g j ˆ H gl . ∂λ jα

(10c)

b. Modified Gaussian relaxation In our approach, we carry out the first step of the relaxation using a single separable multi-dimensional GWP, denoted as g 0 and its widths are optimized using 9



Page 10 of 39

Eq. (4b). The center of the GWP is set at the potential minimum of the ground state PES. The matrix elements of C and Y can be simplified as follows: (α 0) (0β ) Cα , β = ( G00(αβ ) − G00 G00 ) ,

(11a)

Yα = ( H 00(α 0) − G00(α 0) H 00 ) ,

(11b)

where the density and overlap matrices can be omitted as only one single GWP is used and the subscript index 0 in matrix denotes g 0 . In the second step, more GWPs with the same widths as g 0 are added following the scheme used in G-MCTDH/vMCG.24 The first group of these so-called excitation GWPs are placed around g 0 and their positions are chosen such as they have an overlap of 0.8 with g 0 .21 This is obtained by shifting one of the multidimensional centers with a placement according to the overlap while other dimension coordinates are kept unchanged. Once centers of all DOFs are all filled by new GWPs, a second group of GWPs follows with an overlap of 0.64(0.8 × 0.8) with g 0 and so on. Such a choice is completely empirical. Too small an overlap might let to incompleteness of the basis, while too large an overlap could cause numerical instability due to near-zero eigenvalues of the overlap matrix. The initial coefficients of g 0 and the excitation GWPs are set 1.0 and 0.0, respectively; and they evolve in imaginary time by solving the following equation: −1   A& j = −  ∑  S jl  Hlm -E1  Am .  lm 

(12)

The Hamiltonian matrix elements ( Hlm ) are evaluated analytically in terms of Gaussian moments,21 where the local harmonic approximation (LHA) is used for evaluating the potential matrix. This approximation has been used by Heller8, 38 and 10


Page 11 of 39 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


also in G-MCTDH/vMCG,17 in which the PES is expanded quadratically with respect to the time-dependent Gaussian center. As the Gaussian parameters in this step are time-independent so that the restriction matrix

τ jl (Eq. (7)) vanishes in Eq. (4). The

well-known 5th order Runge-Kutta scheme is employed to solve Eq. (4). The energy expectation value, namely the zero-point energy (ZPE) because only the ground vibrational state is considered here, is used as an indicator of convergence. c. Coordinate systems and bases Nuclear motions of a molecule with N atoms in laboratory frame (LF) could be described by a center of mass vector R cm and N-1 relative position vectors,39 thus overall translational motion is only related to R cm and naturally removed. As implemented in Quantics, the current G-MCTDH/vMCG method can only treat rectilinear coordinate systems. To better represent the systems under investigation, we chose Jacobi vectors as they are orthogonal. Using the vibrational kinetic energy operator of a tri-atom system derived in previous work,24 we investigate the classic example of NOCl photodissociation32, 40-43 explicitly in a so-called Jacobi Cartesian coordinate system. Following the same idea, we extended the model to the tetra-atomic system of NH3, which is a prototype for photodissociation of tetra-atomic systems.44-52 i.

NOCl Two independent Jacobi vectors can be used to describe the dynamics of the three

r uuur atoms in NOCl: r = NO describes the vibrational motion of the NO diatom and

11



Page 12 of 39

r uuur R = gCl describes the Cl dissociative motion relative to the mass center, denoted by g, of the NO diatom. The corresponding reduced masses are given below: m=

mN mO , and mN + mO

M =

( mN + mO ) mCl m N + mO + mCl

(13)

The Cartesian coordinates of the two vectors, ( xr , yr , zr ) and ( xR , yR , zR ), depend on how we define the BF frame. Following the same definition of the BF and Cartesian coordinates24 known as two-vector embedded reference frame,39,

53

the

three-dimension vibrational motion can be only related to zR , xr and zr , while others remain to be zero values. So the kinetic energy operator is given as

1  Pˆ 2 pˆ 2 pˆ 2  Tˆvib =  Z + x + z  2 M m m 

(14)

with h ∂ ˆ h ∂ ˆ h ∂ , Px ≡ . PˆZ ≡ Pz ≡ i ∂z R i ∂xr i ∂zr

(15)

Then the Gaussian basis defined in this coordinate is written as:

g ( zR ,xr ,zr ) = g ( zR ) g ( xr ) g ( zr ) ,

(16)

  q − q 2  0 g (q ) = N q exp  −    , q = zR , xr , zr ,   2σ q    

(17)

where Nq is the normalization factor, q0 is the Gaussian center of each dimension and σ q is the corresponding width. ii.

NH3 There are several possible choices of Jacobi vectors for a tetra-atomic molecule,

which may not be equally efficient for every system. A proper choice may lead to better description of the vibration and thus fast convergence. In Ref. 54, Gatti provided 12


Page 13 of 39 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 choice of defining the three Jacobi vectors of ammonia for studying its umbrella motion, which represents the initial reaction coordinate in ammonia photodissociation. To this end, the three Jacobi vectors are given as:

r uuuuuuuur r uuuuuuur r uuuuur R1 = H (1) H (2) , R2 = G12 H (3) , R3 = G123 N ,

(18)

where G12 and G123 are the mass centers of the H (1) H (2) and H (1) H (2) H (3) moieties, respectively. Similarly, the two-vector embedding scheme39 is used to define the BF frame.

r

r

The z-axis of BF is defined paralleled to R2 . With the help of R1 , the BF is defined as:

r R2 r ez = r , R2

r R2 r ey = r ∧ R2

r R1 r r r r , and ex = ey ∧ ez . R1

(19)

r Then the Cartesian coordinates of each vector are given as Ri =1,2,3 = ( xi , yi , zi ) . Similar

restrictions in NOCl were used to exclude the overall rotation, which is written as

y1 = y2 = x2 = 0 . So we can derive the 6-D vibrational kinetic energy operator:  ∂2 ∂2 ∂2 ∂2 ∂2 ∂2  , Tˆ = − h 2  + + + + + 2 2 2 2 2 2   2 µ1∂z1 2 µ1∂x1 2 µ2 ∂z2 2 µ3∂x3 2 µ3∂y3 2 µ3∂z3 

(20)

in which µi =1,2,3 is the corresponding reduced mass of each Jacobi vector. Likewise, the Gaussian basis is given as:

g (q) = ∏ g ( qi ) , with qi = z1 , x1 , z2 , x3 , y3 , z3 . i

III.

Results

a. Zero-point energies and ground state wavefunctions i. NOCl

13


(21)


The ground state PES of Manthe was used.32 In the first step of optimization, we tested different sets of the initial guess widths ( σ Z , σ x , σ z ) for g 0 (shown in Table 1) and all were found to be converged to the same final values after relaxation, as shown in Figure 1, even when some initial guesses were quite far from the destination, e.g., 0.5 compared to the final value of 0.0791. In the second step, the wavefunction is expanded in terms of the g0 and additional excitation GWPs. The expansion coefficients are obtained by solving the coefficient relaxation Eq. (4a) using the widths optimized in the first step. The evolution of the energy expectation value is shown in Figure 2 in imaginary time for different numbers of GWPs. Increasing the number of GWPs eventually leads to convergence of the ZPE, which becomes quite close to the ZPE value obtained using Lanczos algorithm, as shown in Table 2. To illustrate the better convergence behavior of the new method, we also computed the results with fixed widths. As shown in Table 2 (Column 3), the convergence of the latter is inferior. This result suggest both the widths and expansion coefficients have a large impact on the convergence of the ground state wavefunction. To further illustrate the convergence, we calculated the overlap between the exact wavefunction and that obtained using the modified G-MCTDH/vMCG method. The exact wavefunction was obtained by diagonalizing the Hamiltonian with the Lanczos algorithm. As shown in Table 2 (Column 4), the results are quite satisfactory, especially given the use of LHA, which introduces errors to the Hamiltonian. ii.

NH3 14


Page 14 of 39

Page 15 of 39 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 calculation was performed on the ground (X) adiabatic PES obtained from a set of recent developed diabatic potential energy matrix by Yarkony and coworkers.48-49 As in the case of NOCl described above, the widths of the six DOFs in the single Gaussian (g0) are first optimized and the converged values are listed in Table 3. With the parameters in Table 3, the second step relaxation is performed to optimize the expansion coefficients. As expected, adding excitation GWPs helps to improve the agreement with both the ZPE and wavefunction overlaps, as shown in Table 4. The calculated ZPE is lower than the exact one by as much as 0.1 eV, which is not small. Similar problems have been noted in previous calculations of energy levels of water and glycine by Fourier transform of the real time autocorrelation function which is obtained by propagation of coefficients of time-independent Gaussian function bases.55 We believe that these errors can be attributed to the LHA used in the calculations. Despite the discrepancy of ZPE, the overlaps listed in the table suggest that the wavefunction is quite sufficient. This point is confirmed below. b. Absorption spectra i.

NOCl The ground state wavefunction determined in the previous subsection is used to

compute the absorption spectrum using the excited state Hamiltonian. The excited state PES of Schinke and coworkers was used.40 The resulting spectrum is compared with the exact MCTDH result and the G-MCTDH/vMCG ones with fixed Gaussian

15



widths from Ref.

24

Page 16 of 39

. In the latter calculations, multiple Gaussians with the fixed

widths ( σ Z , σ x , and σ z = 0.0794, 0.1142, and 0.1370, respectively) were used. As can be seen from the Figure 3, the result without relaxation (red line) has unphysical oscillations in the blue wing of the absorption spectrum, which is believed to be caused by the poor description of the initial state.24 The oscillations are smoothened when the widths and expansion coefficients of the GWPs are relaxed as they provide a better description of the ground state wavefunction. The peak energy of the absorption spectrum is compared in Table 5, which shows that relaxation also improves its agreement with the exact results. The exact absorption spectrum also has a small bump near 1.3 eV, which was attributed to a recurrence of the autocorrelation function due to fast N-O vibrational motion on the excited state PES.41 This bump is also seen in the absorption spectrum obtained using the G-MCTDH/vMCG method with relaxation, but the amplitude is much smaller. This quantitative difference is attributed to the error introduced by the LHA on the excited state. ii.

NH3 We used 40, 60, and 80 GWPs to calculate the absorption spectrum of ammonia.

The propagation was done on the excited (A) adiabatic state PES. A window function was

used

to

alleviate

( )

the

Gibbs

phenomenon,

which

is

written

as

( ) is the Heaviside function. For the results t

g (t )= cos( 2πTt )Θ 1 − T , where Θ 1 − T t

reported here, the following values are used: T=50 fs, forcing the autocorrelation

(

function to be 0 at T=50 fs. An additional damping exp − ( t / τ ) 16


2

)

was used which

Page 17 of 39 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


is also used in MCTDH calculations.7 Details about the choice of damping function and parameter τ can be found in Ref. 7. The result is compared with the recent work using a grid based method48 and same window and damping function for comparison, which serves as benchmark. As shown in Figure 4, the overall shape of the absorption spectrum obtained using the modified G-MCTDH/vMCG method agrees well with the benchmark despite some discrepancies on low energy peaks, which might be attributed to errors either the initial state and/or propagation. The individual peaks in the absorption spectrum can be assigned to the umbrella (v2) mode of NH3, which is excited due to the pyramidal-to-planar transition.48 These features are labeled as 2n, where n denotes the corresponding quantum number. Due to the double-well nature of the ground state PES, the odd and even quantum numbers correspond to different inversion symmetries. In our G-MCTDH/vMCG calculations, however, the wavefunction is placed in one of the two equivalent wells, thus ignoring the inversion symmetry. The net consequence is the neglect of the tunneling splitting, but it has a small impact on the absorption spectrum. The diffuse peaks in the absorption spectrum suggest strong recurrences of the autocorrelation

function.

This is

very

difficult

to

reproduce within

the

G-MCTDH/vMCG formalism because long-time propagation on the excited state PES typically distorts a wave packet, making a Gaussian representation inaccurate. In the current case, the reproduction of the absorption spectrum is due largely to the fact that the reaction coordinate in the Franck-Condon region is mainly along the umbrella 17



coordinate, which is well represented in our model. This underscores the importance of coordinates in excited state dynamics. IV.

Conclusions In this work, we propose a modified G-MCTDH/vMCG relaxation method for

calculating the ground state wavefunction, used for photodissociation dynamics calculations. To this end, the wavefunction is expanded in terms of frozen GWPs. The Gaussian widths and expansion coefficients are determined by propagation in imaginary time, which minimizes energy. The optimization of both parameters helps to yield more accurate results. Most previous studies all used approximations or other methods to estimate these important parameters. To the best of our knowledge, it is the first time that both the Gaussian widths and expansion coefficients are optimized using rigorous equations. Differing from previous methods, however, the optimization is achieved in two separate steps, in order to avoid optimizing the linear (expansion coefficients) and non-linear (widths) parameters at the same time. Judging from both the energy expectation value and overlap with the exact wavefunction, this new approach is capable of converging the results to a satisfactory level. It is also robust, as evidenced by rapid convergence regardless of initial guesses. This method is demonstrated in two cases. The first example is the direct photodissociation of NOCl, involving three DOFs. The second example is the indirect photodissociation of NH3, which has six DOFs. In both cases, we show that the method described in this work yields good agreement with the exact results. The reproduction of the highly-structured spectrum of NH3 is particularly satisfactory 18


Page 18 of 39

Page 19 of 39 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


because it implies the Gaussian wave packet is capable of handling relatively long dynamics. We believe that this method can certainly improve the performance of G-MCTDH/vMCG dynamics calculations and extendable to nonadiabatic dynamics involving more than one electronic state.

Acknowledgements: S.H. and D.X. thank the National Natural Science Foundation of China (Grant Nos. 21733006 and 21590802 to D.X.) and the National Key Research and Development Program of China (2017YFA0206500 to D.X.). H.G. acknowledges the US Department of Energy (DE-SC0015997) for generous support. We sincerely thank Prof. Graham Worth for several useful discussions.

19



Page 20 of 39

Table 1. Initial guesses of the width (Bohr, a0) parameters and energy expectation values for the NOCl system in comparison with their final optimized values.

Initial Final

σ Z (a0)

σ x (a0)

σ z (a0)

E (eV)

0.5000 0.3000 0.1000 0.0500 0.0791

0.5000 0.3000 0.1000 0.0500 0.0755

0.5000 0.3000 0.1000 0.0500 0.0836

3.9366 0.6453 0.2194 0.2860 0.1966

20


Page 21 of 39 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


Table 2. Convergence of the ZPE with respect to the number of GWPs for the NOCl system. The convergence is also demonstrated by the overlap with the exact wavefunction and the exact value of ZPE obtained using the Lanczos algorithm. No. of GWPs 1 7 13 21 30 40 60 70 80 Lanczos algorithm a.

Optimal widths were used.

b.

All three widths

ZPE (eV)a 0.1965 0.1963 0.1768 0.1705 0.1704 0.1700 0.1693 0.1686 0.1685 0.1688

σ Z , σ x , and σ z

ZPE (eV)b 0.2194 0.2016 0.1866 0.1779 0.1773 0.1745 0.1731 0.1726 0.1725

were set to 0.1.

21


Overlapa 0.8769 0.8777 0.9350 0.9855 0.9857 0.9880 0.9919 0.9963 0.9977


Page 22 of 39

Table 3. Center positions of the six 1D GWPs and their optimized widths (in Bohr, a0) for the NH3 system.

Coordinate Center Width

z2 2.6604 0.1736

x1 3.0719 0.2004

z1 0.0000 0.2527

x3 0.0000 0.0967

22


z3 0.0000 0.0968

y3 0.7323 0.1259

Page 23 of 39 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


Table 4. Energy expectation values and overlap between the G-MCTDH/vMCG wavefunction with the exact one obtained using the Lanczos algorithm for the NH3 system. No. of GWPs 20 40

ZPE (eV) 0.8888 0.8685

Overlap 0.8836 0.9165

60 80

0.8522 0.8515

0.9293 0.9296

120

0.8506

0.9306

200

0.8486

0.9310

Exact

0.9525

1.0

23



Table 5. Comparison of the peak energy of the absorption spectrum calculated using G-MCTDH/vMCG with and without relaxation for the NOCl system. The MCTDH result is included as the benchmark. No. of GWPs 13 21 30 40 MCTDH

E in eV (∆E in cm-1) without relaxation 1.1488 (+135.5) 1.1600 (+225.8) 1.1572 (+203.3) 1.1516 (+158.1) 1.1320

24


E in eV (∆E in cm-1) with relaxation 1.1322 (+1.6) 1.1376 (+45.2) 1.1432 (+90.3) 1.1404 (+67.8)

Page 24 of 39

Page 25 of 39 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


References: 1.

Kosloff, R., Time-dependent quantum-mechanical methods for molecular dynamics. J. Phys.

Chem. 1988, 92, 2087-2100. 2.

Althorpe, S. C.; Clary, D. C., Quantum scattering calculations on chemical reactions. Annu. Rev.

Phys. Chem. 2003, 54, 493-529. 3.

Balint-Kurti, G. G., Time-dependent and time-independent wavepacket approaches to reactive

scattering and photodissociation dynamics. Int. Rev. Phys. Chem. 2008, 27, 507-539. 4.

Nyman, G.; Yu, H.-G., Quantum approaches to polyatomic reaction dynamics. Int. Rev. Phys.

Chem. 2013, 32, 39-95. 5.

Zhang, D. H.; Guo, H., Recent advances in quantum dynamics of bimolecular reactions. Annu.

Rev. Phys. Chem. 2016, 67, 135-158. 6.

Hu, X.; Zhou, L.; Xie, D., State-to-state photodissociation dynamics of the water molecule.

WIREs Comput. Mol. Sci. 2017, e1350. 7.

Meyer, H.-D.; Gatti, F.; Worth, G. A., Multidimensional Quantum Dynamics: MCTDH Theory and

Applications. Wiley-VCH: 2009. 8.

Heller, E. J., Time-dependent approach to semiclassical dynamics. J. Chem. Phys. 1975, 62,

1544-1555. 9.

Heller, E. J., Frozen Gaussians: A very simple semiclassical approximation. J. Chem. Phys. 1981,

75, 2923-2931. 10. Martinez, T. J.; Ben-Nun, M.; Levine, R. D., Multi-electronic-state molecular dynamics: A wave function approach with applications. J. Phys. Chem. 1996, 100, 7884-7895. 11. Ben-Nun, M.; Martı́nez, T. J., Nonadiabatic molecular dynamics: Validation of the multiple spawning method for a multidimensional problem. J. Chem. Phys. 1998, 108, 7244-7257. 12. Ben-Nun, M.; Martinez, T. J., Ab initio quantum molecular dynamics. Adv. Chem. Phys. 2002, 121, 439-512. 13. Curchod, B. F. E.; Martínez, T. J., Ab initio nonadiabatic quantum molecular dynamics. Chem. Rev. 2018, 118, 3305-3336. 14. Shalashilin, D. V.; Child, M. S., The phase space CCS approach to quantum and semiclassical molecular dynamics for high-dimensional systems. Chem. Phys. 2004, 304, 103-120. 15. Shalashilin, D. V., Quantum mechanics with the basis set guided by Ehrenfest trajectories: Theory and application to spin-boson model. J. Chem. Phys. 2009, 130, 244101. 16. Makhov, D. V.; Glover, W. J.; Martinez, T. J.; Shalashilin, D. V., Ab initio multiple cloning algorithm for quantum nonadiabatic molecular dynamics. J. Chem. Phys. 2014, 141, 054110. 17. 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. 18. Meyer, H.-D.; Manthe, U.; Cederbaum, L. S., The multi-configuration time-dependent Hartree approach. Chem. Phys. Lett. 1990, 165, 73-78. 19. Manthe, U.; Meyer, H.-D.; Cederbaum, L. S., Multiconfiguration time-dependent Hartree study of complex dynamics: Photodissociation of NO2. J. Chem. Phys. 1992, 97, 9062-9071. 20. 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. 25



21. Burghardt, I.; Giri, K.; Worth, G. A., Multimode quantum dynamics using Gaussian wavepackets: The Gaussian-based multiconfiguration time-dependent Hartree (G-MCTDH) method applied to the absorption spectrum of pyrazine. J. Chem. Phys. 2008, 129, 174104. 22. Worth, G. A.; Robb, M. A.; Burghardt, I., A novel algorithm for non-adiabatic direct dynamics using variational Gaussian wavepackets. Faraday Disc. 2004, 127, 307-323. 23. Worth, G. A.; Robb, M. A.; Lasorne, B., Solving the time-dependent Schrödinger equation for nuclear motion in one step: direct dynamics of non-adiabatic systems. Mol. Phys. 2008, 106, 2077-2091. 24. Lasorne, B.; Robb, M. A.; Worth, G. A., Direct quantum dynamics using variational multi-configuration Gaussian wavepackets. Implementation details and test case. Phys. Chem. Chem. Phys. 2007, 9, 3210-3227. 25. Mendive-Tapia, D.; Lasorne, B.; Worth, G. A.; Robb, M. A.; Bearpark, M. J., Towards converging non-adiabatic direct dynamics calculations using frozen-width variational Gaussian product basis functions. J. Chem. Phys. 2012, 137, 22A548. 26. Richings, G. W.; Polyak, I.; Spinlove, K. E.; Worth, G. A.; Burghardt, I.; Lasorne, B., Quantum dynamics simulations using Gaussian wavepackets: the vMCG method. Int. Rev. Phys. Chem. 2015, 34, 269-308. 27. Schinke, R., Photodissociation Dynamics. Cambridge University Press: Cambridge, 1993. 28. Wu, Y.; Batista, V. S., Matching-pursuit for simulations of quantum processes. J. Chem. Phys. 2003, 118, 6720-6724. 29. Thompson, A. L.; Punwong, C.; Martínez, T. J., Optimization of width parameters for quantum dynamics with frozen Gaussian basis sets. Chem. Phys. 2010, 370, 70-77. 30. Polyak, I.; Allan, C. S. M.; Worth, G. A., A complete description of tunnelling using direct quantum dynamics simulation: Salicylaldimine proton transfer. J. Chem. Phys. 2015, 143, 084121. 31. Kosloff, R.; Tal-Ezer, H., A direct relaxation method for calculating eigenfunctions and eigenvalues of the Schroedinger equation on a grid. Chem. Phys. Lett. 1986, 127, 223-230. 32. 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. 33. Meyer, H.-D.; Worth, G. A., Quantum molecular dynamics: propagating wavepackets and density operators using the multiconfiguration time-dependent Hartree method. Theo. Chem. Acc. 2003, 109, 251-267. 34. Meyer, H.-D.; Quéré, F. L.; Léonard, C.; Gatti, F., Calculation and selective population of vibrational levels with the Multiconfiguration Time-Dependent Hartree (MCTDH) algorithm. Chem. Phys. 2006, 329, 179-192. 35. Worth, G. A.; Giri, K.; Richings, G. W.; Burghardt, I.; Beck, M. H.; Jäckle, A.; Meyer, H.-D. The QUANTICS Package, Version 1.1, University of Birmingham, Birmingham, U.K, 2015. 36. Beck, M. H.; Jackle, A.; Worth, G. A.; Meyer, H.-D., The multiconfiguration time-dependent Hartree (MCTDH) method: A highly efficient algorithm for propagating wavepackets. Phys. Rep. 2000, 324, 1-105. 37. Balint-Kurti, G. G.; Dixon, R. N.; Marston, C. C., Time-dependent quantum dynamics of molecular photofragmentation processes. J. Chem. Soc. Faraday Trans. 1990, 86, 1741-1749. 38. Lee, S.-Y.; Heller, E. J., Exact time-dependent wave packet propagation: Application to the photodissociation of methyl iodide. J. Chem. Phys. 1982, 76, 3035-3044. 26


Page 26 of 39

Page 27 of 39 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


39. Mladenović, M., Rovibrational Hamiltonians for general polyatomic molecules in spherical polar parameterization. I. Orthogonal representations. J. Chem. Phys. 2000, 112, 1070-1081. 40. 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. 41. Untch, A.; Weide, K.; Schinke, R., The direct photodissociation of ClNO(S1): An exact three-dimensional wave packet analysis. J. Chem. Phys. 1991, 95, 6496-6507. 42. Guo, H.; Seideman, T., Quantum mechanical study of photodissociation of oriented ClNO(S1). Phys. Chem. Chem. Phys. 1999, 1, 1265-1272. 43. van Harrevelt, R.; Manthe, U., Multidimensional time-dependent discrete variable representations in multiconfiguration Hartree calculations. J. Chem. Phys. 2005, 123, 064106. 44. Bonhommeau, D.; Valero, R.; Truhlar, D. G.; Jasper, A. W., Coupled-surface investigation of the photodissociation of NH3(Ã): Effect of exciting the symmetric and antisymmetric stretching modes. J. Chem. Phys. 2009, 130, 234303. 45. Lai, W.; Lin, S. Y.; Xie, D.; Guo, H., Full-dimensional quantum dynamics of A-state photodissociation of ammonia. Absorption spectra. J. Chem. Phys. 2008, 129, 154311. 46. Lai, W.; Lin, S. Y.; Xie, D.; Guo, H., Non-adiabatic dynamics of A-state photodissociation of ammonia: a four-dimensional wave packet study. J. Phys. Chem. A 2010, 114, 3121-3126. 47. Giri, K.; Chapman, E.; Sanz, C. S.; Worth, G., A full-dimensional coupled-surface study of the photodissociation dynamics of ammonia using the multiconfiguration time-dependent Hartree method. J. Chem. Phys. 2011, 135, 044311. 48. Zhu, X.; Ma, J.; Yarkony, D. R.; Guo, H., Computational determination of the Ã state absorption spectrum of NH3 and of ND3 using a new quasi-diabatic representation of the X̃ and Ã states and full six-dimensional quantum dynamics. J. Chem. Phys. 2012, 136, 234301. 49. Ma, J.; Zhu, X.; Guo, H.; Yarkony, D. R., First principles determination of the NH2/ND2(A/X) branching ratios for photodissociation of NH3/ND3 via full-dimensional quantum dynamics based on a new quasi-diabatic representation of coupled ab initio potential energy surface. J. Chem. Phys. 2012, 137, 22A541. 50. Ma, J.; Xie, C.; Zhu, X.; Yarkony, D. R.; Xie, D.; Guo, H., Full-dimensional quantum dynamics of vibrationally mediated photodissociation of NH3 and ND3 on coupled ab initio potential energy surfaces: Absorption spectra and NH2(Ã2A1)/NH2(X̃ 2B1) branching ratios. J. Phys. Chem. A 2014, 118, 11926-11934. 51. Xie, C.; Ma, J.; Zhu, X.; Zhang, D. H.; Yarkony, D. R.; Xie, D.; Guo, H., Full-dimensional quantum state-to-state non-adiabatic dynamics for photodissociation of ammonia in its A-band. J. Phys. Chem. Lett. 2014, 5, 1055-1060. 52. Xie, C.; Zhu, X.; Ma, J.; Yarkony, D. R.; Xie, D.; Guo, H., Communication: On the competition between adiabatic and nonadiabatic dynamics in vibrationally mediated ammonia photodissociation in its A band. J. Chem. Phys. 2015, 142, 091101. 53. Tennyson, J.; Sutcliffe, B. T., The ab initio calculation of the vibrational‐rotational spectrum of triatomic systems in the close‐coupling approach, with KCN and H2Ne as examples. J. Chem. Phys. 1982, 77, 4061-4072. 54. Gatti, F., Vector parametrization of the N-atom problem in quantum mechanics. III. Separation into two subsystems: Application to NH3. J. Chem. Phys. 1999, 111, 7225-7235. 27



55. Skouteris, D.; Barone, V., A new Gaussian MCTDH program: Implementation and validation on the levels of the water and glycine molecules. J. Chem. Phys. 2014, 140, 244104.

28


Page 28 of 39

Page 29 of 39 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


Figure captions: Figure 1. Evolution of the Gaussian widths in imaginary time for different initial guesses for the NOCl system. Figure 2. Convergence of the energy expectation value with respect to imaginary time for expansions with different numbers of GWPs for the NOCl system. Figure 3. NOCl absorption spectra calculated by the G-MCTDH/vMCG method with and without relaxation. The MCTDH result is plotted as the benchmark. Figure 4. Absorption spectrum of NH3 calculated by the modified G-MCTDH/vMCG method using different numbers of relaxed GWPs. The even and odd parities in the exact result are represented by black and red bold lines. The highest of the peak for each parity is scaled to 1.0.

29



Figure 1

30


Page 30 of 39

Page 31 of 39 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


Figure 2

31



Figure 3

32


Page 32 of 39

Page 33 of 39 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


Figure 4

33



TOC Graphic

34


Page 34 of 39

Page 35 of 39 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 1 327x93mm (150 x 150 DPI)



Fig 2 161x123mm (150 x 150 DPI)


Page 36 of 39

Page 37 of 39 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 292x211mm (150 x 150 DPI)



Fig 4 187x123mm (150 x 150 DPI)


Page 38 of 39

Page 39 of 39 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


TOC 35x14mm (300 x 300 DPI)


Modified Gaussian Wave Packet Method for Calculating Initial State

Recommend Documents