Computation of the S1 â S0 Vibronic Absorption Spectrum of

Subscriber access provided by University of South Dakota

Spectroscopy and Excited States 1

0

Computation of the S # S vibronic absorption spectrum of formaldehyde by variational Gaussian wavepacket and semiclassical IVR methods Matteo Bonfanti, Jakob Petersen, Pierre Eisenbrandt, Irene Burghardt, and Eli Pollak J. Chem. Theory Comput., Just Accepted Manuscript • DOI: 10.1021/acs.jctc.8b00355 • Publication Date (Web): 24 Aug 2018 Downloaded from http://pubs.acs.org on August 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 41 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

Computation of the S1 ← S0 vibronic absorption spectrum of formaldehyde by variational Gaussian wavepacket and semiclassical IVR methods Matteo Bonfanti,∗,† Jakob Petersen,‡ Pierre Eisenbrandt,† Irene Burghardt,∗,† and Eli Pollak∗,‡ Institute of Physical and Theoretical Chemistry, Goethe University Frankfurt, Max-von-Laue-Str. 7, D-60438 Frankfurt/Main, Germany, and Chemical and Biological Physics Department, Weizmann Institute of Science, 76000 Rehovot, Israel E-mail: [email protected]; [email protected]; [email protected]

Abstract The vibronic absorption spectrum of the electric dipole forbidden and vibronically allowed S1 (1 A2 ) ← S0 (1 A1 ) transition of formaldehyde is calculated by Gaussian wavepacket and semiclassical methods, along with numerically exact reference calculations, using the potential energy surface of Fu, Shepler and Bowman [J. Am. Chem. Soc. 133, 7957 (2011)]. Specifically, the variational Multiconfigurational Gaussian (vMCG) approach and the Herman–Kluk Semiclassical Initial Value Representation (HK-SCIVR) are compared to assess the accuracy and convergence of these methods, benchmarked against numerically exact time-dependent wavepacket propagation ∗

To whom correspondence should be addressed Institute of Physical and Theoretical Chemistry, Goethe University Frankfurt, Max-von-Laue-Str. 7, D-60438 Frankfurt/Main, Germany ‡ Chemical and Biological Physics Department, Weizmann Institute of Science, 76000 Rehovot, Israel †

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

(TDWP) on the reference potential energy surface. The vMCG calculation is shown to converge quite well with about 100 variationally evolving Gaussian functions and using a local cubic expansion instead of the conventional local harmonic approximation. By contrast, the HK-SCIVR approach with ∼105 trajectories reproduces the vibrationally structured spectral envelope correctly, but yields a strongly broadened spectrum. The comparison of the computed absorption spectrum with experiment shows that the relevant vibronic progressions are reasonably reproduced by all computations, but deviations of the order of 10-100 cm−1 occur, underscoring that both electronic structure calculations and dynamical approaches remain challenging in the calculation of typical small-molecule excited-state spectra by trajectory-based methods.

1

Introduction

Even for small molecular systems, the accurate determination of vibronic absorption spectra in the presence of anharmonic effects is a challenging problem for ab initio methods. When the harmonic approximation turns out to be inadequate, it is necessary to go beyond the description of geometrical effects, i.e., Franck–Condon displacements and Duschinsky rotations, for which the Franck–Condon factors can be computed analytically. 1 One of the most straightforward time–dependent approaches consists in evolving a vertically excited wavepacket on the potential energy surface (PES) of the final electronic state. The spectrum is then computed as the Fourier transform of its autocorrelation function. 2,3 This procedure requires computing the PES within the configuration space regions explored by the wavepacket, and representing it globally as an analytic function of the relevant set of coordinates. This requires, in turn, the calculation of many high–level ab initio data points combined with a multidimensional fitting procedure, using functions that are flexible enough to represent significant variations of the potential surface. Especially when dealing with vibronic spectra this is a tough challenge. Due to the fact that the ground state may be localized in a region which does not overlap the minimum of the excited state and that the 2


Page 2 of 41

Page 3 of 41 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


absorption spectrum covers a range of many thousands of wavenumbers one needs a precise representation of the surface over a rather large range of energies. This is per se a difficult task, although much progress has been made in recent years thanks to the tremendous progress in electronic structure methodologies and sophisticated fitting techniques, such as artificial neural networks 4–6 or more specialized functional forms. 7–9 Once a reliable and complete PES has been obtained, we are faced with the problem of solving the time–dependent Schrödinger equation, which is also a formidable numerical task. Although the Multi–Configuration Time–Dependent Hartree (MCDTH) approach 10,11 and its hierarchical multi–layer extension 12,13 have made it possible to effectively describe tens to hundreds (or more) degrees of freedom, all grid–based methods are affected by the well– known curse of dimensionality, 14 i.e., an exponential scaling problem that sets an inevitable limit to the number of dimensions that can be treated. Furthermore, techniques that are based on a tensor decomposition of the configuration space, like MCTDH, often require a specialized representation of the potential – e.g., in sum-over-products (SOP) form – that poses further difficulties in practical applications. These technical challenges are potentially avoided in so–called direct–dynamics (DD) or on–the–fly methods, where only local information about the PES is needed, which can be computed while propagating the wavepacket without the need of a global representation of the potential. Examples of such dynamical methodologies for use in a DD scheme are the Ab Initio Multiple Spawning (AIMS) 15–17 method, the semiclassical thawed Gaussian approximation (TGA), 18 the variational Multiconfiguration Gaussian 19–22 (vMCG) method and the Herman–Kluk semiclassical initial value representation 23 (HK–SCIVR). In a number of recent studies, these methods have been applied to the simulation of excited state wavepacket dynamics and have shown promising results. 22,24,25 The advantage of these methods is twofold: First, they do not need a global analytic PES and thus can be used without any fitting or interpolation step. Second, they rely on a representation of the wavefunction that is not based on a product grid and thus avoids the exponential scaling with the number

3



of dimensions. In some implementations, a flexible on–the–fly interpolation is coupled to the dynamical method to reduce the number of potential evaluations. 26 Overall, this approach is advantageous over a prior fitting/interpolation step, since no assumptions about the relevant configuration space are required and the construction of the potential is guided by the dynamics. In practice, though, sparse trajectory sets and discontinuities of the PES can pose major problems. As a result, the use of trajectory-based quantum methodologies to obtain spectroscopic rather than mechanistic excited-state information is far from established. In this paper we carry out a comparative study of the vMCG and HK–SCIVR methods for the evaluation of an absorption spectrum that exhibits marked anharmonicities and ranges over an energy scale of ∼10,000 cm−1 and assess their performance against a numerically exact TDWP benchmark that is obtained with a grid based approach for the same PES. We focus on formaldehyde (CH2 O), since it is small enough to be described with the time–dependent wavepacket (TDWP) approach in the (non-rotating) full dimensionality – i.e., in 6D – and exhibits a well characterized symmetry–forbidden and vibronically–allowed S1 ← S0 transition. This n–π ∗ transition from the S0 (1 A1 ) to the S1 (1 A2 ) electronic state borrows intensity through vibronic coupling along the non–totally–symmetric normal modes, in particular the out–of–plane bending mode. The PES along this mode is strongly anharmonic, with a symmetry–broken equilibrium configuration corresponding to a pyramidal geometry of the molecule. As a consequence, the PES shows a typical symmetric double– well feature, that obviously cannot be captured at the harmonic level and thus makes this system an ideal test case for a time–dependent methodology. Experimentally, the UV absorption spectrum of formaldehyde was measured by various groups in the range 240–360 nm that corresponds to the S1 ← S0 transition. 27–29 The peaks of the main vibrational progressions have been assigned and individually analyzed in a number of high–resolution studies, as summarized in a detailed review by Clouthier and Ramsay. 30 Several attempts have been made to reproduce this spectrum with computational methods. Notably, we mention the recent works by Tatchen and Pollak 25 and by Lin et al. 31 In the

4


Page 4 of 41

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


former, the HK-SCIVR method is applied with an approximated unit prefactor, using an on–the–fly evaluation of the potential and forces based on time–dependent density functional theory (TDDFT). A reasonable agreement with experiment is obtained, even though clearly improvable. The present work includes the full Herman-Kluk prefactor, anticipating an improved result. In the work by Lin et al., 31 the spectrum is simulated with a time–independent method by computing the transition rate constant with Fermi’s Golden Rule. The vibronic coupling is treated in a perturbative manner, beyond the harmonic approximation, with input information on the electronic adiabatic states and dipole derivatives obtained from complete active space self-consistent field (CASSCF) calculations. This methodology, specifically tailored to treat the formaldehyde vibronic spectra, yields results that are in good agreement with experiment in the low-frequency range of the absorption spectrum. To make an explicit comparison between the different approaches, we base our calculations on a recent global PES for the S1 state of formaldehyde by Fu et al. 32 which employs a basis of permutationally invariant polynomials developed by the Bowman group. 8 This global S1 potential has been developed to describe the reactive dynamics of formaldehyde on coupled singlet and triplet PESs 32 and was not specifically optimized to describe the Franck-Condon geometry. Hence, we cannot expect quantitative agreement with experiment, as will be further discussed below. Rather than achieving an accurate comparison with the experimental spectrum, our goal is to assess the convergence of the results obtained with the vMCG and HK-SCIVR methods with respect to the TDWP benchmark, using one and the same reference PES. The paper is organized as follows: In Section 2 we discuss the propagation methods and potential surface employed in this comparative study. In Section 3 the vibronic absorption spectra simulated using vMCG and HK–SCIVR are compared with the exact quantum mechanical and experimental spectra. Finally, we provide concluding remarks in Section 4.

5



Page 6 of 41

Several Appendices address aspects of the vMCG and HK-SCIVR calculations.

2

Theory

As explained above, the main purpose of this work is to examine the performance of the vMCG and HK–SCIVR methods in the computation of a prototypical anharmonic spectrum, i.e., the S1 ← S0 vibronic absorption spectrum of formaldehyde, and explore whether these methods can achieve convergence to the numerically exact results. In Sec. 2.1, we start by discussing the general framework for computing the formaldehyde spectrum, and then turn to the main features of the HK–SCIVR (Sec. 2.2) and vMCG (Sec. 2.3) propagation methods. Among the two key challenges of trajectory-based on–the–fly propagation — i.e., the electronic structure accuracy and the dynamics accuracy — we focus on the dynamics side in the present study. To facilitate the propagation procedure, we employ the abovementioned S1 PES by Fu et al., 32 which is based upon extensive multi-reference configuration interaction (MRCI) ab initio data. This PES, which is briefly described in Sec. 2.4, is employed either with pointwise PES evaluations, or else as a global potential for numerically exact benchmark calculations performed in normal-mode coordinates. This permits us to assess the vMCG and HK–SCIVR results against a numerically exact spectrum computed for the same PES using standard quantum dynamical techniques.

2.1

Vibronic Absorption Spectrum

In this study, we focus on the vibronic absorption spectrum for excitation from the electronic ground state in the absence of rovibronic coupling. Within this approximation, the extinction coefficient for excitation with monochromatic light is described as follows by Fermi’s Golden Rule, 3 !

I(ω) = A ω

X i

∆E + Ee,i − Eg,0 −ω , |hie |ˆ µeg | 0g i| δ h ¯ 2

6


(1)

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


where A = (2π)2 NA /(3 log(10)¯ hc(4π0 )), with NA the Avogadro constant. Further, ∆E is the adiabatic electronic energy difference between the excited state (e) and the ground state (g), µ êg is the magnitude of the transition dipole moment between the states g and e which is an operator in the nuclear coordinate space, and |ig i and |ie i are the i–th vibrational eigenˆ g ) and excited-state (H ˆ e ) Hamiltonians with corresponding functions of the ground-state (H eigenvalues Eg,i and Ee,i , respectively. In Eq. (1), we assumed that the initial state is the vibrational ground state |0g i of the electronic ground state. Using the Fourier representation of the Dirac delta function, the absorption spectrum can readily be re-expressed as a Fourier transform 2 X i i Aω Z ∞ dt exp − (∆E − Eg,0 − h ¯ ω) t h0g |ˆ µ|ie i exp − Ee,i t hie |ˆ µ|0g i 2π −∞ h ¯ h ¯ i i Aω Z ∞ dt exp − (∆E − Eg,0 − h = ¯ ω) t χ(t) (2) 2π −∞ h ¯

I(ω) =

where the correlation function χ is given by

ˆ e (t)ˆ χ(t) = h0g |ˆ µeg K µeg |0g i,

(3)

with the time-evolution operator for the excited state,

ˆ e (t) = exp(−iH ˆ e t/¯ K h) =

X i

i |ie i exp − Ee,i t hie | h ¯

(4)

whose expansion in the excited-state eigenfunctions was used in Eq. (2). In the specific case of formaldehyde, the S1 ← S0 transition is dipole forbidden in the C2v symmetry of the ground state, and therefore, the intensity of the vibronic absorption spectrum is primarily due to Herzberg–Teller coupling. 25,33 Thus, we need to approximate the transition dipole moment by using the first order terms of its Taylor expansion

µ êg ≈ µ ˆHT =

N X

µi qg,i ,

i=1

7


(5)


Page 8 of 41

where qg,i is the i–th mass-weighted normal mode coordinate of the ground electronic state and µi = ∂ µ ˆ/∂qg,i is the dipole derivative with respect to qg,i . With this expression for the transition dipole moment, we refer to χ(t) of Eq. (3) as the Herzberg-Teller correlation function. ˆ g can often be accurately represented by a separable The ground vibrational state of H harmonic approximation, i.e.,

E

E

E

, |0g i = 01g ⊗ 02g ⊗ · · · ⊗ 0N g

(6)

E where 0ig is the harmonic approximation of the ground vibrational state corresponding

to the i–th normal mode of the ground electronic state. Thus, |0g i is determined from a normal mode analysis of the optimized molecular geometry in the ground electronic state. In this context we note that the accuracy of such a normal mode approximation is reasonable (see Table 1 below). The error in the harmonic approximation however, does not have a major effect on the absorption spectrum to be presented below for two reasons. One is that although the energy of the ground state used is not accurate, it is in any case not important, since the ab-initio computation does not provide the correct energy gap between the ground and excited states. Therefore, we will present the absorption spectra relative to the lowest spectroscopically allowed transition frequency. Secondly, the accurate shape of the ground state wavefunction will influence intensities of the excited state progressions, but not their energetic location. As will be seen later on, the harmonic ground state is sufficient for a good representation of the overall intensity of the absorption spectrum. Assuming that the dipole derivatives µi have been obtained numerically by finite differencing, the remaining challenge is to determine the time evolution of µ êg |0g i on the PES ˆ e (t)ˆ associated with the excited state, i.e., the K µeg |0g i portion of the correlation function of Eq. (3). In the following, we summarize the two methods in the focus of the present work to evaluate this time evolution.

8


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


2.2

Herman–Kluk Semiclassical Initial Value Representation

Semiclassical theories 34 are based on classical trajectories and are therefore most straightforwardly applicable in an on–the–fly fashion. Semiclassical methods can effectively describe non–classical phenomena, such as quantum interference, and thus are in principle capable of providing a good description of molecular spectra. First examples are the simulation of the absorption spectrum of formaldehyde 25 and the simulation of the vibrational spectrum for the ground state of CO2 by Ceotto and coworkers. 35 We focus here on the well-known Herman–Kluk Semiclassical Initial Value Representation (HK-SCIVR) frozen Gaussian propagator, 23,34 which is a local propagator that is able to simulate vibrational spectra, especially if employed in conjunction with convergence acceleration techniques such as cellularization or filtering. 36–38 In the HK-SCIVR formulation, the time evolution operator is given by

ˆ K(t) =

Z

i dpdq S (pt , qt ) |gΓ (pt , qt )ihgΓ (p, q)|, f Rt (p, q) exp h ¯ (2π¯h)

(7)

where p and q are mass weighted f –dimensional momentum and position variables, respectively, and St is the action integral

St (p, q) =

Z 0

t

pT pτ dτ τ − V (qτ ) , 2 "

#

(8)

along the classical trajectory moving from the phase space point (p, q) to (pt , qt ) in time t, with V the PES in question. The HK prefactor appearing in Eq. (7) is defined as s

R(pt , qt ) =

1 i 1 1 1 1 1 1 1 1 Γ 2 Mqq Γ− 2 + Γ− 2 Mpp Γ 2 − i¯hΓ 2 Mqp Γ 2 + Γ− 2 Mpq Γ− 2 det 2 h ¯

,

(9)

where the so-called monodromy matrices are given by

Mab =

∂a (pt , qt ) , ∂b 9


(10)


Page 10 of 41

where a and b may either be p or q, and Γ is a time–independent f × f width matrix of the coherent state |gΓ (pt , qt )i, which in the coordinate representation is given by

hx|gΓ (pt , qt )i =

det Γ πf

!1/4

i 1 exp − (x − qt )T Γ (x − qt ) + pTt (x − qt ) . 2 h ¯

(11)

When inserting the integral expression of Eq. (7) into Eq. (3), the absolute value squared overlap |hgΓ (p, q)|ˆ µeg |0g i|2 serves as the Monte–Carlo sampling function for generating random initial conditions of the classical trajectories. The trajectories and associated monodromy matrices are integrated using a generalized Leap-Frog algorithm. 39 During the time– evolution, the branch of the square root function in Eq. (9) is chosen in such a way that continuity of the HK prefactor is maintained. 40 The width matrix Γ is chosen to be diagonal with elements corresponding to the S0 normal mode frequencies, i.e., Γii = ωi /¯h.

2.3

Variational Multiconfigurational Gaussian Method

Quantum dynamical techniques based on moving Gaussian wavepacket (GWP) basis sets have attracted much attention over recent years as a promising technique to describe systems with a large number of degrees of freedom in a fully quantum mechanical setting. These methodologies have the advantage of relying on basis functions that are intrinsically local and thus naturally lead to an on–the–fly implementation, as exemplified by the Ab Initio Multiple Spawning (AIMS) 15–17 method and the Direct Dynamics variational Multiconfigurational Gaussian (DD-vMCG) method. 19–22 Furthermore, GWP functions are naturally suited to represent the dynamics of high–dimensional harmonic baths, and this fact was exploited in several methods based on a hybrid Gaussian–grid representations of the wavefunction, as in the Local Coherent State Approximation (LCSA) 41 and the Gaussian-based Multiconfiguration Time-Dependent Hartree (G-MCTDH) method. 42–45 Indeed, the vMCG method was developed as a variant of the more general G-MCTDH scheme. In general, methods based on moving GWPs can be divided into two main categories.

10


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


In the G-MCTDH, vMCG and LCSA approaches, the evolution of all the parameters of the wavefunction ansatz – including the GWP parameters – is dictated by the Dirac– Frenkel/McLachlan Variational Principle 46,47 and is therefore fully quantum. In AIMS, in the Coupled Coherent States (CCS) method 48,49 and in the related Multiconfigurational Ehrenfest (MCE) approach, 50 the evolution of the GWP centers is classical and follows Hamilton’s equations, whereas the expansion coefficients are determined variationally. Although the wavefunction propagation for the latter class of methods is in practice numerically cheaper and more stable, their convergence to the exact solution is complicated by issues regarding the non–conservation of energy 51 and the choice of the initial conditions. 52 For an accurate comparison as necessary in the present study, we therefore choose to work with a fully variational method that has already been applied in an on–the–fly context, i.e., the vMCG approach. Within the vMCG method, 19–22 the f –dimensional wavefunction is expanded as a linear combination of N GWPs, according to the ansatz

Ψ (x, t) =

N X

An (t) h x | gn (t) i

(12)

n=1

where An (t) is a vector of time-evolving linear coefficients and the time-dependent Gaussian wavepackets h x | gn (t) i are defined similarly to the coherent state of Eq. (11), except for a phase factor,

h x | gn (t) i =

(n) (n) h x | gΓ (pt , qt ) i exp

i (n) T (n) (p ) qt h ¯ t

(13)

With a diagonal width matrix Γ as in Eq. (11), the above form corresponds to the standard Frozen Gaussian (FG) wavepackets employed in the vMCG method. The use of FGs rather than Thawed Gaussians (TG) avoids the known numerical instabilities that appear in the variational dynamics of the width parameter. In the following, the set of time-dependent (n)

parameters will be denoted ξt

(n)

= Γ(qt

(n)

(n)

(n)

+ ipt ), labeling a phase-space point (qt , pt ). 11



Page 12 of 41

For completeness, Appendix A compares the present notation with the standard notation employed in the context of the vMCG method. Using the Dirac–Frenkel/McLachlan Variational Principle, we obtain equations of motion for the time dependent An coefficients and parameters ξ (n) , 42 in the form of two coupled first–order differential equations. The linear coefficients An evolve according to ˙ = 1 S−1 (H − ı¯hτ ) A A ı¯ h

(14)

˙ collect the coefficients and their first order derivatives, respecwhere the vectors A and A tively, S is the overlap matrix of the GWPs, i.e., Snm = hgn |gm i, H is the Hamiltonian matrix in the basis of GWPs, Hnm = hgn |H|gm i, and τ is the differential overlap, i.e., τnm = hgn |g˙ m i, which determines the change of S in time and couples the dynamics of the linear coefficients to the evolution of the GWP parameters. The variational dynamics of the ξ (n) parameters is given by the following equation: 1 ξ˙ = C−1 Y ı¯h

(15)

where the matrix C is given by

Cnj,mk = A∗n hgn |xj ( 1 − P0 ) xk |gm i Am

(16)

while the vector Y is defined as

Ynk =

X

A∗n hgn |xk (1 − P0 ) H|gm i Am

(17)

m

where P0 =

P

nm

|gn i(S−1 )nm hgm | is the non-orthogonal projector on the GWP basis.

Particular care needs to be taken regarding the numerical inversion of the matrices C and S, which are frequently ill–conditioned because of two main reasons: first, linear dependencies

12


Page 13 of 41 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


in the GWP basis set – typically occurring when two GWP centers get close in phase space – and second, small populations of the GWP configurations, in particular at the beginning of the simulations when a large set of unoccupied functions are included in the basis set. For a detailed description of these inversion problems and other details of vMCG implementation, we refer to a recent review. 22 Other than these matrix inversions, the working equations of vMCG, Eqs. (14) and (15), involve simple matrix–matrix and matrix–vector multiplications of matrix elements obtained by Gaussian integration. These integrals can be computed with elementary analytic formulas when the potential is expressed as a polynomial function of the coordinates. One of the main advantages of the GWP basis functions is that any arbitrary PES can be cast in this form thanks to their localization in phase–space. More specifically, for the potential matrix ¯ nm (usually the mean elements hgn |V (x)|gm i, we choose a reference point of coordinates x position of the right gm (ket) GWP 42 or the midpoint between the centers of gn (bra) and ¯ nm , so gm (ket) 22 ). We assume that the integrand is non-negligible only in the vicinity of x that the most significant contributions to the integral come from the first few terms of the Taylor expansion of the PES. Using the local expansion of the potential X ∂V 1 X ∂ 2 V ¯ + (xi − x (xi − x ¯i ) + ¯i )(xj − x¯j )+ V (x) ≈ V (x) 2! ij ∂xi ∂xj x¯ i ∂xi x ¯ 1 X ∂ 3V (xi − x ¯i )(xj − x¯j )(xk − x¯k ) + . . . 3! ∂xi ∂xj ∂xk ijk

(18)

¯ x

usually truncated at 2nd or 3rd order, the integral hgn |V (x)|gm i is evaluated analytically as a sum of moments of the Gaussian product function gn∗ gm . The coefficients of this summation ¯ nm , that need to be computed at each step of the propagation are given by the derivatives in x as the GWP centers move in configuration space. In addition to simplifying the computation of Gaussian integrals, the local approximation thus forms the basis for the direct–dynamics use of vMCG. Even though the local harmonic approximation (LHA), i.e., a second order truncation of

13



Eq. (18), is most often used in the context of vMCG, we found it necessary in the present application to adopt a third order expansion of the potential, i.e., a local cubic approximation (LCA). In this way, we were able to minimize the error and the energy oscillations associated with the local approximation error and achieved better convergence to the exact TDWP benchmark. For a detailed comparison between LHA and LCA results, we refer to Sec. B of the Appendix.

2.4

Potential Energy Surfaces

The computation of the absorption spectrum according to Eqs. (2)-(3) requires specifying the vibrational ground-state of the electronic ground-state (S0 ) PES, preceding propagation on the excited-state S1 PES. Therefore, we briefly comment on a harmonic approximation to the S0 state, before addressing the S1 PES by Fu et al. 32 2.4.1

Ground state (S0 ) PES

As discussed above it is sufficient to consider the ground state PES within a harmonic approximation. Normal-mode frequencies and reduced masses are obtained using density functional theory (DFT) with the BP86 functional and def2-TZVP basis set as implemented in the TURBOMOLE 53 package. The normal-mode frequencies compare favorably with the experimental frequencies of Ref. [ 54], as detailed in Table 1. As mentioned above, these frequencies and the associated geometries serve to construct the initial wavepacket that subsequently evolves on the S1 PES. The equilibrium geometry of the S0 PES is reported in Table 2. The results show that the functional that we adopted for the ground state potential reproduces the geometrical parameters in very good agreement with the work by Carter and Handy, who obtained accurate results by an extensive fitting of the rovibrational spectra of H2 CO and D2 CO. 55

14


Page 14 of 41

Page 15 of 41 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


2.4.2

Excited-state (S1 ) PES

The time–dependent calculations presented in the following are based on the S1 PES of formaldehyde computed by Fu, Shepler and Bowman, 32 which was constructed as a permutationally invariant polynomial of Morse–like variables and was fitted to ∼18,500 MRCI(Q)/ccpVTZ ab initio points. According to Ref. [ 32], the vertical excitation energy S1 ← S0 is 3.97 eV (91.5 kcal/mol). In our analysis, the lowest vibronic transition will be taken as reference energy. In the following, we will give a brief assessment of the quality of the PES by comparing some static properties to the available experimental results. A key feature of the S1 potential is the non–planar equilibrium geometry, with a bending angle of approximately 34.01◦ . 56,57 In Table 3 we list the geometric parameters of this symmetry–broken minimum, obtained by numerical optimization of the PES. The results are in excellent agreement with the experimental results of Ref. [ 57] for the bond distances, whereas the bending coordinate ρbend and the HCH angle αHCH do not share the same accuracy. Table 4 reports the harmonic frequencies at the two equivalent minima and at the saddle point located between these minima. At this transition state geometry, the molecule is planar and much closer to the Frank–Condon geometry. The comparison with experimental frequencies shows that some modes are significantly shifted, particularly the high–frequency Table 1: Normal mode frequencies ωDFT and reduced masses mDFT of the S0 state obtained by geometry optimization in C2v symmetry using DFT (BP86 functional and def2-TZVP basis set as implemented in TURBOMOLE 53 ), as compared with experimental frequencies ωEx. taken from Ref. [ 54]. The dipole derivatives µDFT were determined numerically via finite differencing. 25 Mode description ν1 sym C-H stretch ν2 C-O stretch ν3 CH2 scissor ν4 out-of-plane bend ν5 asym C-H stretch ν6 CH2 rock

ωDFT [cm−1 ] 2787 1763 1487 1154 2831 1225

ωEx. [cm−1 ] 2782 1746 1500 1167 2843 1249

15


mDFT [g/mol] µDFT [a.u.] 1.042 7.288 1.112 1.370 4.86 · 10−3 1.122 3.18 · 10−3 1.345 0.90 · 10−3


Table 2: Geometric parameters of the equilibrium geometry of the singlet ground state S0 . Geometric parameter BP86/def2-TZVP Expt. Ref. 55 rCO 1.20834 ˚ A 1.2031 ± 0.0005 ˚ A ˚ rCH 1.11843 A 1.1003 ± 0.0005 ˚ A ◦ αOCH 122.12 121.62 ± 0.05 ◦ Table 3: Geometric parameters of the equilibrium geometry of first excited singlet state S1 . Geometric parameter PES by Fu et al. rCO 1.327 ˚ A rCH 1.096 ˚ A αHCH 113.31◦ ρbend 41.75◦

Expt. Ref. 57 1.3232 ˚ A 1.1028 ˚ A 118.111◦ 34.01◦

stretching modes of the C–H bonds. However, it should be noted that the comparison is more favorable for the CO stretching mode, which is most relevant for the simulation of the absorption spectrum, since the most intense lines reflect the ν2 and ν4 vibrational progressions. 27 Since mode ν4 breaks the planar symmetry and drives the molecule to the bent minima of the PES, the potential along this coordinate is highly anharmonic and has a typical double–well appearance that cannot be described by a harmonic frenquency. Instead, a useful parameter is the height of the energy barrier that connects the two minima, which is computed as 382.54 cm−1 . The barrier is comparatively low, such that all ν4 vibrational states except for the lowest, spectroscopically dark state, lie above the barrier. 30 The computed barrier height is in very good agreement with the value of 350.33 cm−1 that has been extracted by fitting a double well potential to describe the vibrational progression of ν4 . 57 However, despite this agreement we will see in the simulated spectra that there is a significant mismatch in the vibrational quanta of the ν4 normal mode, indicating that the overall shape of the potential along this mode is not accurate enough to reach quantitative agreement with experiment. With regard to the accuracy of the PES, the MRCI level of electronic structure theory is satisfactory as such. However, the PES of Ref. [ 32] was not developed for high–

16


Page 16 of 41

Page 17 of 41 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: Normal modes of the S1 potential energy surface by Fu et al., 32 as compared with experimental data. 30 Frequency / cm−1 (saddle point) CH2 symm. stretching 3177 CO stretching 1162 CH2 scissoring 1363 out–of–plane bending 449i CH2 asymm. stretching 3302 CH2 rocking 919

Mode description ν1 ν2 ν3 ν4 ν5 ν6

Frequency / cm−1 (minimum) 3083 1234 1405 628 3133 1005

Expt. Ref. 30 2846 1183 1293 2968 904

resolution spectroscopic applications but for a mechanistic study of the photodissociation of the molecule, 32 and thus extends much beyond the Franck–Condon region. As a consequence, the PES has a rather sparse sampling of the geometries that are most relevant for the present study; only approximately 2,500 out of the 18,424 total MRCI points lie in the energy range of the absorption spectrum. Moreover, the accuracy of the fitting reported by the authors is 1.1 kcal/mol = 385 cm−1 , which is of course reasonable when considering the extension of the PES, but insufficient for spectroscopic accuracy. Therefore, we anticipate that agreement with the experimental absorption spectrum will be limited, despite the high level of the electronic structure calculations and the advanced fitting technique.

3

Results and Discussion

We now proceed to calculate the Herzberg–Teller absorption spectrum according to Section 2.1, i.e., as the Fourier transform of wavepacket correlation functions. In the following, we first give details about the Herzberg–Teller correlation functions (Sec. 3.1), and then turn to a numerically exact benchmark calculation (Sec. 3.2). Following this, we assess the performance of the HK-SCIVR and vMCG methods (Sec. 3.3).

17



3.1

Page 18 of 41

Details of Herzberg-Teller correlation function

According to Section 2.1, the Herzberg–Teller absorption spectrum is given as the Fourier transform of the autocorrelation function of the wavepacket |Ψi =

P6

i=1

µi qî |0g i, that evolves

in time on the S1 excited state PES. Equivalently, we can use the linearity of the scalar product to express the correlation function as a sum of 36 time–dependent terms,

χ(t) =

6 X 6 X

ˆ e (t)ˆ µi µj h0g |ˆ qi K qj |0g i ≡

i=1 j=1

6 X 6 X

χij (t)

(19)

i=1 j=1

each containing a cross–correlation function χij (t) obtained by projecting the state qˆj |0g i evolved on the excited-state PES onto the intial state qî |0g i. In practice, however, this expression can be considerably simplified. First, three of the dipole derivatives µi vanish due to symmetry and only the non totally–symmetric ν4 , ν5 and ν6 modes contribute to the summation in Eq. (19). Second, in our case we have verified that the cross–terms ˆ e (t)µj qˆj |0g i, i 6= j, are small compared to the autocorrelation terms and can thus h0g |µi qî K be neglected. In view of the above, the Herzberg–Teller spectrum is determined by only three contributions coming from the autocorrelation functions of the S0 ground–state multiplied by one of the ν4 , ν5 and ν6 normal modes,

X

χ(t) ≈

ˆ e (t)ˆ µ2i h0g |ˆ qi K qi |0g i ≡

i∈{4,5,6}

X

χii (t).

(20)

i∈{4,5,6}

Orientational averaging preserves this summation, while dividing the above expression by a factor of 3, i.e., χav = 1/3

P

i∈{4,5,6}

χii . This is not generally valid, though, for Herzberg-

Teller spectra, and this simple relation holds due to the high symmetry of the system and the neglect of cross-terms. The separation into independent terms in Eq. (20) will be employed both in the TDWP benchmark calculations and in the vMCG calculations. We further approximate |0g i by a simple product of one-dimensional Gaussian functions with width parameters given by

18


Page 19 of 41 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 S0 normal mode frequencies, where we refer to Table 1 for the ground-state harmonic frequencies and dipole derivatives µ3 , µ4 and µ5 .

3.2

Benchmark calculations

To compare the approximate techniques under study – HK-SCIVR and vMCG – with a quantum dynamical benchmark, we first obtained the spectrum with a standard TDWP technique. To this end, we computed the three autocorrelation functions of Eq. (20) by numerically exact full–grid propagation, using the Heidelberg MCTDH package. 10,11,58 In these calculations, the wavefunction was expanded with a discrete variable representation (DVR), obtained combining a one–dimensional Gauss–Hermite DVR for all the normal modes except for ν4 , for which we chose a uniform grid that better suits the double well feature of the potential. The PES, expressed as a function of the S0 normal modes coordinates, was directly evaluated at the DVR grid points. Results are shown in Fig. 1 along with the experimental spectrum of Meller and Moortgat 28,59 and the transitions collected by Clouthier and Ramsay. 30 An exponential damping of the autocorrelation functions was applied, multiplying the data in the time–domain by the Gaussian function exp[−t2 /(2∆t)2 ] with ∆t = 125 fs (a time constant that corresponds to a frequency broadening of 42 cm−1 ). Note, though, that the rovibronic structure of the experimental spectrum leads to a more pronounced broadening. The simulated and the experimental spectra are in reasonably good agreement. In particular, they present a similar line structure, especially at low energy, composed of a series of doublets that result from the combination of the ν2 and ν4 progressions. However, at a quantitative level, the agreement is not as satisfactory: closer inspection of the first few lines immediately shows that there is a significant shift in the ν4 progression, whereas the agreement is better for the vibrational excitation of ν2 . This observation suggests that the fitted potential still needs substantial improvement to reach a quantitative comparison between experiment and simulations, in particular with regard to the description of the double–well 19



1

TDWP Experiment

a) Intensity

0.5

0

−0.5 −1

0

1

1

Intensity

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 41

2

3

4

TDWP Experiment

410 441

6

7 210 430

210 410

b) 0.5

5

8

9

430 210 420 610

420 610

420

10 220 430

220 410

510

0 0

1

2 3

3

−1

ω − ω00 [10 cm ] Figure 1: Numerically exact TDWP theoretical absorption spectrum (solid black line) based on the PES of Fu, Shepler, and Bowman 32 as compared with the experimental spectrum (solid red line) adapted from Meller and Moortgat. 28,59 Note that the experimental spectrum illustrates rovibronic fine structure which is not treated in our study. For better comparison, the origins of both spectra are shifted to zero. (a) Full spectral width of the absorption spectrum, (b) low-frequency part of the spectrum where the vibronic transitions reported by Clouthier and Ramsay 30 are marked by vertical blue lines. For a given mode νn (see the assignment of Table 1), experimental vibrational transition energies l ← k are labeled as nlk , and likewise for combination transitions. As can be seen from the low-frequency part of the spectrum, the dominant progression involves the out-of-plane bending mode (ν4 ) together with excitation of the CO stretching mode (ν2 ). Minor participation is observed for the CH2 asymmetric stretching mode (ν5 ) and the CH2 rocking mode (ν6 ).

20


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


feature of the excited state. It should be noted that our simulations do not include the rotational motion of the molecule and cannot reproduce the full rovibrational structure of the absorption spectrum, which is illustrated in the experimental spectrum by Meller and Moortgat 28,59 shown in Fig. 1. However, the comparison between the vibrational lines extracted from the experiment by Clouthier and Ramsay 30 and the full spectrum suggests that centrifugal couplings have a minor effect and cannot account for the disagreement between theory and experiment. Apart from these observations, we stress that the present investigation is concerned with the comparison between different approximate propagation methods, such that the limited agreement with experiment is not an impediment for our study as such.

3.3

HK-SCIVR and vMCG results

We now compare the results of the HK-SCIVR and vMCG computations with the quantum– dynamical benchmark. For both methods, in-house codes were employed. Potential evaluations along the classical paths (in HK-SCIVR) or along the quantum variational paths (in vMCG) are carried out in Cartesian coordinates while the propagation is subsequently carried out in ground-state normal mode coordinates (as is also often done in a genuine on-the-fly propagation, see Ref. [ 22]). In addition to removing rotations and translation, this set of coordinates allows us to construct the initial wavefunction as a simple product of GWPs, as in Eq. (6), and permits a direct comparison with the TDWP benchmark results. As described above, vMCG and HK–SCIVR also require the evaluation of first and second potential derivatives. For this purpose, numerical differentiation by means of finite differences was applied. The derivatives computed in this way were tested numerically and appeared to be stable and robust with respect to the discretization and the order of the finite–difference formula, as could be expected from the smoothness of the analytic potential that is a polynomial of Morse type functions. The HK-SCIVR spectrum was obtained by using up to ∼3·105 classical trajectories to 21



evaluate the phase space integral, where the initial conditions were sampled from the initial wavepacket |Ψi =

P6

i=4

µi qî |0g i, and each trajectory was integrated for 1 ps. On the other

hand, the vMCG spectrum was computed similarly to the TDWP setup, combining the three autocorrelation contributions of Eq. (20). Each of these functions resulted from a 1 ps time propagation of the initial state qî |0g i, with i = 4, 5 and 6.

Intensity

1

TDWP vMCG HK

a)

0.5

0 0 1

Intensity

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 41

1

2

3

4

5

6

7

8

9

10

b)

0.5

0 0

1

2 3

3 −1

4

5

ω − ω00 [10 cm ] Figure 2: (a) Exact quantum mechanical absorption spectrum (TDWP) based on the PES of Fu, Shepler, and Bowman 32 (black) as compared with the approximate spectra computed using the vMCG (red) and HK-SCIVR (blue) methods. The origins of the three spectra are shifted to zero. (b) Expanded representation of the low-energy spectral region. In the vMCG simulations, we adopted GWP functions with the same width as the initial condition, with centers initially arranged on a coordinate product grid (see Sec. B of the Appendix for further details). A larger number of GWPs were used along those coordinates that are more strongly coupled in the dynamics. We employed the maximal number of functions that allowed a stable propagation, without incurring numerical problems related to linear dependencies of the basis set or the escape of GWPs to non-physical attractive PES 22


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


regions resulting from fitting artifacts. Specifically, we used 108 GWP basis functions for the ν4 and ν6 correlation functions, and 60 GWPs for the ν5 correlation function. Further details on the convergence of the calculations with the basis set size are given in Sec. B of the Appendix. The vMCG correlation functions were damped with the same Gaussian function used in the TDWP calculations, i.e. exp[(−t2 /(2∆t)2 ] with ∆t = 125 fs. In contrast, the HK correlation functions, which exhibit an extremely rapid intrinsic decay, were set to zero beyond 150 fs, to remove residual noise. The results of these calculations are shown in Fig. 2, along with the benchmark of Fig. 1. First of all, it turns out that the HK-SCIVR spectrum is not a very accurate approximation of the TDWP benchmark, even though the overall structure is quite well represented. Especially the resolution of the HK spectrum is rather poor with a significant background, and the main peak positions are in general slightly shifted towards higher frequency values. Furthermore, the intensity of these main peaks do not correspond accurately to the TDWP intensities. The resolution is only slightly improved by increasing the number of trajectories. In fact, the 105 sample is quite well converged, such that the broadening of the HK-SCIVR spectrum is not a result of poor convergence. Also, tunneling effects, which pose a major challenge for HK-SCIVR propagation, should not be a dominant factor since the barrier in the symmetry-breaking mode is very shallow. Rather, we attribute the broadening to the lack of unitarity of the HK-SCIVR approximation and the ensuing loss of normalization, due to the dimensionality of the system. The loss of unitarity acts like a damping function with a decay time of ∼ 30 fs. Renormalization would likely improve these results. As shown in Sec. C of the Appendix, using a unit prefactor enhances the norm loss. By contrast, the vMCG results are very close to the TDWP benchmark: the main peaks of the absorption spectrum are well resolved and have the expected intensity. The method, however, does not perform as well for the least intense transitions, because an accurate description of the spectral fine structure is prevented by the broadening of the peaks and

23



the noise of the background. These effects give rise to the presence of small negative peaks, especially in the low frequency region. These errors are due to the basis set truncation and also to the fact that the local harmonic (or cubic) approximation may introduce an unphysical time dependence of the potential, depending on the number and positions of the GWPs. As shown in Sec. B of the Appendix, going from the LHA to the LCA substantially reduces this error. 0.6

a)

0.5

|χ44 |

0.4 0.3 0.2 0.1 0 0.6 TDWP vMCG

b)

0.5

|χ55 |

0.4 0.3 0.2 0.1 0 0.6

c)

0.5 0.4

|χ66 |

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 24 of 41

0.3 0.2 0.1 0 0

50

100

150

200

250

300

350

400

t [fs] ˆ e (t)ˆ Figure 3: Autocorrelation functions |h0g |ˆ qi K qi |0g i| based on TDWP and vMCG are compared; here, the initial wave packet qî |0g i has been renormalized to unity.

24


Page 25 of 41 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 reason for this limitation in resolution becomes clearer when we examine the autocorrelation functions in the time domain, as shown in Fig. 3, where the three vMCG functions ˆ e (t)ˆ |h0g |ˆ qi K qi |0g i| are plotted along with the corresponding TDWP benchmark. It is evident that vMCG reproduces the first few recurrences extremely well, but agreement with the benchmark dynamics then deteriorates. After 100-150 fs, the quality of the results decreases until the method completely fails to reproduce the autocorrelation function. This obviously happens at longer times, when the small intensity and the irregular position of the peaks indicate that the wavefunction has spread over a large portion of the vibrational spectrum of the PES. At this point, we expect the wavefunction to have a complex phase–space structure that is difficult to represent with a limited number of FG GWP functions. The improvement of the basis–set flexibility by a simple increase of the number of GWPs turns out challenging, because the additional functions decrease the numerical stability of the propagation at short time (causing linear dependence problems or increasing the number of unoccupied states). Furthermore, a larger basis–set would result in a significant increase of the computational cost, especially in the context of direct dynamics. Given that this additional cost would be unnecessary for the short–time dynamics, where the method appears to be reasonably well converged, an adaptive basis set would be desirable.

4

Concluding remarks

Following up on earlier efforts to compute the formaldehyde S1 absorption spectrum by onthe-fly semiclassical methods, 25 we have carried out a comparative study of the HK-SCIVR and vMCG methods using the S1 potential by Fu, Shepler, and Bowman. 32 The vibronically allowed S1 absorption spectrum is strongly anharmonic due to the key role of the symmetrybreaking out-of-plane bending mode, which is coupled to several other modes. At the same time, the spectral structures of this 6D system are fairly regular in the low-frequency region, such that we would expect favorable convergence properties for time-dependent methods.

25



To make sure that any errors are not related to insufficient Monte Carlo sampling we have used the fitted form of the excited state potential. As a result the propagation of trajectories was fast and the sample sizes large enough to assure convergence. Our study shows that (i) the HK-SCIVR approach correctly describes the spectral envelope and the main progressions, but a considerable broadening persists even for trajectory ensembles exceeding ∼ 105 trajectories, while (ii) the vMCG approach produces good results with about 100 variationally evolving basis functions. To obtain satisfactory results, vMCG was employed in conjunction with a cubic local approximation to the PES, rather than the standard local harmonic approximation. The lack of accuracy of the HK-SCIVR approximation, which is evident despite the fact that HK prefactor is correctly taken into account in the present study in contrast to the earlier implementation of Ref. [ 25], poses on the one hand a serious limitation to the applicability of the method. On the other hand it does imply that in practice, for the vibronic absorption spectrum of formaldehyde the prefactor is not necessarily needed to compute a low-resolution vibronic spectrum and that one can “get away” with a computation with unit prefactors, which converges with far fewer trajectories and is less sensitive to classically chaotic structures. Is this true in general? Certainly not for ground state dynamics where the typical energy range spanned in the computation is not very large and the HK prefactor is essential for accounting for normalization. However, we speculate that for vibronic transitions which cover a very broad energy range, inclusion of the HK prefactor will not be too important. This does need further clarification. We do note that the good performance in terms of peak positions and intensity encourages the search for techniques designed to accelerate the convergence of the HK–SCIVR approximation, such as the recently devised Multiple Coherent States SC–IVR approach. 60–64 As for the vMCG approach, the cubic numerical scaling 19,22,42 ∼ (N f )3 as a function of the number of GWPs (N ) and the dimensionality (f ) precludes moving much further than in the current application. Indeed, the present vMCG calculations were more expensive

26


Page 26 of 41

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


than the numerically exact TDWP reference calculations. Furthermore, numerical issues make convergence beyond the first few recurrences difficult. In order to apply variational GWP methods to larger systems, more advanced methods are necessary that make use of the mean-field structure of the G-MCTDH approach 42–44 and its hierarchical multi-layer variant. 45 Since these methods typically require a sum-over-products (SOP) representation of the PES, this implies in the context of on-the-fly applications that suitable schemes need to be devised that match the electronic structure information to an SOP representation. Furthermore, the use of adaptive basis sets is desirable at all levels of treatment. Since one of the main purposes of the present study was to assess different semiclassical and GWP-based methodologies, we relied on an existing fitted PES that was employed by Fu et al. 32 for a mechanistic study of the photodissociation dynamics of formaldehyde rather than truly performing on-the-fly computations. Using this surface, we were able to obtain numerically exact quantum results with standard TDWP calculations, as a benchmark to assess the quality of the more approximate HK–SCIVR and vMCG methods. The comparison with the experimental spectrum shows that the quality of the PES is not sufficient, though, to reach quantitative agreement in spectroscopic applications. In line with the RMSD of the PES fit by Fu et al., 32 estimated by the authors to be of the order of hundreds of wavenumbers, we can assume that an important factor contributing to this discrepancy has to do with the accuracy of the fitting. The PES by Fu et al. 32 was constructed for the study of the reactive dynamics of formaldehyde over a broad range of energies and configuration space, and was not optimized for the Franck-Condon region. In this sense, the PES employed in the present study is not representative of state-of-theart potential surfaces which are explicitly constructed to achieve spectroscopic accuracy – sometimes with an accuracy below 1 cm−1 for small molecular systems. 6 In summary, both electronic structure calculations and dynamical approaches remain challenging in the calculation of typical small-molecule spectra by trajectory-based methods. However, it is also clear that semiclassical and GWP-based approaches are able to go

27



Page 28 of 41

far beyond mechanistic on-the-fly studies, and will further benefit from parallel current developments in electronic structure theory, fitting procedures, and quantum dynamics using local basis sets.

Acknowledgement We thank the German-Israeli Foundation for Scientific Research and Development for support of this project under grant number GIF I-1337-302.5/2016. We gratefully acknowledge Joel M. Bowman and Bina Fu for providing the formaldehyde S1 potential energy surface and for their kind assistance in the analysis of the results. We also thank one of the referees for helpful comments regarding the treatment of orientational averaging. M.B. acknowledges fellowship support by the Alexander von Humboldt Foundation.

Appendix A

Standard vMCG notation

In standard vMCG notation, 19–22 the GWP functions of Eq. (13) are written as a product of one–dimensional FGs,

h x | gn (t) i =

f Y

(n)

exp (a

)i x2i

+

(n) (ξt )i xi

+

(n) (ηt )i

!

(A.1)

i=1

(n)

(n)

with negative real (a(n) )i = Γi /2 and complex time-dependent (ξt )i and (ηt )i . As ex(n)

(n)

(n)

plained in the main text, the set of complex parameters (ξt )i = −2(a(n) )i (qt )i + i(pt )i (n)

(n)

label a phase-space point (qt , pt ). Using the above FG form, the set of time-dependent parameters is denoted Λ = {ξ, η}.

28


Page 29 of 41 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 variational dynamics of the GWP parameters is given as 19,22,42 ˙ = 1 C−1 Y Λ ı¯h

(A.2)

˙ is the vector of time derivatives of the parameters, and the general form of the where Λ matrix C and vector Y are specified in Refs. [ 19,22,42]. As detailed in Ref. [ 42], it turns (n)

out that the parameters ηi

are undetermined; indeed, these complex phase parameters

are redundant since they are equivalent to the linear coefficients An of the expansion Eq. (14). Consequently, their dynamics has to be fixed a priori by choosing specific phase and normalization conditions for the GWPs. In particular, the real part is fixed such as to guarantee normalization of the GWPs, while the imaginary part is set to zero in the present study. With these conventions, Eq. (A.2) is subsequently solved for the reduced parameter ˜ vector Λ(t) = {ξ(t)} as explained in the main text.

B

Details on vMCG convergence

The good accuracy of the vMCG calculations was made possible by two aspects: first, the use of a rather large GWP basis and, second, by the use of a local cubic approximation (LCA) for the potential evaluation. We give here some more information on the convergence of the spectra, specifically with respect to these two crucial technical aspects. We fix our attention ˆ e (t)ˆ on the q4 component of the autocorrelation function, i.e. χ44 = µ24 h0g |ˆ q4 K q4 |0g i, that gives the largest contribution to the absorption spectrum. However, similar results were found for the other two components as well. In Fig. 4, the results of vMCG calculations with a varying number of GWP configurations are plotted along with the TDWP benchmark. In all simulations, we arranged the initial centers of the GWPs along a coordinate product grid. The total number of functions is then given by the number of points chosen along each normal mode. Specifically, we indicate by N111411 a calculation in which the grid of the initial centers is constituted by 4 points along 29



TDWP vMCG N111411 vMCG N151411 vMCG N333411

0.6

|χ44 |

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 30 of 41

0.4

0.2

0 0

25

50

75

100

125

150

175

200

225

250

t [fs] ˆ e (t)ˆ Figure 4: Autocorrelation functions χ44 = µ24 h0g |ˆ q4 K q4 |0g i computed from different vMCG calculations with a variable number of GWP basis functions (red, blue and brown lines) are compared to the TDWP result (black line). For a detailed description of the different vMCG setups, we refer to the main text. ν4 and 1 point in all the other coordinates. Likewise, in the N151411 calculation the grid is the product of 4 points along ν4 , 5 points along ν2 and 1 points along the other coordinates, for a total number of 20 basis functions. In the last calculation, 3 points are used for each of the total–symmetric coordinates, with a total number of 108 GWPs. The comparison between these three calculations clearly shows that by increasing the number of GWPs the autocorrelation function gets more structured and converges to the benchmark results. Specifically, while the first recurrence is qualitatively well described even by the smaller basis set (except for a small mismatch in the peak intensity), a larger number of basis functions is required to describe the dynamics at longer time. After the first 4 or 5 recurrences, even the calculation with the largest 108 GWP basis starts to differ significantly from the benchmark autocorrelation function. In Fig. 5 we plot results obtained with the largest basis set N333411 using both LHA and LCA to compute the Hamiltonian matrix elements in an on-the-fly fashion. By comparison with the TDWP benchmark we see that the LHA results are qualitatively correct up to 150 fs. At longer time, the error due to the second–order truncation of the potential builds up and results in an evident degradation of the quality of the autocorrelation function. This long time behavior of the LHA results translates, by Fourier transformation, to a much increased 30


Page 31 of 41

0.6 TDWP vMCG (LCA) vMCG (LHA)

0.4

|χ44 |

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


0.2

0 0

25

50

75

100

125

150

175

200

225

250

275

300

t [fs] ˆ e (t)ˆ q4 K q4 |0g i computed from vMCG calFigure 5: Autocorrelation functions χ44 = µ24 h0g |ˆ culations in which the potential matrix elements are estimated via the LHA (blue line) and the LCA (red line). The graph also reports the TDWP result (black line) for comparison. noise level of the spectrum. In this regard, the LCA results are significantly improved, even though not fully converged.

C

HK-SCIVR results with/without prefactor

As described in the main text, the HK-SCIVR results are characterized by a strong broadening of the spectral peaks. Here, we present a comparison of the full HK-SCIVR results as compared with results obtained with a unit prefactor. We use a slightly reduced trajectory ensemble with 1.5·105 trajectories. Panel (a) of Fig. 6 shows the absolute value of the total autocorrelation function of Eq. (19) normalized to 1 at t = 0, i.e., |χ(t)/χ(0)|. Results are compared for the full HKSCIVR propagator (see Eq. (7)) and the approximate propagator with a unit prefactor, i.e., Rt (p, q) = 1, benchmarked against the TDWP reference. The semiclassical results – both with and without prefactor – correctly reproduce the initial decay of the autocorrelation function, and correspondingly have a good description of the global envelope of the absorption spectrum. The positions of the initial recurrences of the autocorrelation function match the fully quantum TDWP results.

31



TDWP HK UP

a)

|χ(t)/χ(0)|

0.5

0 0

50

100

150

200

t [fs] 1

Intensity

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 32 of 41

b)

0.5

0 0

1

2

3

4

5

6 3

7

8

9

10

−1

ω − ω00 [10 cm ] Figure 6: (a) Absolute value of the normalized autocorrelation function |χ(t)/χ(0)| and (b) the corresponding Fourier transformed spectra obtained from TDWP (black), full HK (red) and HK (UP) with the prefactor approximated to unity (blue). Spectra were obtained from a damped correlation functions in the TDWP case (with ∆t = 125 fs), whereas the HK correlation functions were set to zero beyond 150 fs, to remove residual noise.

32


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


However, both semiclassical autocorrelation functions exhibit a rapid decay in time. In the case of the calculation with a unit prefactor, this significantly reduces the first recurrence and almost completely mask the following ones. This fast decay has been observed before in various semiclassical applications and can be attributed to the lack of norm conservation of the propagator with the approximated prefactor. This behaviour is in principle improved by the use of the full HK prefactor of Eq. (9). In our case, some improvement is indeed observed in the full HK-SCIVR calculations, as illustrated in Fig. 6, but is insufficient to reproduce the autocorrelation function on the ∼200 fs time scale. This is reflected in the corresponding absorption spectra (lower panel of Fig. 6): the HK-SCIVR spectrum has an improved resolution as compared with the unit prefactor calculation, but it is far from fully resolving the spectral peaks. As further discussed in the main text, this is likely due to a significant loss of norm due to the dimensionality of the system.

33



References (1) Avila Ferrer, F. J.; Santoro, F. Comparison of vertical and adiabatic harmonic approaches for the calculation of the vibrational structure of electrronic spectra. Phys. Chem. Chem. Phys. 2012, 14, 13549. (2) Heller, E. J. Photofragmentation of symmetric triatomic molecules: Time dependent picture. J. Chem. Phys. 1978, 68, 3891–3896. (3) Tannor, D. J. Introduction to Quantum Mechanics: A Time-Dependent Perspective; University Science Books, 2006. (4) Behler, J. Neural network potential-energy surfaces in chemistry: a tool for large-scale simulations. Phys. Chem. Chem. Phys. 2011, 13, 17930. (5) Koch, W.; Zhang, D. H. Communication: Separable potential energy surfaces from multiplicative artificial neural networks. J. Chem. Phys. 2014, 141, 21101. (6) Manzhos, S.; Dawes, R.; Carrington, T. Neural network-based approaches for building high dimensional and quantum dynamics-friendly potential energy surfaces. Int. J. Quantum Chem. 2015, 115, 1012–1020. (7) Bowman, J. M.; Carrington, T.; Meyer, H.-D. Variational quantum approaches for computing vibrational energies of polyatomic molecules. Mol. Phys. 2008, 106, 2145– 2182. (8) Braams, B. J.; Bowman, J. M. Permutationally invariant potential energy surfaces in high dimensionality. Int. Rev. Phys. Chem. 2009, 28, 577–606. (9) Peláez, D.; Meyer, H.-D. The multigrid POTFIT (MGPF) method: Grid representations of potentials for quantum dynamics of large systems. J. Chem. Phys. 2013, 138, 014108.

34


Page 34 of 41

Page 35 of 41 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


(10) Meyer, H.-D.; Manthe, U.; Cederbaum, L. The multi-configurational time-dependent Hartree approach. Chem. Phys. Lett. 1990, 165, 73–78. (11) Beck, M. H.; Jäckle, A.; Worth, G. A.; Meyer, H.-D. The multiconfiguration timedependent Hartree (MCTDH) method: a highly efficient algorithm for propagating wavepackets. Phys. Rep. 2000, 324, 1–105. (12) Wang, H.; Thoss, M. Multilayer formulation of the multiconfiguration time-dependent Hartree theory. J. Chem. Phys. 2003, 119, 1289–1299. (13) Manthe, U. A multilayer multiconfigurational time-dependent Hartree approach for quantum dynamics on general potential energy surfaces. J. Chem. Phys. 2008, 128, 164116. (14) Bellman, R. E. Dynamic Programming; Princeton University Press: Princeton NJ, 1957. (15) Mart´ınez, 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. (16) Ben-Nun, M.; Mart´ınez, T. J. Electronic absorption and resonance Raman spectroscopy from ab initio quantum molecular dynamics. J. Phys. Chem. A 1999, 103, 10517–10527. (17) Virshup, A. M.; Punwong, C.; Pogorelov, T. V.; Lindquist, B. A.; Ko, C.; Mart´ınez, T. J. Photodynamics in complex environments: ab initio multiple spawning quantum mechanical/molecular mechanical dynamics. J. Phys. Chem. B 2009, 113, 3280–3291. (18) Wehrle, M.; Oberli, S.; Van´ıˇcek, J. On-the-fly ab initio semiclassical dynamics of floppy molecules: absorption and photoelectron spectra of ammonia. J. Phys. Chem. A 2015, 119, 5685–5690. (19) Worth, G. A.; Burghardt, I. Full quantum mechanical molecular dynamics using Gaussian wavepackets. Chem. Phys. Lett. 2003, 368, 502–508. 35



(20) Worth, G. A.; Robb, M. A.; Burghardt, I. A novel algorithm for non-adiabatic direct dynamics using variational Gaussian wavepackets. Faraday Discuss. 2004, 127, 307. (21) 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. (22) Richings, G.; Polyak, I.; Spinlove, K.; Worth, G.; Burghardt, I.; Lasorne, B. Quantum dynamics simulations using Gaussian wavepackets: the vMCG method. Int. Rev. Phys. Chem. 2015, 34, 269–308. (23) Herman, M. F.; Kluk, E. A semiclasical justification for the use of non-spreading wavepackets in dynamics calculations. Chem. Phys. 1984, 91, 27–34. (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. (25) Tatchen, J.; Pollak, E. Semiclassical on-the-fly computation of the S0→S1 absorption spectrum of formaldehyde. J. Chem. Phys. 2009, 130, 041103. (26) Frankcombe, T. J.; Collins, M. A.; Worth, G. A. Converged quantum dynamics with modified Shepard interpolation and Gaussian wave packets. Chem. Phys. Lett. 2010, 489, 242–247. (27) Strickler, S. J.; Barnhart, R. J. Absolute vibronic intensities in the 1 A2 ← 1 A1 absorption spectrum of formaldehyde. J. Phys. Chem. 1982, 86, 448–455. (28) Meller, R.; Moortgat, G. K. Temperature dependence of the absorption cross sections of formaldehyde between 223 and 323 K in the wavelength range 225-375 nm. J. Geophys. Res. Atmos. 2000, 105, 7089–7101.

36


Page 36 of 41

Page 37 of 41 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


(29) Gratien, A.; Nilsson, E.; Doussin, J.-F.; Johnson, M. S.; Nielsen, C. J.; Stenstrøm, Y.; Picquet-Varrault, B. UV and IR absorption cross-sections of HCHO, HCDO, and DCDO. J. Phys. Chem. A 2007, 111, 11506–11513. (30) Clouthier, D. J.; Ramsay, D. A. The spectroscopy of formaldehyde and thioformaldehyde. Annu. Rev. Phys. Chem. 1983, 34, 31–58. (31) Lin, C.-K.; Li, M.-C.; Yamaki, M.; Hayashi, M.; Lin, S. H. A theoretical study on the spectroscopy and the radiative and non-radiative relaxation rate constants of the S0 1 A1 –S1 1 A2 vibronic transitions of formaldehyde. Phys. Chem. Chem. Phys. 2010, 12, 11432. (32) Fu, B.; Shepler, B. C.; Bowman, J. M. Three-state trajectory surface hopping studies of the photodissociation dynamics of formaldehyde on ab initio potential energy surfaces. J. Am. Chem. Soc. 2011, 133, 7957–7968. (33) Small, G. J. Herzberg–Teller vibronic coupling and the Duschinsky effect. J. Chem. Phys. 1971, 54, 3300–3306. (34) Child, M. S. Semiclassical Mechanics with Molecular Applications; Oxford University Press: Oxford, 2014. (35) Ceotto, M.; Atahan, S.; Shim, S.; Tantardini, G. F.; Aspuru-Guzik, A. First-principles semiclassical initial value representation molecular dynamics. Phys. Chem. Chem. Phys. 2009, 11, 3861. (36) Walton, A. R.; Manolopoulos, D. E. Application of the frozen Gaussian approximation to the photodissociation of CO2 . Chem. Phys. Lett. 1995, 244, 448–455. (37) Walton, A. R.; Manolopoulos, D. E. A new semiclassical initial value method for FranckCondon spectra. Mol. Phys. 1996, 87, 961–978.

37



(38) Thoss, M.; Wang, H. Semiclassical description of molecular dynamics based on initialvalue representation methods. Annu. Rev. Phys. Chem. 2004, 55, 299–332. (39) Brewer, M. L.; Hulme, J. S.; Manolopoulos, D. E. Semiclassical dynamics in up to 15 coupled vibrational degrees of freedom. J. Chem. Phys. 1997, 106, 4832–4839. (40) Kay, K. G. Integral expressions for the semiclassical time–dependent propagator. J. Chem. Phys. 1994, 100, 4377–4392. (41) Martinazzo, R.; Nest, M.; Saalfrank, P.; Tantardini, G. F. A local coherent-state approximation to system-bath quantum dynamics. J. Chem. Phys. 2006, 125, 194102. (42) 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. (43) Burghardt, I.; Nest, M.; Worth, G. A. Multiconfigurational system-bath dynamics using Gaussian wave packets: Energy relaxation and decoherence induced by a finitedimensional bath. J. Chem. Phys. 2003, 119, 5364–5378. (44) Burghardt, I.; Giri, K.; Worth, G. A. Multimode quantum dynamics using Gaussian wavepackets: The Gaussian-based multiconfiguration time-dependent Hartree (GMCTDH) method applied to the absorption spectrum of pyrazine. J. Chem. Phys. 2008, 129, 174104. (45) Römer, S.; Ruckenbauer, M.; Burghardt, I. Gaussian-based multiconfiguration timedependent Hartree: A two-layer approach. I. Theory. J. Chem. Phys. 2013, 138, 064106. (46) Lubich, C. From Quantum to Classical Molecular Dynamics: Reduced Models and Numerical Analysis; European Mathematical Society Publishing House: Zuerich, Switzerland, 2008.

38


Page 38 of 41

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


(47) Bonfanti, M.; Burghardt, I. Tangent space formulation of the Multi-Configuration Time-Dependent Hartree equations of motion: The projector-splitting algorithm revisited. Chem. Phys. 2018, https://doi.org/10.1016/j.chemphys.2018.05.029. (48) Shalashilin, D. V.; Child, M. S. Description of tunneling with the help of coupled frozen Gaussians. J. Chem. Phys. 2001, 114, 9296–9304. (49) 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. (50) Shalashilin, D. V. Nonadiabatic dynamics with the help of multiconfigurational Ehrenfest method: Improved theory and fully quantum 24D simulation of pyrazine. J. Chem. Phys. 2010, 132, 244111. (51) Habershon, S. Linear dependence and energy conservation in Gaussian wavepacket basis sets. J. Chem. Phys. 2012, 136, 014109. (52) Shalashilin, D. V.; Child, M. S. Basis set sampling in the method of coupled coherent states: Coherent state swarms, trains, and pancakes. J. Chem. Phys. 2008, 128, 054102. (53) TURBOMOLE V6 2009, A development of University of Karlsruhe and Forschungszentrum Karlsruhe GmbH, 1989-2007, TURBOMOLE GmbH: Karlsruhe, 2007; available from http://www.turbomole.com. (54) Hadad, C. M.; Foresman, J. B.; Wiberg, K. B. Excited states of carbonyl compounds. 1. Formaldehyde and acetaldehyde. J. Phys. Chem. 1993, 97, 4293–4312. (55) Carter, S.; Handy, N. C. The geometry of formaldehyde. J. Mol. Spectrosc. 1996, 179, 65–72. (56) Job, V. A.; Sethuraman, V.; Innes, K. K. The 3500 transition 1 A2 -1 A1 of formaldehydeH2 , D2 , and HD. J. Mol. Spectrosc. 1969, 30, 365–426. 39



(57) Jensen, P.; Bunker, P. R. The geometry and the inversion potential function of formaldehyde in the and electronic states. J. Mol. Spectrosc. 1982, 94, 114–125. (58) Worth, G. A.; Beck, M. H.; Jäckle, A.; Meyer, H.-D. The MCTDH Package, Version 8.5.5 (2016). See http://mctdh.uni-hd.de. (59) Keller-Rudek, H.; Moortgat, G. K.; Sander, R.; Sörensen, R. The MPI-Mainz UV/VIS spectral atlas of gaseous molecules of atmospheric interest. Earth Syst. Sci. Data 2013, 5, 365–373. (60) Ceotto, M.; Atahan, S.; Tantardini, G. F.; Aspuru-Guzik, A. Multiple coherent states for first-principles semiclassical initial value representation molecular dynamics. J. Chem. Phys. 2009, 130, 234113. (61) Ceotto, M.; Tantardini, G. F.; Aspuru-Guzik, A. Fighting the curse of dimensionality in first-principles semiclassical calculations: Non-local reference states for large number of dimensions. J. Chem. Phys. 2011, 135, 214108. (62) Buchholz, M.; Grossmann, F.; Ceotto, M. Mixed semiclassical initial value representation time-averaging propagator for spectroscopic calculations. J. Chem. Phys. 2016, 144, 94102. (63) Buchholz, M.; Grossmann, F.; Ceotto, M. Application of the mixed time-averaging semiclassical initial value representation method to complex molecular spectra. J. Chem. Phys. 2017, 147, 164110. (64) Buchholz, M.; Grossmann, F.; Ceotto, M. Simplified approach to the mixed timeaveraging semiclassical initial value representation for the calculation of dense vibrational spectra. J. Chem. Phys. 2018, 148, 114107.

40


Page 40 of 41

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


Graphical TOC Entry

41


Computation of the S1 â S0 Vibronic Absorption Spectrum of

Recommend Documents

Computation of the S1 â S0 Vibronic Absorption Spectrum of