Multidimensional quantum dynamics and infrared spectroscopy of hydrogen bonds

Multidimensional quantum dynamics and infrared spectroscopy of hydrogen bonds

Physics Reports 430 (2006) 211 – 276 www.elsevier.com/locate/physrep Multidimensional quantum dynamics and infrared spectroscopy of hydrogen bonds K...

2MB Sizes 0 Downloads 36 Views

Physics Reports 430 (2006) 211 – 276 www.elsevier.com/locate/physrep

Multidimensional quantum dynamics and infrared spectroscopy of hydrogen bonds K. Giesea , M. Petkovi´ca, b , H. Naundorfa , O. Kühna,∗ a Institut für Chemie und Biochemie, Freie Universität Berlin, Takustr. 3, D-14195 Berlin, Germany b Faculty of Physical Chemistry, University of Belgrade, Studentski trg 12, 11 000 Belgrade, Serbia

Accepted 21 April 2006 editor: S. Peyerimhoff

Abstract Hydrogen bonds are of outstanding importance for many processes in Chemistry, Biology, and Physics. From the theoretical perspective the small mass of the proton in a hydrogen bond makes it the primary quantum nucleus and the phenomena one expects to surface in a particular clear way are, for instance, zero-point energy effects, quantum tunneling, or coherent wave packet dynamics. While this is well established in the limit of one-dimensional motion, the details of the multidimensional aspects of the dynamics of hydrogen bonds are just becoming accessible to experiments and numerical simulations. In this review we discuss the theoretical treatment of multidimensional quantum dynamics of hydrogen-bonded systems in the context of infrared spectroscopy. Here, the multidimensionality is reflected in the complex shape of linear infrared absorption spectra which is related to combination transitions and resonances, but also to mode-selective tunneling splittings. The dynamics underlying these spectra can be unravelled by means of time-resolved nonlinear infrared spectroscopy. As a fundamental theoretical ingredient we outline the generation of potential energy surfaces for gas and condensed phase nonreactive and reactive systems. For nonreactive anharmonic vibrational dynamics in the vicinity of a minimum geometry, expansions in terms of normal mode coordinates often provide a reasonable description. For reactive dynamics one can resort to reaction surface ideas, that is, a combination of large amplitude motion of the reactive coordinates and orthogonal harmonic motion of the remaining coordinates. For isolated systems, dynamics and spectroscopy follow from the time-dependent Schrödinger equation. Here, the multiconfiguration time-dependent Hartree method is shown to allow for describing the correlated dynamics of many degrees of freedom. Classical trajectory based methods are also discussed as an alternative to quantum dynamics. Their merits and shortcomings are scrutinized in the context of incorporating tunneling effects in the calculation of spectra. For the condensed phase, reduced density operator based approaches such as the quantum master equation are introduced to properly account for the energy and phase relaxation processes due to the interaction of the hydrogen bond with its surroundings. © 2006 Published by Elsevier B.V. PACS: 33.20.Ea; 03.65.Sq; 03.65.Xp; 03.65.Yz; 82.20.−w; 82.30.Rs; 82.53.−k Keywords: Quantum dynamics; Intramolecular energy redistribution; Tunneling; Energy and phase relaxation

∗ Corresponding author.

E-mail address: [email protected] (O. Kühn). 0370-1573/$ - see front matter © 2006 Published by Elsevier B.V. doi:10.1016/j.physrep.2006.04.005

212

K. Giese et al. / Physics Reports 430 (2006) 211 – 276

Contents 1. Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2. Gas phase hydrogen bond Hamiltonian . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.1. Anharmonic coupling . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.2. Model potentials for hydrogen transfer . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.3. Reaction surface Hamiltonian . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.3.1. Reaction coordinates . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.3.2. A full-dimensional all-Cartesian Hamiltonian . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.3.3. Fixed reference case . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.3.4. Flexible reference case . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.3.5. Selection of relevant modes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.3.6. Atomic reaction coordinates . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.3.7. Collective reaction plane coordinates . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.3.8. Reduced normal modes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3. Coherent quantum dynamics . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.1. Mean-field vs. multiconfiguration methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.2. Laser-driven dynamics and intramolecular vibrational energy redistribution . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.2.1. The case of a single minimum potential . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.2.2. The case of an asymmetric double minimum potential . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.3. Mode-specific tunneling splittings . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4. Semiclassical and classical trajectory-based methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.1. Semiclassical approximation for the calculation of spectra . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.2. Classical trajectory approach for tunneling . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.2.1. The Makri–Miller model (MM) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.2.2. The extended Makri–Miller model (EMM) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.2.3. Application of the MM and EMM methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5. Condensed phase Hamiltonian . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6. Dissipative quantum dynamics of hydrogen bonds . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.1. Lineshape theories and nonlinear response functions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.2. Quantum master equation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.3. Vibrational energy cascading in an intramolecular hydrogen bond . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.4. HOD in D2 O and H2 O . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.5. HAT and PT reactions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.6. Multidimensional IR spectroscopy . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7. Sundry topics . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7.1. Geometric isotope effects . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7.2. Excited state PT . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7.3. Laser control of hydrogen transfer . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8. Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Acknowledgment . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Appendix A. Determination of the SMC Hamiltonian parameters . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Appendix B. Calculation of multidimensional eigenstates . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Appendix C. Parities of motion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . References . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

212 215 215 219 221 221 224 226 226 226 227 228 228 229 229 231 231 234 237 243 243 245 245 247 249 253 255 255 257 259 261 263 263 265 265 266 267 268 269 269 269 270 270

1. Introduction Hydrogen bonds (HBs) are of pivotal relevance for the understanding of the equilibrium and nonequilibrium properties of many molecular systems. The scope of HB research is really cross-disciplinary encompassing physics, chemistry, and biology [1–12]. The first published record of the concept of a HB as a new type of weak bond goes back to a paper by Latimer and Rodebush in 1920 [13,14]. But it was only in the 1930s that the interest revived, most notably in the context of proton and hydroxyl ion mobility in water [15]. Their unusually high mobility had been known since the work of Grotthus in the early 19th century, but the first attempt for a microscopic understanding in terms of Born–Oppenheimer potential energy surfaces (PES) for the nuclear motion was provided by Huggins in 1936 only [16]. The breakthrough in recognition was probably the discourse upon HBs in Pauling’s “ The Nature of the Chemical Bond” [17]. It was also in the 1930s that the signatures of the formation of HBs in stationary infrared (IR) absorption

K. Giese et al. / Physics Reports 430 (2006) 211 – 276

213

spectra had been realized (see, e.g., Ref. [18]). Soon thereafter hydrogen-bonding found its central place in biology with its role in protein structure [19] and in particular for understanding the DNA [20]. There is no space to cover the many facets of the development in the years that followed. The state-of-affairs as of 1976 had been summarized in Ref. [1] and an up-to-date compilation is forthcoming [12]. A new perspective on HBs emerged with the development of ultrafast time-resolved spectroscopies that started in the 1980s. Having a time-resolution in the pico- and femtosecond domain it became possible to unravel the dynamical processes underlying broad and sometimes structured stationary absorption spectra in particular in the condensed phase. First, dictated by experimental constrains, it has been proton transfer taking place after electronic excitation which was on the agenda [21,22]. During the last years, however, laser sources with the appropriate time resolution in the range between 500 to 4000 cm−1 became available to nonlinear IR spectroscopy, opening the road to the observation of real-time HB dynamics in the electronic ground state and thus to the differentiation between the various mechanisms contributing to the lineshape [23]. It goes without saying that this development on the experimental side has triggered substantial theoretical efforts. Apart form the general complexity of condensed phase systems, theory is challenged by two specific features of HBs which are already relevant in the gas phase: their quantum nature which often comes along with multidimensionality, that is, the coupled motion of many degrees of freedom (DOF). Interestingly, an early account on the multidimensional nature of HB dynamics comes from the 1936 paper by Huggins [16]. Fig. 1 shows the (empirical) potential energy curve developed in Ref. [16] for the motion of a proton between two oxygen atoms as a function of the distance between the latter. Apparently, the shape of this potential changes from

Fig. 1. Empirical potential energy curves for the motion of a proton between two oxygen atoms in a HB (the proton coordinate is defined with respect to the center of the O–O distance vector). Upon decreasing the O–O distance the shape of the potential changes from a double to a single minimum (reprinted with permission from Ref. [16]. Copyright 1936 American Chemical Society).

214

K. Giese et al. / Physics Reports 430 (2006) 211 – 276

(a)

energy

energy

energy

a double minimum for long distances to a single minimum for short ones. This has, of course, dramatic consequences for the proton wave function and the associated properties (for a respective recent simulation, see also Ref. [24]). Now, let us assume a situation where an O · · · O motion has been initiated, for instance, by some laser excitation [23] or is driven by the fluctuations of the environment [25]. For the proton wave function this gives rise to a time-dependent potential such that, e.g., the respective energy levels are modulated in time. The actual effect on the IR spectrum might range from a broadening to the appearance of some substructure, the details will depend on the time scale of the modulation as compared to the time scale of the O–H vibrational motion (∼ 10–20 fs) as well as on the type and strength of anharmonic coupling. Thus the presence of a HB, A.H · · · B, is reflected in the change of a number of properties as compared with the case of a free A–H bond. This fact lends itself to the classification into weak, medium-strong, and strong HBs [4]. Taking the perspective of IR spectroscopy, one observes a red-shift, a broadening, and the development of an often complex substructure of the band associated with the vibration of the bridging H atom. Although the actual numbers differ in ˚ with the literature (compare, e.g., Refs. [4,26]) weak HBs can be characterized by a HB distance exceeding 2.8–3.2 A A–H stretching band being in the range between 3300 and 3500 cm−1 . Medium strong or moderate bonds have a HB ˚ and an absorption in the range of 2600–3300 cm−1 , while in strong HBs the bond length length of about 2.5–3.2 A ˚ is smaller than 2.5 A and the absorption is found in the range from 2100 cm−1 down to about 1000 cm−1 . Specific examples will be discussed in the following. The strength of the HB comes along with a typical shape of the PES for the motion of the H atom. For weak HBs having a relatively large HB distance one often finds a symmetric (A· · ·A) or asymmetric (A · · · B) double minimum potential as shown for the symmetric case in Fig. 2a. In particular in the symmetric case one has besides the local vibrational motion to account for the tunneling between the different wells. Viewed from the time-domain perspective this gives a contribution to the H transfer probability as indicated, e.g., by the deviation of the temperature dependence of the rate coefficients from the Arrhenius behavior [27]. In the frequency domain, tunneling gives rise to a doubling of vibrational states and shows up, e.g., in corresponding microwave [28], infrared [29], or fluorescence excitation [30] spectra. The tunneling splitting  is related to the tunneling period T via  = h/T . Medium strong HBs have a much lower barrier and therefore a larger tunneling splitting is to be expected. The case of an asymmetric bond is shown in Fig. 2b. Here, the dynamics is mostly that of a rather anharmonic vibration. However, due to the coupling between the A–H and A · · · B motions, this dynamics will be multidimensional as well. For strong HBs the barrier becomes small or negligible and the zero-point energy (ZPE) of the anharmonic vibrational motion may exceed the barrier height as shown in Fig. 2c. In general one distinguishes hydride, H atom (HAT), and proton (PT) transfer depending on whether two, one or no electron charge is transferred along with the positively charged nucleus [31]. The former two cases can be viewed as proton-coupled electron transfer reactions and, indeed, often the distinction is not very clear and one has to be aware that, for instance, the transfer of the proton can be coupled to the simultaneous transfer of the respective negative charge along a different pathway, e.g., in small organic molecules through the conjugated system [32]. This mechanistic feature has, of course, drastic consequences for the actual reaction especially in a polar environment [25]. Leaving aside these electrostatic aspects and returning to the multidimensional picture of nuclear dynamics one can envisage

reaction coordinate

(b)

reaction coordinate

(c)

reaction coordinate

Fig. 2. Potential energy curves typical for weak symmetric (a), moderate asymmetric (b), and strong symmetric (c) HBs. Vibrational energy levels in of the potentials are also sketched.

K. Giese et al. / Physics Reports 430 (2006) 211 – 276

215

two fundamental types of couplings between the HAT or PT reaction coordinate and motions of the surrounding [33]. First, there will be promoting modes which modulate the potential along the reaction coordinate such as to reduce the barrier (cf. Fig. 1). The second type of modes are those whose motion is required to change the molecular structure from the reactant to the product configuration, e.g., accounting for the rearrangement of single and double bonds. These are also called reorganization modes [33]. The effect of these couplings can be observed as mode-specific tunneling splitting, for instance, in electronic excitation spectra of tropolone which has an intramolecular HB [34]. It has been associated with the promotion of H atom tunneling along ammonia chains [35] and in a broader sense it also plays a role, e.g., in hydride transfer in enzyme catalysis [36,37]. The multidimensionality of HB quantum dynamics and its spectroscopic signatures is the central theme of this review. First, we will consider the problem of defining and calculating multidimensional PES in the gas phase in Section 2. This section is divided into the cases of nonreactive (Section 2.1) and reactive (Sections 2.2 and 2.3) PES. Special emphasis is put on the detailed derivation of a Cartesian reaction surface (CRS) Hamiltonian, which provides a straightforward road to a full-dimensional simulation of about 100–200 DOF by ab initio quantum chemical methods. Having defined the Hamiltonian, the solution of the time-dependent multidimensional nuclear Schrödinger equation by means of the multiconfiguration time-dependent Hartree (MCTDH) method is briefly presented in Section 3.1. It is emphasized that by construction the CRS Hamiltonian is ideally suited for the combination with MCTDH wave packet propagations. First, this is illustrated by exemplary calculations of ultrafast IR laser-driven dynamics initiating intramolecular vibrational energy redistribution (IVR) in a single and a double minimum HB system in Section 3.2. Next, the issue of vibrationally mode-specific tunneling is addressed in Section 3.3. Using the case of intramolecular HAT in tropolone as an example, the determination of an appropriate CRS Hamiltonian is exercised in some detail. The favorable scaling of classical molecular dynamics with the number of DOF keeps triggering efforts aimed to incorporate approximate quantum effects into such simulations. In Section 4 we first briefly review semiclassical approaches which are built on the WKB formalism. This sets the stage for the subsequent discussion of different types of tunneling trajectories whose performance is tested for the case of double HAT in carboxylic acid dimers. Condensed phase processes are addressed in Sections 5 and 6. First, the system–bath approach to condensed phase problems is introduced in Section 5. Next, we summarize different line shape models of HB systems as well as the response function formalism, vital for the analysis of ultrafast nonlinear IR spectroscopies (Section 6.1). The response functions can be simulated by means of a quantum master equation (QME) which is discussed in Section 6.2. As an application we give results for the ultrafast relaxation dynamics of an intramolecular HB system which highlights its multidimensional nature by virtue of a solvent-assisted energy cascading after OH-stretch excitation (Section 6.3). Subsequently, we give a brief summary of the vibrational dynamics of HOD in D2 O and H2 O which has been rather extensively studied by means of nonlinear IR spectroscopy in recent years. After a short account on reactive PT and HAT dynamics in the condensed phase has been given in Section 6.5 we focus on multidimensional IR spectroscopy as a novel means to unravel complex spectra in the condensed phase. In particular it is shown how 2D-IR spectroscopy applied to a generic dissipative HAT model can be used to scrutinize details of the system–bath coupling. While this review is focussed on ground state dynamics and IR spectroscopy we have summarized closely related topics in Section 7. First, we address the properties of the vibrational ground state wave function itself and discuss how multidimensionality appears as the interdependence of primary and secondary geometric isotope effects (Section 7.1). Second, we briefly quote some results on electronic excited state dynamics where extremely fast PT is initiated by photoexcitation (Section 7.2). Third, we give a perspective on laser control of PT and HAT in Section 7.3. Finally, we provide a synopsis in Section 8. 2. Gas phase hydrogen bond Hamiltonian 2.1. Anharmonic coupling In the following we will focus on the properties of adiabatic PES for the electronic ground state assuming the validity of the Born–Oppenheimer approximation for the decoupling of electronic and nuclear DOF [38]. In passing we note that there are recent promising attempts to treat electrons and protons on the same footing by separating off the coordinates of the heavy atoms only (for an overview, see Ref. [39]). The vibrational properties of molecules in the vicinity of a stationary point on the Born–Oppenheimer PES are often quite accurately described by means of normal mode coordinates. Denoting the Cartesian coordinate vector for the Nat atoms as x = (x1 , . . . , xNat ) with xi = (xi,x , xi,y , xi,z )

216

K. Giese et al. / Physics Reports 430 (2006) 211 – 276

√ and the mass-weighted Cartesian vector as X with Xi, = xi, mi ( = x, y, z), the normal mode vector Q is defined as a linear combination of Cartesian displacements from some reference geometry Xref by Q = L(X − Xref ),

(1)

where L is the unitary transformation matrix containing the eigenvectors of the matrix of the second derivatives of the potential (Hessian) at the reference geometry. At a minimum of the PES, one has 3Nat − 6 eigenvectors describing vibrational motion and, respectively, three eigenvectors for infinitesimal rotation and translation. The correlation between the HB length and the shape of the potential for the H motion expressed in Fig. 1 strongly suggests that the dynamics of HBs and in particular their IR spectra can only be understood in terms of multidimensional PES for the coupled DOFs. In other words, the harmonic approximation for the motion of the set of normal mode coordinates, {Qi }, will not be adequate for the present purpose. Nevertheless, it provides a starting point for developing a principal understanding in terms of simple models. Early qualitative ideas have been developed in the 1940s [40] and were reviewed by Sheppard [41] as well as Bratos and Hadži [42]. The first quantitative assessment was given by Witkowski and Marechal for the IR spectra of carboxylic acid dimers [43,44]. The standard anharmonic model Hamiltonian [26,45] for discussing the IR spectra of weak HBs consists of an A–H stretching oscillator, Q , coupled to an A · · · B oscillator, Q , via a cubic anharmonicity K : H = T + T + 21 2 Q2 + 21 2 Q2 + K Q2 Q ,

(2)

where T and T denote the respective uncoupled kinetic energies containing the conjugate momenta {Pj }. Notice that this treatment neglects any effect of HAT or PT. Next, one assumes that the dynamics of Q is much faster than that of Q . As a consequence it is reasonable to rewrite Eq. (2) by introducing an effective frequency for the fast coordinate:  2K eff (Q ) =  1 + Q . (3) 2 To simplify matters this expression is considered in a Taylor expansion up to first order in Q only which gives eff (Q ) ≈  +

K Q . 

(4)

Making use of the time scale separation between the two modes one can invoke the so-called second Born–Oppenheimer (adiab) separation. This amounts to introducing the adiabatic quantum states, AH (Q ), of the fast mode Q which depend parametrically on the slow mode coordinate Q . For the respective adiabatic energies one has E(adiab) (Q ) = (AH + 21 )eff (Q ). AH This in turn gives the following shifted oscillator PES for the motion of the slow mode    1 2 2 1 K VAH (Q ) =  Q + AH +  + Q . 2 2 

(5)

(6)

Notice that it is to be expected that the HB contracts upon excitation since the elongation of the A–H bond in the excited state of the anharmonic PES will lead to a strengthening of the HB. This will lead to a shortening of the A · · · B distance (0) which is reflected in the minima of the oscillator PES in Eq. (6) which are located at Q = −(AH + 1/2)K / 2 . The situation is sketched in Fig. 3 where we also introduced the quantum states of the slow mode. From this figure it becomes clear that the problem has been reduced to a model which is well-known from electron-vibrational spectroscopy [38]. Therefore, it is not surprising that the understanding of the related IR spectra heavily borrows from respective theories for electronic transitions. For instance, based on Fig. 3 one would anticipate that the IR spectrum shows a Franck–Condon type progression of lines with the intensity pattern depending on the relative shift of the PES. This holds in particular if it is assumed that the dipole moment function depends on the fast coordinate only, i.e., (Q , Q ) ≈ 0 +

j Q . jQ

(7)

Upon increasing the HB strength the IR spectrum may become broader and more structured, which signals that additional anharmonicities have to be included. Besides the anharmonicity of the A–H stretching vibration, it is the so-called Fermi

K. Giese et al. / Physics Reports 430 (2006) 211 – 276

217

energy

νAH =1

νAH = 0

A .... B Fig. 3. The adiabatic separation between the fast A–H oscillator coordinate and the slow A · · · B coordinate allows for the definition of shifted oscillator PES for each quantum state of the fast coordinate.

resonance coupling between the A–H stretching fundamental and the A–H bending first overtone transition which is important. Denoting the bending coordinate as Q we have the modified Hamiltonian H = T + T + T + 21 2 Q2 + K Q3 + 21 2 Q2 + 21 2 Q2 + K Q2 Q + K Q Q2 .

(8)

Now there are two fast modes and the adiabatic separation scheme needs to be modified. This is most conveniently done by introducing appropriate zero-order states for the uncoupled stretching and bending motions. Denoting the latter by |AH  and |AH , respectively, and expanding the adiabatic wave function in terms of the product basis {|AH |AH } yields a Hamiltonian matrix for the fast mode where the zero-order states are coupled via matrix elements of the type AH |K Q Q2 |AH . Diagonalization of this matrix will yield eigenstates where the stretching fundamental and the bending overtone zero-order states are mixed. This usually leads to oscillator strength borrowing of the bending overtone from the stretching fundamental which gives a typical double peak structure in the spectrum. The region between the two peaks is commonly called Evans window. A second quantization description of this approach has been given in Refs. [46,47] (cf. the review in Ref. [48]). In principle this strategy can be extended to encompass strong HBs as well, e.g., by incorporating additional modes or by including higher-order terms in the PES and dipole moment surface [26]. In any case, one should be aware that the picture given so far may fail due to oversimplification when it comes to a quantitative comparison with the experiment based on ab initio quantum chemical anharmonic force constants. This does not mean that force constants treated as empirical parameters would not have given a reasonable fit to some experimental data. The importance of higher-order terms has been nicely illustrated for the moderate HB in the acetic acid dimer in Ref. [49], where force constants up to 6th order have been necessary to explain the fine structure of the spectrum in the OH-stretching region. The expansions of the PES in terms of normal mode coordinates discussed so far are in fact only certain approximations of a more general PES V (Q1 , . . . , Q3Nat −6 ). Expressed in normal mode coordinates the general form of the

218

K. Giese et al. / Physics Reports 430 (2006) 211 – 276

Fig. 4. IR stick spectrum of malonaldehyde as calculated from a 4D model PES including the OH-stretching as well as the O · · · O mode shown in the upper part of the figure. The PES does not include the possibility of HAT. (For more details, see also Ref. [52].)

Hamiltonian has been given by Watson [50]. It reads H=

1 2

 , =x,y,z

(J −  )I (J − ) −

3Nat −6 h¯ 2  1  I + Pi2 + V (Q1 , . . . , Q3Nat −6 ). 8 =x,y,z 2

(9)

i=1

Here, J are the rotational operators, the 3 × 3 matrix I is the inverse of the moments of inertia tensor, and  = Nat    ij ij Qi Pj is the so-called vibrational angular momentum. It contains the Coriolis coupling constants ij = k=1 × Lk ,i Lk ,j with  being the unit antisymmetric tensor and L has been given in Eq. (1) (note that in the present case the composite index k (k = 1, . . . , Nat , = x, y, z) is used). Being interested in intramolecular HB vibrations of larger molecules one can restrict the discussion to the case of no rotation, J = 0. Furthermore, since the first two terms of Eq. (9) are proportional to the inverse of the moments of inertia they are usually assumed to be negligible. Notice, however that especially for small molecules this will not be the case (see, e.g., Ref. [51]). In principle the PES entering Eq. (9) can be calculated on a numerical grid which circumvents the calculation of higher-order derivatives. An example is given in Fig. 4 where we show the IR absorption spectrum of malonaldehyde calculated on the basis of a full 4D PES containing the OH-stretching, OH-in- and out-of-plane bending as well as a low-frequency O · · · O mode [52]. Most interesting in terms of H-bonding is the region between 2700 and 3100 cm−1 . Here, one finds three almost equally spaced lines which would suggest a Franck–Condon progression in the OHstretching vibration with respect to the O · · · O mode. However, the spacing is not fitting the O · · · O frequency and the

K. Giese et al. / Physics Reports 430 (2006) 211 – 276

219

fact that the two spacings are equal is a coincidence. Analysis of the respective eigenstates reveals the strong mixing in particular with the in-plane OH-bending vibration. This situation resembles Eq. (8), although the PES is substantially more anharmonic [52]. With increasing number of DOFs, the effort to calculate the full-dimensional PES increases exponentially. Nevertheless, since the importance of correlations between modes is likely to decrease with the number of participating modes it is reasonable to expand the PES with respect to such correlations [53]:    V (Q1 , . . . , Q3Nat −6 ) = V (1) (Qi ) + V (2) (Qi , Qj ) + V (3) (Qi , Qj , Qk ) + . . . , (10) i

i
i
where the so-called n-mode correlation potentials, V (n) , have been introduced and the summations run over distinct pairs, triples, etc. In passing we note that any force field PES based on a Taylor expansion as discussed above has this form. In principle, normal mode coordinates can also be used for reactive problems. For this purpose, it is advantageous to use the transition state geometry as the reference for the Taylor expansion. This strategy has been applied to e.g. the calculation of the infrared spectrum of H3 O− 2 in Ref. [54]. It is clear, however, that normal mode coordinates are not likely to give a very compact representation of the PES away from the reference point. Here, proper large amplitude coordinates usually will perform much better as will be discussed below. Finally, we remark on the possibility of extracting normal modes and associated frequencies from harmonically driven classical molecular dynamics which avoids the calculation of second derivatives altogether [55]. 2.2. Model potentials for hydrogen transfer There is a widely used set of two-dimensional model PES which captures different coupling types between large and small amplitude motions (see, e.g., Refs. [56–59]). The dynamics and the associated spectra depend, of course, on the symmetry of the PES. In the following we will consider symmetric potentials appropriate for A–H · · · A situations. This will provide the background for the discussion of the tunneling splittings in Sections 3.3 and 4. The asymmetry required to model the case A–H · · · B can be easily incorporated by adding, e.g., a term linear in the reaction coordinate. Starting point is the most simple polynomial expression of a symmetric double-well potential given by the square-quartic Hamiltonian, V (1) (X) = − 21 aX 2 + 41 cX 4 ,

(11)

where only √ the large amplitude coordinate X is considered; a and c are constants that determine the two minima at Xmin = ± a/c and the barrier height E = a 2 /4c. The saddle point is at X = 0. A small amplitude vibration along the mass-weighted normal coordinate Q is modeled by the harmonic oscillator PES, V (1) (Q) = 21 2 Q2 ,

(12)

where is the frequency of the oscillator. If the two coordinates are coupled, the dynamics is governed by the Hamiltonian H=

PQ2 PX2 + + V (X, Q) 2 2

(13)

with the PES having the form, Eq. (10), i.e. V (X, Q) = V (1) (X) + V (1) (Q) + V (2) (X, Q).

(14)

The large amplitude coordinate is anti-symmetric with respect to the molecular symmetry transformation T (permutation and subsequent rotation) while the small amplitude coordinate Q can either be symmetric or anti-symmetric with respect to T. Thus, if there is a minimum at (X0 , Q0 ) there must be an equivalent minimum either at (−X0 , Q0 ) or at (−X0 , −Q0 ). The case V (2) (X, Q) =  X 2 Q corresponds to the symmetrical mode coupling (SMC) with  denoting the coupling constant and Q is assumed to be symmetric with respect to the molecular symmetry transformation T.

220

K. Giese et al. / Physics Reports 430 (2006) 211 – 276

The potential, VSMC , can be written as a displaced harmonic oscillator coupled to a double-well,  VSMC = − 21 aX 2

+

1 ˜ 4 4 cX

+

1 2 2

Q+



2

2 X

2

,

(15)

with a new constant c˜ = c − 2. The displacement of the oscillator Q(0) = −X2 / 2 is proportional to the square of X. The saddle point and minima are at (0, 0) and (± a/c, ˜ −a/c ˜ 2 ), respectively. The barrier height is E = a 2 /4c. ˜ Fig. 5 (top) shows an SMC potential with typical parameters. Two characteristic paths are also given: the intrinsic reaction path (IRP), i.e. the path of steepest descent (in mass-weighted coordinates) from the barrier top and a straight line path connecting the minima. The SMC Hamiltonian depends on four parameters, a, c, , and . The set of parameters can be reduced without restrictions by introducing dimensionless positions (x, ˜ q) ˜ by [56] X = xm x, ˜ (16) Q = xm q, ˜ (17)  where xm = a/c˜ is the minimum position of x, and by dividing the SMC Hamiltonian by 8E = 2a 2 /c. ˜ The new SMC Hamiltonian H˜ SMC reads   2 1 g 2 j2 j2 2  2 2 2 + + ( x ˜ − 1) + ( x ˜ − 1) , (18) ( x ˜ + 1) q ˜ + H˜ SMC = − 2 jx˜ 2 jq˜ 2 8 2 2 where the constant 1/8 is added and the oscillator DOF is shifted by /2 . The new parameters are defined by √ g = h¯ c/a ˜ 2a = h¯ x /8E ,

(19)

√  = / 2a = / x , (20) √ = /2 a c, ˜ (21) √ where x = 2a is the harmonic frequency corresponding to the minima of the original SMC Hamiltonian. It depends on three parameters only. Formally, the parameter g corresponds to the reduced Planck constant of the new SMC Hamiltonian. Moreover, it is related to the dimensionless mass m by m = g −1/2 where h¯ = 1. The remaining parameters,  and , have the same meaning as the original and , respectively. In Ref. [60] a simple method was proposed to determine the SMC parameters. It will be addressed in Appendix A (cf. Section 2.3.1) and an application is given in Section 4.2.3. The case V (2) (X, Q) = XQ

(22)

corresponds to the anti-symmetric mode coupling where the small amplitude coordinate Q is assumed to be antisymmetric with respect to T. The saddle point is at (0, 0); the minima are at Xmin = ± a/c + 2 /c 2 and Qmin = −Xmin / 2 . A typical PES is shown in Fig. 5 (middle). The case V (2) (X, Q)= X2 Q2 is called squeezed coupling. The small √ amplitude coordinate may either be symmetric or antisymmetric. The saddle point and minima are at (0, 0) and (± a/c, 0), respectively. A typical instance of the PES is shown in Fig. 5 (bottom). The PES can be written as an oscillator with X-dependent frequency (X), VSQZ = V (1) (X) + 21 2 (X)Q2 , (X) =



2 + 2X 2 ,

(23) (24)

where for  > 0 the mode is weakened upon approaching the saddle point. The squeezed coupling case is typically realized for out-of-plane modes in, e.g., tropolone [61].

K. Giese et al. / Physics Reports 430 (2006) 211 – 276

221

2.0

Q

0.0

-2.0

-4.0

-6.0 6.0 4.0

Q

2.0 0.0 -2.0 -4.0 -6.0 1.5 1.0

Q

0.5 0.0 -0.5 -1.0 -1.5

-1.5

-1.0

-0.5

0.0

0.5

1.0

1.5

X Fig. 5. Contour plots of instances of PES showing three coupling types (see text). The contour line spacing is 1/5 of the respective barrier height. The IRP (thick black line; schematic) and straight line paths (dashed) are indicated. Top: SMC (Eq. (15), a = c˜ = 1, = 0.25, and  = 3/32). Middle: Eq. (22), a = c = 1, = 0.25, and  = 0.0884. Bottom: Eq. (23), a = c = 1, = 0.25, and  = 0.2. The straight line path coincides with the IRP.

2.3. Reaction surface Hamiltonian 2.3.1. Reaction coordinates The generic Hamiltonians discussed in the previous section describe a large amplitude reaction coordinate coupled to an orthogonal small amplitude harmonic vibration. While the parameters of such types of Hamiltonians can be easily obtained from quantum chemistry calculations as indicated, it is desirable to have a more rigorous method for obtaining

222

K. Giese et al. / Physics Reports 430 (2006) 211 – 276

a Hamiltonian which does not make too many a priori assumptions concerning the dimensionality or functional form of the coupling between large and small amplitude coordinates. There are two general concepts in this respect. First the description of vibrations and in particular of IR spectra in terms of normal mode vibrations [62]. Second, the modeling of chemical reactions by means of a reaction path [63]. Both concepts have been blended in the Miller–Handy–Adams reaction path Hamiltonian [64]:  1 2 1 (ps − kl Qk Pl Bkl (s))2 + V0 (s) + (Pk + 2k (s)Q2k ). (25) H=  2 2 (1 + k Qk Bk,3Nat −6 (s)) 2 k

Here, s is the arc length coordinate along the IRP and Qk (k = 1, . . . , 3Nat − 7) are the orthogonal harmonic DOF. (The ps and Pk are the respective conjugate momenta.) The potential along the reaction coordinate s is given by V0 (s) and the kinetic coupling between the different DOF can be calculated from the knowledge of the normal mode transformation vector (cf. Eq. (1)) for mode k, Lk (s), as well as its derivative with respect to s, Lk (s), via Bkl (s) = Lk (s) · Lk (s) (here Lk denotes the kth column vector of L, see Eq. (1)). 1/2

 2 , a vibrationally adiabatic approximation is reasonable. It was For small curvature, (s) = k Bk,3Nat −6 (s) shown that the effect of small IRP curvature can be accounted for by an effective s-dependent mass eff (s) in the adiabatic Hamiltonian [65], H (ps , s) =

1 p 2 + V0 (s) + E0 (s), 2eff (s) s

(26)

where E0 (s) is the ZPE of the harmonic modes [64]. If the reaction path curvature is too large then the reaction swath has to be taken into account. Moreover, the orthogonal DOF Q may become nonunique for distances away from the reaction path that exceed the radius of curvature. For the case of large reaction path curvature a straight line reaction path Hamiltonian has also been suggested in Ref. [66]. There have been a number of applications of the reaction path concept to the determination of vibrational spectra (e.g. Refs. [67–70]), even including anharmonic corrections for the orthogonal modes [71], or rate coefficient calculations (e.g. Ref. [72]). A large reaction path curvature indicates large couplings among the reactive and the orthogonal DOF. This is reminiscent of the behavior of general heavy-light-heavy reactions [73]. To cope with this situation it has been suggested to use more than one reaction coordinate. In the reaction surface Hamiltonian approach (and variations thereof), instead of the intrinsic reaction coordinate, few internal coordinates (i.e., bond lengths r) serve as large amplitude DOF while the remaining coordinates (Q) are typically optimized such as to minimize the energy for a given point on the reaction surface. Their motion is then again accounted for within harmonic approximation [74–77]. The kinetic energy operator of the aforementioned approach has a rather complicated form making a numerical treatment difficult. Using Wilson’s G-matrix notation one has [75,76]    1 pr Grr GrQ T = (pr , PQ ) , (27) GQr GQQ PQ 2 where pr and PQ are the vectors of momenta for the large and small amplitude coordinate, respectively. An assessment of various approximations to Eq. (27) in terms of the vibrational spectrum has been given recently in Ref. [78]. A much simpler structure of the kinetic energy operator is obtained by adopting a formulation based on Cartesian coordinates such as the CRS Hamiltonian by Ruf and Miller [79]. The CRS Hamiltonian relies on the selection of few atomic Cartesian coordinates that perform large amplitude displacements during the reaction. For instance, for planar HAT systems the in-plane coordinates of the reacting hydrogen are a reasonable choice. The remaining Cartesian DOF are considered to perform only small amplitude displacements and the full PES is approximated by a second order Taylor expansion with respect to the small amplitude DOF. However, there may be significant couplings between small amplitude collective coordinates and large amplitude atomic coordinates just because a complete separation between these two does not necessarily yield the most compact representation of the PES. For the intramolecular HAT in tropolone, the situation is illustrated in Fig. 6. Panel (a) shows an overlay of the minimum geometry (Cs symmetry) and saddle point geometry (C2v symmetry). A significant motion of the heavy atoms upon approaching the saddle point geometry is visible; especially the oxygen-oxygen distance is shortened. The chemical bonding has significantly changed between the two considered geometries, i.e., the PES seen by the two oxygens differs also quite strongly. Thus,

K. Giese et al. / Physics Reports 430 (2006) 211 – 276

H

H

O C

H O

O

O

Z

C

C

223

O

O

C

C

C

Y (a)

(b)

(c)

Fig. 6. The OCCO–H fragment of tropolone. (a) Overlay of minimum XR (filled circles) and saddle point XTS (open circles) geometry. Coordinate origin is the center of mass. (b) Atom displacement corresponding to direction w1 (anti-symmetric). The displacement of the remaining atoms is small. (c) Same for direction w2 (symmetric).

the choice of the atomic (y, z) coordinates leads to a strong coupling of large amplitude (H-atom) and small amplitude coordinates. All geometries X (mass-weighted) of the molecule span a N-dimensional vector space, where N is the total number of DOF. Thus, instead of Cartesian coordinates of individual atoms one may alternatively choose any orthonormal set of N-dimensional vectors to describe large amplitude changes of the molecular geometry. A suitable choice of such vectors has been proposed by Takada and Nakamura [60]. Let XR be the 3 · Nat dimensional mass-weighted vector corresponding to the geometry of one minimum (denoted as right). The center of mass is at the origin. Likewise, let XL and XTS correspond to the other minimum (denoted as left) and the saddle point geometry, respectively. The left (XL ) and right (XR ) minimum is only unique up to an arbitrary rotation. This arbitrariness is removed by rotating XL around XR in order to minimize the distance |XR − XL | [60]. Correspondingly, the saddle point geometry XTS is rotated around the center geometry, XC = (XR + XL )/2, in order to minimize |XC − XTS |. Then, the two vectors, XR − XL , |XR − XL | XC − XTS , w2 = |XC − XTS | w1 =

(28) (29)

span the so-called reaction plane [80]. The vector w1 corresponds to the direct tunneling direction, i.e., the straight line connecting both minima and the center geometry; the vector w2 points from the saddle point towards the center geometry. Thus, for a symmetric double well system the displacement vector XR − XTS is partitioned into the two symmetry components, i.e. w1 and w2 for the anti-symmetric and symmetric motion, respectively. The two directions w1 and w2 are orthogonal, w1 · w2 = 0, because they transform according to different irreducible representations. In panels (b) and (c) of Fig. 6 we show the vectors w1 and w2 defined in Eqs. (28) and (29) for the case of tropolone. Notice, that w1 [panel (b)] is similar to the y coordinate, while w2 [panel (c)] describes a concerted motion of the hydrogen together with the oxygens. Notice that the reaction plane concept was used previously by Yagi et al. [80] as a guide for the relevant region of configuration space in their full-dimensional treatment of malonaldehyde. Moreover, the idea was implicitly used by other authors in the context of an anharmonic expansion of the PES for the same molecule [81]. The relevance of the reaction plane definition can be appreciated in the context of reaction paths proposed by Garrett and Truhlar [82]. They considered a family of paths parameterized by a parameter  (0  1), which for the case of two equivalent minima coincides with the two turning points on the IRP, ±s0 . The reaction path is then given by ˜ X(s; ) = (1 − )X(s) + [X(s0 ) + (s − s0 )(X(s0 ) − X(−s0 ))/2s0 ],

(30)

where X(s) is the coordinate representation of the IRP. For  = 0 this path coincides with the IRP and corresponds to the tunneling path in the small-curvature limit; for  = 1 it is a straight line connecting the two turning points and corresponds to the large-curvature limit. It was shown that the least action principle applied to this path yields

224

K. Giese et al. / Physics Reports 430 (2006) 211 – 276

reasonable semiclassical tunneling splittings for molecules like malonaldehyde [83,84], tropolone [84], and various carboxylic acid dimers [85]. In passing we emphasize, however, that the family of Garrett–Truhlar paths covers a curved but two-dimensional submanifold of the full-dimensional configuration space only. The direction w1 is essentially the straight line part of the family of Garrett–Truhlar paths. Based on these observations, the physical relevance of the reaction plane can be established, when it is possible to show that the IRP (the other extremum of the family of Garrett–Truhlar paths) lies approximately in the reaction plane. To be specific, consider the projection of the IRP onto the reaction plane, X(s) = [(X(s) − XTS ) · w1 ]w1 + [(X(s) − XTS ) · w2 ]w2 + XTS .

(31)

The difference between the IRP and its projection can be expressed by the root mean squared atomic displacement, 1 (x(s) − x(s))2 , (32) (s) = √ Nat where lower case letters refer to non-mass-weighted coordinates. The smaller (s) the closer is the IRP to the reaction plane. An example for a comparison between a full-dimensional and a reaction plane IRP is given in Fig. 7 for tropolone. Apparently, (s) is rather small indicating the relevance of the reaction plane definition for this case. 2.3.2. A full-dimensional all-Cartesian Hamiltonian In the following we will give a derivation of a full-dimensional reaction surface Hamiltonian in terms of generic reaction coordinates v1 and v2 with corresponding vectors v1 and v2 . As discussed in Section 2.3.1 there are at least two reasonable choices for reaction coordinates (cf. Fig. 6). First, one may choose (v1 , v2 ) to equal certain atomic coordinates, for instance, the position (y, z) of the reactive hydrogen atom in tropolone. This is identical to the choice of Ruf and Miller [79]. Second, one may choose (v1 , v2 ) to equal the reaction plane coordinates (w1 , w2 ) discussed in Section 2.3.1 [86]. Note, that in the following there is no inherent restriction as far as the number of reactive DOFs is concerned. The reaction coordinates (v1 , v2 ) are assumed to be orthogonal to the 3 infinitesimal translational vectors as well as to the three infinitesimal rotational vectors corresponding to the saddle point geometry XTS . For two generic reaction coordinates there is a (N − 2)-dimensional subspace of the space of all geometries X, where N is the number of DOF, i.e., N =3 Nat for a molecule with Nat atoms. A basis for this subspace can be obtained, e.g., by considering a projection of the full-dimensional Hessian at the saddle point geometry onto the subspace. Let U (X) be the full-dimensional PES, then the full-dimensional Hessian at the saddle point is given by j2 U (f) Kij (XTS ) = . (33) jXi jXj XTS The projector onto the subspace is given by I − P, where I is the identity matrix and P is the projector onto the space spanned by v1 and v2 , P = v1 v1T + v2 v2T .

(34)

A diagonalization of the projected matrix, (I − P)K(f) (XTS )(I − P),

(35)

yields N eigenvectors; the set of eigenvectors is denoted as {v1 , v2 , e1 , . . . , eN−2 }. There are eight eigenvectors with vanishing eigenvalue: two correspond to the reaction coordinate vectors v1 and v2 , the others correspond to the 3 infinitesimal translations, and to the three infinitesimal rotations with respect to the saddle point geometry. For convenience, the set of eigenvectors is ordered such, that the 6 translational/rotational vectors are eN−7 , . . . , eN−2 . The remaining N − 8 vectors ej are the normal modes of the saddle point with respect to the subspace. Coordinates corresponding to the eigenvectors ej including rotation and translation are denoted as Qj . For systems with two symmetrically equivalent minima, the eigenvectors ej and the associated coordinates Qj are either symmetric or anti-symmetric with respect to the molecular symmetry transformation T.

K. Giese et al. / Physics Reports 430 (2006) 211 – 276

225

Fig. 7. Comparison between a reaction plane and a full-dimensional IRP for the intramolecular HAT in tropolone using a DFT/B3LYP (6−31+G(d)) PES. (a) Contour plot of the potential cut along the reaction plane spanned by w1 and w2 (in a0 amu1/2 ). The contour line spacing is 500 cm−1 and the maximum contour line is at 6000 cm−1 . The thick black line is the projection of the IRP onto the reaction plane X(s), and the barrier height is 2161 cm−1 . (b) Root mean squared atomic displacements Eq. (32) for the full-dimensional IRP geometries on the approximate PES (solid line) and for the quantum chemically determined PES (filled squares). For more details, see Ref. [86].

Any geometry X can be expanded in terms of the set of eigenvectors discussed so far, X = XTS + v1 v1 + v2 v2 +

N−2 

Q j ej ,

(36)

j =1

where the saddle point geometry serves as a reference, i.e., the coordinates (v1 , v2 , Q1 , . . . , QN−2 ) describe displacements from the saddle point geometry. In the CRS framework small displacements X = X − X0 with respect to a so-called reaction surface X0 are considered, X =

N−2  j =1

(ref)

(Qj − Qj

)ej .

(37)

226

K. Giese et al. / Physics Reports 430 (2006) 211 – 276

The displacements are orthogonal to the space spanned by v1 and v2 and the reaction surface is parameterized by the (ref) (ref) reaction coordinates, i.e., Qj = Qj (v1 , v2 ). The Taylor expansion of the PES with respect to the reaction surface reads, U (X) = U (X0 ) + G(f) (X0 ) · X + 21 XT K(f) (X0 )X + · · · ,

(38)

(f)

where Gi = jU/jXi is the full dimensional gradient and the Hessian is given in Eq. (33). For small displacements of the Q a truncation of the Taylor series after the harmonic term is reasonable. The truncated PES is denoted V in the following. The three translational DOF can be discarded without loss of generality; the three rotational DOF can be discarded under the assumption, that the infinitesimal rotational vectors do only slightly change on the plane spanned by the reaction coordinates. The PES is then expressed as, V = V (v1 , v2 , Q1 , . . . , QN−8 ). The potential values, gradients, and the Hessians on the reaction surface can be obtained by standard quantum chemistry calculations. 2.3.3. Fixed reference case (ref) The choice Qj ≡ constj corresponds to a reaction surface that is independent of (v1 , v2 ) and was termed the fixed reference case [79]. (Note that this kind of reference has to be distinguished from the choice of XTS as reference in (ref) Eq. (36).) Formally, due to this choice, the reaction surface becomes a plane. The special case Qj ≡ 0 is important for the formulation of the Cartesian reaction plane Hamiltonian; in this case, the PES reads V (v1 , v2 , Q) = U (v1 , v2 ) −



Fi (v1 , v2 )Qi +

i

1 Kij (v1 , v2 )Qi Qj , 2

(39)

ij

where U (v1 , v2 ) ≡ U (X0 ), and the force Fi acting on mode Qi and the Hessian Kij are related to the full-dimensional quantities by, respectively, Fi (v1 , v2 ) = −G(f) (X0 ) · ei , Kij (v1 , v2 ) = eiT K(f) (X0 )ej .

(40) (41)

2.3.4. Flexible reference case (ref) (ref) With a flexible reference, Qj = Qj (v1 , v2 ), X0 describes a generally non-planar surface. This facilitates the selection of the region of the configuration space that is relevant for the reaction. It could be obtained, for instance, by partial geometry optimization for a each given point (v1 , v2 ). Inserting Eq. (37) into Eq. (38) yields the PES of the CRS Hamiltonian with respect to the flexible reference, V (v1 , v2 , Q) = U˜ (v1 , v2 ) −

 i

1 F˜i (v1 , v2 )Qi + Kij (v1 , v2 )Qi Qj , 2

(42)

ij

with abbreviations, U˜ = U +



(ref)

Fj Q j

+

j

F˜i = Fi +



1 (ref) (ref) Kij Qi Qj , 2

(43)

ij

(ref)

Kij Qj

,

(44)

j

where all quantities are considered to be (v1 , v2 )-dependent and evaluated with respect to X0 . 2.3.5. Selection of relevant modes The PES of the CRS Hamiltonian of the previous Section and the exact PES agree only up to second order terms with respect to Q, but the CRS Hamiltonian is full-dimensional. For an efficient numerical treatment, it is necessary to select certain modes out of the set {Qj } that are especially important and to introduce a further approximation for

K. Giese et al. / Physics Reports 430 (2006) 211 – 276

227

the remaining modes. The selection procedure depends on the choice that was made for the reaction coordinates. The goal is to formulate a reduced dimensional CRS Hamiltonian (h¯ = 1), HCRS = −

n  1 j2 1 j2 1 j2 − − + V (v1 , v2 , q), 2 jv12 2 jv22 2 jqk2 k=1

(45)

where q = (q1 , . . . , qn ) is a set of relevant coordinates. The choice of Q-modes is formally arbitrary, i.e., one may switch to another -equivalent- set of modes by a linear transformation, Q = AQ,

(46)

where A is an orthogonal transformation matrix that does not mix modes with different symmetry. It is assumed in the following, that the vector Q can be expressed as   q , (47) Q = S where the vector S = (S1 , . . . , Sn ) with N = n + n + 8 accounts for the set of irrelevant modes. The set of reaction coordinates plus the relevant modes, {v1 , v2 , q1 , . . . , qn }, is called model coordinates; the irrelevant modes are denoted as spectator modes. 2.3.6. Atomic reaction coordinates The simplest choice for reaction coordinates are atomic positions as discussed in Section 2.3.2. One may, for instance, identify (v1 , v2 ) with the position of the reactive hydrogen atom in the plane of the molecule [79]. For atomic reaction coordinates the choice of a flexible reference is preferred since otherwise for large displacement from the reference geometry one would obtain large forces indicating a substantial deviation from some IRP. Thus, the starting point is Eq. (42) which can be written in a more convenient form as follows: let Q(0) = Q(0) (v1 , v2 ) be such that V (v1 , v2 , Q(0) ) is minimal among all Q for any fixed value of the reaction coordinates, i.e., Q(0) satisfies,  (0) Kij Qj . (48) F˜i = j

With this definition the linear and harmonic term of the PES can be joined in a displaced harmonic term, 1  (0) (0) Kij (x, y)(Qi − Qi )(Qj − Qj ), V (v1 , v2 , Q) = U˜ (v1 , v2 ) − E Q (v1 , v2 ) + 2

(49)

i,j

where the so-called reorganization energy, 1  (0) (0) Kij (v1 , v2 )Qi Qj , E Q (v1 , v2 ) = 2

(50)

i,j

is introduced. In analogy to Eq. (47), the vector Q(0) can be divided according to  (0)  q (0) . Q = S(0)

(51)

Formally, the coupling between model coordinates and spectator modes is assumed to be negligible. Then, the number of actual DOF of the PES can be reduced by setting, S ≡ S(0) (v1 , v2 ).

(52)

This choice is preferred among all other possible choices, because even when there is a small coupling that is not strictly negligible, the energetics of the reduced dimensional PES is still equivalent to the full-dimensional one, i.e., the barrier height is the same.

228

K. Giese et al. / Physics Reports 430 (2006) 211 – 276

2.3.7. Collective reaction plane coordinates The minima and the transition state geometries are points on the reaction plane. Therefore, the choice of a fixed reference, Q(ref) ≡ 0, is reasonable. This implies, X0 (w1 , w2 ) = XTS + w1 w1 + w2 w2 . Let Q denote the projection operator onto the n-dimensional subspace spanned by the relevant modes q. Diagonalization of the matrix, (I − P − Q)K(f) |X0 (I − P − Q),

(53)

for each fixed value of (w1 , w2 ) yields n nonvanishing eigenvalues 2i (w1 , w2 ) for the spectator modes as a function of the reaction plane coordinates. The PES in terms of the model coordinates is given by Eq. (39) with S ≡ 0 (cf. Eq. (47)) V (w1 , w2 , q) = U (w1 , w2 ) −

n  i



n n 1  1  h Fi (w1 , w2 )qi + Kij (w1 , w2 )qi qj + ¯ i (w1 , w2 ). 2 2 ij

(54)

i

The last term of this expression corresponds to the ZPE of the spectator modes. The ad hoc inclusion of this term into the reduced-dimensional PES is motivated by intuition; the change of frequencies of many modes may lead to a contribution to the potential, that would have been neglected otherwise. More rigorously, such a term appears within an adiabatic separation of the relevant coordinates and spectator modes. Here the spectator modes are relaxed to their equilibrium positions which depends on the relevant coordinates. This introduces not only the ZPE into the potential but it also gives a contribution to the kinetic energy [87,88]. A recent systematic study in the context of tunneling splittings in the formic acid dimer can be found in Ref. [89]. For all three extremal points on the reduced PES, q = 0 holds. Furthermore, neglecting the ZPE term, the energetics of the reduced PES, e.g., the barrier height, is identical to the full PES by definition of the reaction plane. This is a very important feature of the present formulation, which distinguishes it from the choice of atomic reaction coordinates according to Ref. [79], because one does not need to include relaxed spectator modes in order to yield the same PES energetics as in the full dimensional case. Previously, Yagi et al. [80] used the reaction plane as a starting manifold for the generation of points for a fulldimensional treatment of malonaldehyde. In their approach, the global PES is spanned by a modified Sheppard scheme. The main advantage of the present approximate method, however, is the fact that the PES can be written as a sum of products of function of the reaction coordinates (w1 , w2 ) times the remaining coordinates q. This form makes its use in numerical propagation schemes for the solution of the Schrödinger equation which are based on a product representation of the wave function very efficient (see Section 3.1). Finally, we note that the dipole moment surface of the reduced model can be approximately expressed in a similar manner. For numerical convenience, only the first derivative with respect to the reaction plane can be employed, while the dipole function is treated numerically exact on the reaction plane. The approximation reads  j(full)  (full)  (w1 , w2 , q) =  (w1 , w2 , q)|q=0 + (55) qk , jqk k

(full)

where  vector.

q=0

is the full-dimensional dipole function and  = X, Y, Z is the Cartesian component of the dipole moment

2.3.8. Reduced normal modes In this section generic reaction coordinates (including atomic and reaction plane coordinates) are considered. The model coordinates—except v1 and v2 —are arbitrary in the sense, that any linear combination of them will yield exactly the same results. There are three sets of coordinates, however, that are unique, these are—as in the full-dimensional system—the normal modes of the reduced (n + 2)-dimensional system at the two symmetrically related minima and at the transition state. These normal modes are called reduced normal modes. By neglecting changes of the zero-point energy, which are usually small if all important coordinates are included into the relevant part of the Hamiltonian, one can compute, e.g., the reduced normal modes of the right minimum geometry by diagonalizing the Hessian RK(f) (XR )R,

(56)

K. Giese et al. / Physics Reports 430 (2006) 211 – 276

229

where R is the projector onto the space spanned by the model coordinates (v1 , v2 , q1 , . . . , qn ). The corresponding n+2 ([n+2]D) with k = 1, . . . , n + 2 in the following. Recall that the vectors v1 and v2 reduced normal modes are denoted Yk are both N-dimensional. Likewise, the vectors ej are N-dimensional according to their definition by diagonalization of Eq. (35). Thus, the (n+2)-dimensional model coordinates define a (n+2)-dimensional reduced space that is embedded in the full N-dimensional space. While for setting up the Hamiltonian it is most convenient to use model coordinates, all physically relevant quantities should be related to the unique reduced normal modes. Additionally, to analyze the connection of the reduced dimensional model with the full-dimensional system, one has to investigate overlaps of reduced normal modes with full normal modes. This task can be achieved by defining a N × (n + 2)-matrix B that transforms from the reduced space (v1 , v2 , q1 , . . . , qn ) back into the full N-dimensional Cartesian space: B = (v1 v2 q1 . . . qn ),

(57)

where the vectors v1 , etc., i.e., those vectors that span the so-called reduced space, constitute the columns of the matrix. This matrix always exists by definition and the property, BT B = 1, holds because the constituting vectors are orthonormal. Since molecular IR spectra are normally analyzed in terms of harmonic transition frequencies corresponding to normal modes it is necessary to establish a link between the modes obtained in harmonic approximation and for the reduced dimensional model. In other words, a reasonable reduced model should be characterized by overlaps (f) with certain full normal modes being close to one. Let Yj with j = 1, . . . , N − 6 be the full normal modes corresponding to the right configuration XR , and let Yk be the reduced normal modes corresponding to the right mini(min) (min) mum configuration (v1 , v2 , q = 0). Then the projection of reduced normal mode k onto full normal mode j is given by (f)

pj k = (Yj )T BYk .

(58)

An application of this procedure will be discussed in Section 3.3. 3. Coherent quantum dynamics 3.1. Mean-field vs. multiconfiguration methods In the following we will address the question how the time-dependent nuclear Schrödinger equation, ih¯

j (q; t) = H (q; t). jt

(59)

can be solved efficiently for a multidimensional HB system with general coordinates q. The exact basis set representation of a multidimensional wave function reads (q; t) =

n1 

1 =1

···

nN 

N =1

(N) C1 ,...,N (t)(1) 1 (q1 ) · · · N (qN ),

(60)

(k)

where C is a time-dependent coefficient matrix and k (qk ) are time-independent basis functions. The coefficient matrix scales like O(n¯ N ) with respect to the number of DOF N, where n¯ is a typical number of basis functions for one dimension. Unfortunately, this exponential scaling hampers any attempt of a direct propagation of a multidimensional (say, with N > 4) wave function. In order to motivate a solution to avoid this dimensionality bottleneck consider a reaction coordinate X coupled to a harmonic vibrational mode Q where the interaction is written in the form of Eq. (14), i.e. q = (X, Q). Suppose the total wave function is of Hartree-product form (X, Q; t) = (X; t)(Q; t)

(61)

230

K. Giese et al. / Physics Reports 430 (2006) 211 – 276

with the single particle functions (X; t) and (Q; t). Equations of motion for the single particle functions can be derived by the Dirac–Frenkel variational principle [90] |H − ih¯

j | = 0. jt

(62)

One obtains, neglecting unimportant phase factors, ih¯

j (X; t) = [TX + VSCF (X; t)](X; t) jt

(63)

ih¯

j (Q; t) = [TQ + VSCF (Q; t)](Q; t). jt

(64)

and

The coupling between the two coordinates in this time-dependent self-consistent field (TDSCF) approach is realized by the so-called mean-field potentials, e.g.,

(1) VSCF (X; t) = V (X) + dQ∗ (Q; t)V (2) (X, Q)(Q; t) (65) and similar for VSCF (Q; t). While this approach works well for many molecular applications [91], it is likely to perform rather poorly for the case of HAT as first pointed out by Makri and Miller [92]. The reason for a possible failure is obvious: If the single particle functions are delocalized along the respective coordinates all details of the interaction potential are washed out by the integration in Eq. (65). Recently, this has been nicely illustrated in a model simulation of PT along a water chain [93]. However, even in cases where one deals with a single minimum potential only, the limited flexibility of the Hartree ansatz may surface in the improper description of IVR (see, Section 3.2). To solve this issue several multiconfiguration schemes have been proposed [92,94], but the breakthrough especially in terms of general applicability came only with MCTDH approach developed by Meyer and coworkers [95–97]. Here the ansatz for the wave function reads (q; t) =

n1 

1 =1

···

nN 

A1 ,...,N (t)1 ,...,N (q; t),

(66)

N =1

where a single configuration is given by a Hartree product of N time-dependent single particle functions, (N) 1 ,...,N (q; t) = (1) 1 (q1 ; t) · · · N (qN ; t).

(67)

The integers nj in Eq. (66) refer to the number of single particle functions corresponding to a certain DOF qj . Equations (j ) of motion for the coefficient matrix A1 ,...,N (t) and the single particle functions j (qj ; t) can also be derived by the Dirac–Frenkel variational principle, Eq. (62) [95]. For the numerical integration of the equations of motion the single particle function are often expressed in a discrete variable representation grid [98,99]. For more details in particular concerning the implementation within the Heidelberg MCTDH code [100], see Ref. [96]. The time-dependent basis (i.e., a set of time-dependent single particle functions) can follow the wave packet during the propagation, making the ansatz Eqs. (66)–(67) more efficient than using the same number of time-independent basis functions (cf. Fig. 8). The efficiency can be further increased by combining certain modes [101], which in fact had been suggested earlier to incorporate correlations in the framework of a Hartree product ansatz [102]. Let q(ck ) with 1k r denote a subset of the full coordinate vector q. There are r N mutually disjoined subsets in total. A modified single configuration is now given by a product of functions corresponding to the individual subsets of coordinates, (c1 ) (cr ) ˜ 1 ,...,r (q; t) = (1)  ; t) · · · (r) ; t). 1 (q r (q

(68)

˜ and N → r in Eq. (66). For convenience the basis functions The full MCTDH ansatz is obtained by replacing  →  appearing in Eq. (68) are also denoted single particle functions.

K. Giese et al. / Physics Reports 430 (2006) 211 – 276

231

ϕ (2)

t=0

t>0

q2

Ψ

time q1 ϕ (2)

ϕ (1)

evolution Ψ

ϕ (1) Fig. 8. Illustration of the MCTDH method. The initial wave packet (t = 0) moves and spreads during the time evolution. The time-dependent single particle functions (k) can follow the motion of  (only one single particle function per DOF is shown).

Provided that a separation between slow and fast coordinates has been performed and that the fast coordinates are described by a set of quantum states {|}, the MCTDH wave function can be written as [103] ()

()

 (Q1 , . . . , Qf , t) =

n1  () j1 =1

n

...

f  () jf =1

A

() () () (t) j1 ,...,jf

f 

=1

(,)

j  (Q , t), 

(69)

where the Q are the slow coordinates. Eq. (69) corresponds to the so-called multi-set formulation, where different sets of single particle functions have been defined for each quantum state. There are several important extensions of the MCTDH approach outlined so far. For instance, in the multilayer formulation each single particle function itself is treated recursively as a MCTDH wave function [104]. On the other hand, it has been shown that one can reduce the effort for representing the single particle functions considerably by using parameterized functions such as Gaussians [105]. The MCTDH approach is also capable of diagonalizing a multidimensional Hamiltonian [97]. First, the Hamiltonian is diagonalized within the basis of certain initial single particle functions. This yields expansion coefficients corresponding to each eigenstate. The coefficients corresponding to the desired eigenstate are used to build the quantities appearing in the equations of motion for the single particle functions. Then, a small integration step is performed for imaginary time (t → i) and a new optimized set of single particle function is obtained. The procedure, which is called improved relaxation, is repeated until convergence is achieved [97]. It is important to notice that an efficient implementation of this MCTDH scheme requires to have at hand a Hamiltonian which is of product form in the different coordinates. This way mean-field integrals like Eq. (65) factorize and can be straightforwardly calculated. There is a scheme for transferring general potentials into product form [106], but this is not feasible for multidimensional situations where also the very task of calculating the PES becomes impossible. On the other hand, the CRS Hamiltonian, Eq. (45), is of product form with respect to most coordinates which makes it particularly well-suited for studying multidimensional HB quantum dynamics using the MCTDH approach. 3.2. Laser-driven dynamics and intramolecular vibrational energy redistribution 3.2.1. The case of a single minimum potential As a first example let us consider the HB dynamics triggered by excitation of either the OH or the OD bond in phthalic acid monomethylester (PMME) shown in Fig. 9. PMME has been the first case for which coherent vibrational dynamics of a HB on a subpicosecond time scale could be observed using ultrafast 130 fs IR pump–probe spectroscopy [107,108]. Monitoring the decay of an initially excited OD stretching vibrational state an oscillatory component corresponding to a 100 cm−1 vibration was found in the signal which survived for about 1.5 ps (decay time about 500 fs). The excited

232

K. Giese et al. / Physics Reports 430 (2006) 211 – 276

Fig. 9. Upper panel: the equilibrium structure of deuterated PMME as obtained using DFT/B3LYP with a 6-31 + G(d, p) basis set [114]. Lower panel: the potential along the O1-D coordinate (with respect to the equilibrium value) is shown together with the lowest eigenfunctions.

state itself decayed with a time constant of T1 = 400 fs, whereas the overall vibrational cooling time was about 20 ps. Subsequently, a similar behavior was found for excitation of the OH vibration in the undeuterated species where merely the T1 time decreased to 200 fs [109]. Furthermore, photon echo measurements revealed a homogeneous dephasing time shorter than 100 fs [110]. The rather short T1 time seems to be a more general feature of this type of HBs, e.g., a 200 fs decay was observed for deuterated 2-(2 -hydroxophenyl)benzothiazole in Ref. [111]. In order to unravel the pathways for vibrational energy flow two-color pump–probe experiments have been performed for PMME as well which will be discussed in Section 6.3 [112,113]. ˚ i.e. it should be classified as a medium-strong The HB is nonlinear (see, Fig. 9) with the O–O distance being 2.56 A, HB. Apart from an isoenergetic enantiomer [115] there is only a single minimum PES which is rather anharmonic as can be seen from Fig. 9 where the potential is shown as a function of the OD bond distance for the deuterated species. The anharmonic shift for the OD = 1 → 2 transition amounts to −180 cm−1 in this one-dimensional PES cut. Apart from this shift there is a strong Fermi-resonance with the OD bending overtone [116]. For the normal species, on the other hand, one observes a more significant substructure and broadening in the IR spectrum as seen in Fig. 31 below. Here, the OH = 0 → 1 transition is calculated at 3036 cm−1 and the anharmonic shift for the OH = 1 → 2 transition is −430 cm−1 [112]. Notice, however, that the assignment of the overtone transition must be taken with care, because these zero-order states are usually rather diluted among several eigenstates. Overall, these stationary spectra are indicating the underlying multidimensional HB dynamics. A global view on the dynamics can be obtained by combining a CRS description with a TDSCF propagation. Given a dipole moment surface , the interaction with the external laser field can be included via the Hamiltonian: Hfield (t) = −E(t) cos(t),

(70)

where  is the carrier frequency of the pulse and E(t) is its envelope usually assumed to be of Gaussian E(t) = E0 exp(−2 ln(2)t 2 /2 )

(71)

or sin2 shape (0 t ) E(t) = E0 sin2 (t /),

(72)

where E0 is the amplitude and  characterizes the pulse duration. Notice that for simplicity we assumed that the field is directed along the dipole moment vector. Using the OH or OD-bond length as the reaction coordinate an essentially full-dimensional description of the coupling to the molecular scaffold has been given in Ref. [114]. There it was shown that indeed the periodic modulations of the

K. Giese et al. / Physics Reports 430 (2006) 211 – 276

233

Fig. 10. PMME pump–probe signal (in arbitrary units) calculated on the basis of a full-dimensional CRS-Hamiltonian and a TDSCF propagation. The pump pulse excites the OH-stretch transition and the probe pulse detects the excited state absorption (pump =3150 cm−1 , probe =2750 cm−1 ,  = 330 fs (Eq. (72))). The oscillations are due to a low-frequency HB mode which modulates the O · · · O distance. For more details see Ref. [114].

Fig. 11. MCTDH coordinate expectation value with respect to the equilibrium configuration for a 19D-CRS model of the HB dynamics in PMME-D: √ (a) OD-bond distance, (b) Low-frequency HB mode (in units of aB a.m.u.). (For more details, see Ref. [117].)

pump–probe signal [108] can be traced back to the coupling between the OH-stretching vibration and a low-frequency vibration influencing the strength of the HB (cf. Fig. 3). Taking a simplified definition of the pump–probe signal as being equal to the energy absorbed by a weak probe pulse starting from the state prepared by the pump pulse one obtains the signal shown in Fig. 10 [114]. The oscillatory signal reflects one aspect of the experimental results, i.e. the periodic modulation of the HB strength by coherent excitation of an anharmonically coupled low-frequency mode. However, in contrast to the experiment it does not decay. In principle this observation could have at least two reasons: First, the TDSCF wave function does not provide enough flexibility for intramolecular energy flow. Second, the interaction with the solvent is a key for the understanding of the relaxation process. In order to scrutinize the first point a reduced dimensional model has been defined in Refs. [115,117] by considering only the most important substrate modes. The selection was based on the reorganization energies (cf. Section 2.3) as well as on the TDSCF dynamics [114]. For this model wave packet propagations for the reaction coordinate and 8 and 10 harmonic modes being treated within the MCTDH and the TDSCF scheme, respectively, have been performed in Refs. [115,117]. In Fig. 11 we show the expectation value of the OD bond distance as well as of the low-frequency HB coordinate for an excitation with a  = 300 fs sin2 -shaped pulse, Eq. (72). The laser-driving at the resonant frequency is reflected in the fast oscillations in panel (a). Superimposed one notices in particular a slow modulation with a period of about 500 fs. This is a direct consequence of the anharmonic coupling to the low-frequency HB mode as can be seen from panel (b). The maxima of Qlow  coincide with the minima

234

K. Giese et al. / Physics Reports 430 (2006) 211 – 276

Fig. 12. Expectation value of the uncoupled energy of the OD bond, EOD /hc, for the 19D model of PMME-D (see, Fig. 11). (a) MCTDH propagation, (b) TDSCF approximation. (For more details, see Ref. [117].)

Fig. 13. Enol-keto tautomerism of 2-hydroxybenzaldimine compounds.

of the envelope in panel (a). This implies that the elongation along Qlow leads to a compression of the HB and thus a smaller amplitude of the OD vibration. The energy exchange between the OD bond and the substrate modes can be monitored by calculating the expectation value of the reaction coordinate part of the Hamiltonian neglecting any coupling to the orthogonal normal modes. In Fig. 12 we show EOD  for the MCTDH (panel (a)) and the TDSCF (panel (b)) cases. Apart from the rather similar rapid oscillations we notice that in the MCTDH case EOD  is (on average) slightly decreasing after the initial rise until about 0.25 ps while in the TDSCF case it stays constant. This can be interpreted as the onset of energy randomization, i.e. the irreversible (at least in the considered time interval) energy flow from the OD-vibration into the rest of the molecule. This effect is only captured by the MCTDH approach. This finding clearly highlights the importance of incorporating multiconfiguration effects even for single minimum systems at moderate energies. From the decay of EOD  after about 0.25 ps (averaging over the rapid oscillations) we can estimate the associated time constant which is of the order of about 20 ps. We will return to this example in Section 6.3 where the role of the solvent for understanding the experimentally observed [108] much more rapid decay will be addressed. 3.2.2. The case of an asymmetric double minimum potential The anharmonicity of single minimum systems such as PMME is still moderate since no HAT takes place which would require a change in the conjugation pattern of the molecule as shown, for instance, for the enol–keto tautomerism in 2-hydroxybenzaldimine compounds in Fig. 13. Here, not only the incorporation of large amplitude motion but also accounting for multiconfiguration effects in the quantum dynamics is mandatory. Let us consider salicylaldimine (SA) as the simplest example (R&H) for which a number of quantum chemical investigations exist [32,118–121]. At the DFT/B3LYP(6-31+G(d,p)) level the barrier height is 1940 cm−1 and the enol–keto asymmetry is 1308 cm−1 [121].

K. Giese et al. / Physics Reports 430 (2006) 211 – 276

235

Fig. 14. PES for the 2D planar HAT in salicylaldimine, where x and y give the position of the H atom with respect to the center of mass and the moments of inertia axes of the enol configuration as well as a relaxed molecular scaffold (cf. Section 2.3.4). Also shown are the probability densities for the first three states (a–c) and for the first state which is localized in the keto minimum (d). (For more details, see Ref. [122].)

Fig. 15. Displacement vectors of the normal modes which enter the 7D model of HAT in SA–H/D [121,122]. The harmonic frequencies are: 4 = 336, 6 = 452, 14 = 861, 26 = 1333, 30 = 1513 cm−1 .

Restricting ourselves to a planar model, a possible choice for the large amplitude reaction coordinates are the in-plane x and y positions of the moving H atom with respect to the center of mass and the moments of inertia axes of the enol configuration (cf. Section 2.3.6). The 2D PES for the bare H motion shown in Fig. 14a apparently is rather anharmonic. Moreover, the reorganization energy for the harmonic modes (cf. Eq. (50)) of 3573 and 3112 cm−1 at the transition state and keto configuration, respectively, points to a substantial coupling. In Refs. [121,122] a 7D model has been suggested for SA–H and SA–D, respectively, which includes the five strongest coupled modes shown in Fig. 15. Notice that these modes can be roughly classified as being either of the symmetric (4 , 14 ), i.e. promoting, or the antisymmetric (6 , 26 , 30 ), i.e. reorganization, coupling type. Fig. 14 also shows some of the lowest eigenstates of the reaction coordinates for normal mode coordinates fixed at the enol minimum. They have been obtained from the solution of the Schrödinger equation for the wave functions  (x = (x, y)) ≡ x|, [Tx + U (x)] (x) = E  (x).

(73)

236

K. Giese et al. / Physics Reports 430 (2006) 211 – 276

Fig. 16. Absorption spectrum, Eq. (77), of the 7D model of SA–H (black) and SA–D (gray) calculated using the MCTDH approach. The dashed lines show the pulse spectra used for the laser-driven dynamics in Fig. 17. (For more details, see also Ref. [123].)

These states form a crude diabatic basis, {|}, [38] which can be used for a diabatic representation of the CRS Hamiltonian:  H diab = 21 P2 + [H (Q) + H (Q)(1 −  )]| |. (74) ,

Here, we introduced the diabatic PES H (Q) = E − F Q + 21 QK Q,

(75)

and the coupling between the diabatic states H (Q) = −F Q + 21 QK Q.

(76)

This representation has at least two advantages: First, it offers a convenient means to analyze the spectrum and associated wave packet dynamics in terms of the (zero-order) diabatic states. For instance, Fig. 14 shows the ground state as well as states dominated by bending and stretching fundamental excitations. Second, since the diabatic states are rather localized, the grid on which the first and second derivatives of the PES have to be known is limited. An overall characterization of the model can be obtained by calculating the linear IR absorption spectrum which is conveniently done via the Fourier transformation of the dipole–dipole autocorrelation function [38]:

∞ 4 nmol  I () = Re dtW (t)eit 0 |[ (t),  (0)]|0 , (77) 3h¯ c =x,y,z 0 where |0  is the ground state, nmol is the number density of molecules,  (t) is the dipole operator in the Heisenberg picture with spatial direction , and W (t) is a certain window function to account for the limited time interval on which the correlation function is known. In Fig. 16 we show the IR absorption spectrum for SA–H and the deuterated SA–D which has been evaluated according to Eq. (77) using the MCTDH method [123] in the multi-set formulation, Eq. (69). SA–H shows a broad and structured band covering the range from 2200–2800 cm−1 . Upon deuteration (SA–D) the substructure almost disappears and the band shifts into the range from 1800–2050 cm−1 . Experimental results are available only for N-phenylsalicylaldimine (R–Ph) in CCl4 solution [124]. Apart from the solvent-induced broadening one observes two bands around 2600 and 2750 cm−1 what could be taken as a confirmation of the magnitude of the calculated red-shift. The difference between SA–H and SA–D in Fig. 16 can be traced to the reduced resonant couplings between the reaction

K. Giese et al. / Physics Reports 430 (2006) 211 – 276

237

Fig. 17. Energy change for the 7D model of SA–H (left) and SA–D (right) after ultrafast excitation of the OH and OD dominated absorption band (see pulse spectra in Fig. 16). The data correspond to (from top at t = 750 fs): Exy (gray), EQ26 , EQ14 , EQ30 , EQ6 , and EQ4 (left, OH), and EQ26 , EQ14 , Exy (gray), EQ30 , EQ4 , and EQ6 (right, OD). (For more details, see also Ref. [123].)

Fig. 18. Tautomerism of tropolone.

coordinates and the substrate modes in the deuterated case. This effect is nicely reflected in the dynamics upon OH and OD stretching excitation. In Fig. 17 we compare the energy expectation values of the reaction coordinate and the uncoupled harmonic modes for SA–H and SA–D after excitation with a 260 fs sin2 -shaped pulse. The IVR dynamics of SA–H is characterized by energy randomization on a time scale of about 700 fs, a number which correlates reasonably with available experimental data [124]. Furthermore, none of the harmonic modes seems to play a particularly important role. For the deuterated species no such fast decay is observed. Instead there is a pronounced energy exchange with modes 14 and 26 which persists at least up to 5 ps [122]. Hence we can conclude that for the case of a rather anharmonic double minimum potential, i.e. for typical HAT systems, a subpicosecond decay of the OH-stretching vibration can be caused solely by intramolecular anharmonic interactions. 3.3. Mode-specific tunneling splittings An important example for multidimensional HB dynamics is the dependence of tunneling splittings on the excitation of particular vibrational modes in the molecule. This feature has been observed spectroscopically for electronic ground [125] and excited [30,34] states. As discussed in the introduction mode-specific tunneling splittings are just another manifestation of, e.g., the promoting role environmental motions can have in HAT and PT reactions. Before addressing approximate methods for calculating tunneling splittings in Section 4 we will discuss the quantum route to IR spectra and related tunneling splittings for the exemplary tautomerization of tropolone, Fig. 18. Tropolone is among the most intensively studied examples for mode-specific tunneling splitting although this effect

238

K. Giese et al. / Physics Reports 430 (2006) 211 – 276

has been observed in other system such as, e.g., the formic acid dimer [29]. Studying the vibrational progressions in vibronic spectra of tropolone, Alves and Hollas [126,127] as well as the groups of Mikami [34] and Sekiya [30] found a pronounced mode-specificity in the S1 -state tunneling. For the electronic ground state Tanaka et al. measured the lowest tunneling doublet to be separated by 0 = 0.974 cm−1 [128]. Redington and coworkers [125,129] attributed a doublet at 741 and 751 cm−1 to the fundamental transition of what was called the nascent skeletal tunneling coordinate [130]. In the region of the OH-stretch fundamental transition a doublet is observed at 3102 and 3121 cm−1 [129], i.e. the tunnel splitting is as small as 1 ≈ 20 cm−1 . A similar conclusion concerning the magnitude of 1 was reached in Ref. [131], although based on a different assignment. The controversial assignments can be taken as an indication that the identification of tunneling doublets in a rather structured spectrum is not straightforward. Given the complex line shape of the OH-stretching band one may also ask whether it is reasonable to assign a transition to a particular normal mode whose definition is meaningful in the vicinity of the minima on the PES only. General aspects of mode selective tunneling were addressed theoretically by using two-dimensional model Hamiltonians adjusted to mimic tropolone [60,61]. On the other hand, Redington put forward a two-dimensional adiabatic model which included the OH-stretch and the skeletal tunneling coordinate corresponding to a vibration at around 750 cm−1 [130]. The remaining modes contributed only their ZPE. This model has been rather illuminating and after properly choosing a tunneling path it could even reproduce the experimental splittings quantitatively. In contrast, in Ref. [86] using a full-dimensional analysis of the PES in the vicinity of the minimum we have shown that the intrinsic reaction path approaches the PES minimum not along this skeletal tunneling coordinate but along a combination of the two lowest frequency in-plane modes (∼ 360 cm−1 ). Overall, in view of the highly structured IR spectrum in the OH-stretch region one may argue whether a two-dimensional description is rigorously justified, not to mention the adiabatic approximation. The CRS approach offers a way to investigate the nature of the IR spectrum based on a multidimensional quantum treatment [86]. The reaction plane coordinates have already been discussed in Fig. 6. In principle, all remaining orthogonal DOF could be incorporated in harmonic approximation into a 39-dimensional model. In practice, the numerical effort for the solution of the full-dimensional nuclear Schrödinger equation can be considerable and a reduction of the dimensionality along the lines discussed in Section 2.3.5 is advisable. Here, one is guided by the process which one intends to model. If one puts the focus on a proper description of the low-frequency part of the IR spectrum associated with the reaction coordinate as well as on the OH-stretching region, overlap and resonance criteria suggest an overall twelve-dimensional model [132]. In Fig. 19 we show the reduced normal modes of the right configuration for the 12D model and in Table 1 we give the overlap between the reduced and the full normal modes according to Eq. (58) as well as with the reaction plane coordinates. These overlaps convey some information concerning the target to which the model coordinates have (12D) (12D) (12D) been adapted. There are three modes, Y1 , Y9 , Y12 , with a comparatively large overlap with the reaction (12D) (f) (12D) almost coincides with the OH stretch full mode Y1 . The reduced mode Y9 has plane. The reduced mode Y1 (f) a significant overlap with the full mode Y15 , which has OH bend character and shows the largest IR intensity. Finally, (12D) (f) (f) (f) the reduced mode Y12 has lowest frequency and shows significant overlaps with the full modes Y33 , Y36 , and Y37 , i.e. it represents those modes associated with the approach of the minimum along the IRP as discussed above. Thus (12D) one could call Y12 the reaction mode of the 12D model. This is in contrast to the assignment of Ref. [130] which (12D) (f) would correspond to mode Y10 , i.e. a mode which has a perfect overlap with the full mode Y28 , but at the same time has only a very small overlap with the reaction plane coordinates. The ground state tunneling splitting for this model is obtained as 0 = 2.7 cm−1 , what is about three times the (exp) experimental value of 0 = 0.974 cm−1 [128], which indicates most likely an underestimation of the barrier by density functional theory [86] (cf. Fig. 7). An overview of the IR spectrum, Eq. (77), is given in Fig. 20. There are three distinct regions in the spectrum, marked A, B, and C. The 12D model Hamiltonian has been designed to account for regions A and C only. First, let us focus on region A which contains the information about the reaction mode. Here, it turns out that the resolution which can be reasonably obtained from a time propagation of the autocorrelation function is not sufficient for a detailed analysis of the spectrum. Low-lying eigenstates can also be obtained by using the improved relaxation method as discussed in Section 3.1. There are also other MCTDH based diagonalization schemes as used, for instance, by Coutinho-Neto et al. [133] to determine the ground state tunneling splitting of malonaldehyde within a full-dimensional model potential.

K. Giese et al. / Physics Reports 430 (2006) 211 – 276

239

Fig. 19. The reduced normal modes of the 12D model of tropolone. Frequencies are indicated.

Table 1 Comparison of normal modes of the full-dimensional system with the reduced normal modes of the 12D model for tropolone Full No. 1 7 8

12D

(f) k /2 c

Inten.

Overlap

Mode

(12D) /2 c k

w1

3329

133.7

0.97

Y(12D) 1

3227

−0.64

0.31

0.87

Y(12D) 2

1662

−0.01

0.00

0.89

Y(12D) 3

1655

−0.15

0.04

1593

0.06

−0.01

1674 1663

202.4 2.82

w2

9

1617

104.6

0.91

Y(12D) 4

10

1539

124.5

0.94

Y(12D) 5

1535

0.04

0.00

1500

−0.26

0.07

11

1522

141.9

0.85

Y(12D) 6

12

1480

212.0

0.83

Y(12D) 7

1461

0.17

−0.07

1400

−0.04

0.00

13

1454

27.0

0.59

Y(12D) 8

14 15

1349 1327

57.7 227.5

0.65 0.72

Y(12D) 9

1271

0.60

−0.11

753

−0.11

0.03

28

752

14.5

1.00

Y(12D) 10

30

694

7.2

1.00

Y(12D) 11

695

−0.03

0.00

33 36

449 375

15.2 5.3

0.41 0.70

Y(12D) 12

392

0.32

0.94

37

364

5.2

0.56

Frequencies k (in cm−1 ), IR intensities (in km/mol), overlaps of the reduced normal modes with full normal modes [cf. Eq. (58)], and overlaps (12D) (12D) of the reduced normal modes with the reaction plane directions, |w1 · Yk | (“w1 ”) and |w2 · Yk | (“w2 ”), are given.

240

K. Giese et al. / Physics Reports 430 (2006) 211 – 276

3.00 B

Intensity (arb. units)

2.50 2.00 1.50 1.00

C

A

0.50 0.00 250

750

1250

1750

2250

2750

 / 2 c (cm-1) Fig. 20. Overview of the computed IR spectra of the 12D model of Trn(OH) at T = 5 K. Important regions are marked by A–C. A 3ps MCTDH propagation has been used together with a Gaussian window function W (t) ∼ exp[−t 2 /2 ],  = 1 ps. Table 2 Excitation energies (εj ) and IR intensities for the 12D [“H-12”] and a 4D [“H-4”] model of tropolone [86] Final state f

IR intensity (km/mol)

+ 0 → f

j /hc (cm−1 )

− 0 → f

H-12

H-4

H-12

H-4

H-12

H-4

1) 121+ (4+ 1 1) 12− (4−

357

354

7.1

6.7

14.0

10.6

371

372

8.8

5.2

6.2

6.4

2) 122+ (4+ 2) 122− (4−

697

690

0.6

0.6

0.7

0.9

737

744

0.0

0.1

0.5

0.3

101+ (31+ ) 101− (31− )

752

752

0.7

1.7

14.6

20.7

750

749

15.6

21.5

0.6

1.1

Parenthesized labeling refers to the 4D model.

The results (i.e., eigenenergies and IR intensities) of a still another diagonalization procedure outlined in Appendix B are shown in Table 2 (indicated as “H-12”). The symmetry of the states is either gerade (“+”) or ungerade (“−”). The labeling of states is nk± , where n is the reduced normal number and k is the number of quanta in that mode. The eigenenergies are given with respect to the gerade ground state level + 0 . A pronounced mode selectivity is found for the (12D) k states 12 corresponding to mode Y12 : the splitting increases to 14 cm−1 for the fundamental (121 ) and to 41 cm−1 for the first overtone (122 ). This is not surprising since according to Table 1 this mode most closely resembles the (12D) reaction mode, i.e., that mode, which vibrates in the direction of the IRP. On the other hand, excitation of mode Y10 which corresponds to the skeletal tunneling coordinate of Ref. [130] influences the tunneling splitting only marginally. The origin of the mode-selectivity can also be scrutinized by inspecting the respective wave functions. In Ref. [86] a four-dimensional model of this reaction had been studied which focussed particularly on the low-frequency region (compare the respective energies and IR intensities in Table 2). In Fig. 21 we show 2D representations of the 4D (12D) (12D) and Y10 in the 12D model. Apparently, the densities for states corresponding to the excitation of modes Y12 mode-selectivity can be traced to the symmetry of the mode, i.e. there is an enhancement of the tunneling splitting for (12D) (12D) the symmetric promoting-type mode Y12 and no effect for the asymmetric coupling mode Y10 . Next, we focus on the OH-stretching region C whose separate contributions from the gerade and ungerade initial ground state are shown enlarged in Fig. 22. The major components of the spectrum can be used to assign a tunneling splitting of 1 ≈ (20±8) cm−1 . Although the absolute position of the OH-stretch band is shifted by about 250 cm−1 with respect to the experimental value [129,131], the order of magnitude of the splitting is rather well reproduced [86].

K. Giese et al. / Physics Reports 430 (2006) 211 – 276

241

Fig. 21. Densities ||2 of selected states of the 4D model of tropolone [86]. (The other DOF in this 2D representation have been integrated out.) (12D) The excitation of a symmetric normal mode (approx. Y12 ) has substantial structure in the reaction plane. In the right panel qas is that reduced (12D)

model mode which correspond to an asymmetric normal mode (approx. Y10

0.50

I () (arb.units)

0.40

∆1

).

Ψ+0 Ψ0-

0.30 0.20 0.10 0.00

2600 2700 2800 2900 3000 3100  / 2 c (cm-1)

Fig. 22. The OH stretch region computed for the 12D model of tropolone. Transitions with gerade and ungerade initial state are shown separately. Lines are drawn to guide the eye. The OH stretch tunneling splitting 1 is indicated.

It should be noted, that the OH-stretching excitation is energetically above the reaction barrier. Thus, this suggests that the splitting is due to dynamical tunneling [134]. A similar conclusion was drawn for the HAT in 3,7-dichlorotropolone in Ref. [135]. The complex structure of the IR spectrum in Fig. 22 deserves a more detailed study to pinpoint the contribution of specific modes of the 12D model. This can be done by studying the dynamics underlying this broad band. One concept (0,R) which is frequently used in the theory of IVR is that of zero-order states [136]. In the present case, the states n (0,L) (n ) corresponding to the harmonic approximation to the right (left) minimum form a complete basis (cf. Appendix B). Moreover, in the absence of anharmonic couplings, these states would be the eigenstates of the Hamiltonian which is the motivation for their common use in the interpretation of IR spectra. These states are eigenstates of the right/left-hand normal mode Hamiltonian    1 (R/L) (R/L) H0 h , (78) = + ¯ j nˆ j 2 j

(R/L)

where nˆ j (R) H0

is the number operator of harmonic mode j located on the right/left-hand side. The expectation value of

is defined by the real part of the half-space integral [137]

(R) (R) H0 R = Re dr ∗ H0 . x>0

(79)

242

K. Giese et al. / Physics Reports 430 (2006) 211 – 276

Fig. 23. IVR dynamics for the 12D model of tropolone. The initial wave packet is the harmonic oscillator state corresponding to a singly excited OH stretch. (a) Energy flow for specific modes (indicated) as well as sums of the harmonic ER+L (t) and anharmonic EA (t) contributions. (b) Same as in (a), but for remaining modes (from bottom to top): mode no. 11, 2, 8, 10, 4, 3, 6, 5, and 7. The curves are shifted in steps of 50 cm−1 (except for curves corresponding to mode no. 11, 2, and 8); the energy is E = 0 at t = 0 for all curves.

(L)

Equivalently, H0 L can be defined for the left half-space. The expectation value of the full Hamiltonian H can be divided according to, (R)

(L)

H  = H0 R + H0 L + HA ,

(80)

where the Hamiltonian HA contains the anharmonic part of the full Hamiltonian. In Eq. (78) together with Eq. (79) one can identify the energy expectation value without ZPE of right mode j, as (R) (R) Eˆ j R = h ¯ j nˆ j R

(81)

(L) and analogous for the left mode energy Eˆ j L . In Fig. 23 we show the dynamics after initial preparation of the single excited local OH stretch mode localized at (0,R) the right-hand side: 11 . The time evolution of this state with respect to the full Hamiltonian H leads to preferential energy flow into those modes, that are coupled to the local OH stretch mode. Thus, the energy flow pattern yields information about the coupling among modes. From Fig. 23a we notice that during the first 100 fs the energy quickly redistributes from the OH-stretch to the OH-bend and to the low-frequency reaction mode. Furthermore, there is some partitioning between energy in the harmonic and anharmonic part of the Hamiltonian which changes only slowly after the first 100 fs. Notice that there is only a 250 cm−1 energy contribution flowing from the anharmonic part into the

K. Giese et al. / Physics Reports 430 (2006) 211 – 276

243

harmonic part. Compared to the initial harmonic energy this is less than 8%. This observation justifies the use of harmonic modes for the present analysis. As far as the total 12D model is concerned, Fig. 23b demonstrates that the energy is rapidly redistributed over all modes during the first 500 fs. Closer analysis of the initial decay rates in Fig. 23 allows the conclusion that there are essentially two decay channels, first via the OH bending vibration and second via all remaining modes. In other words, the detailed form of spectral pattern as well as the rather small tunneling splitting in Fig. 22 is a consequence of the coupling between all modes of the 12D model. A similar fast IVR dynamics has been observed in Ref. [137] for a 14D CRS-model of 3,7-dichlorotropolone. Together with the results for the asymmetric double minimum potential this suggest, that rapid subpicoscond IVR after OH-stretching excitation could be a common feature of medium–strong intramolecular HBs having a double minimum PES. 4. Semiclassical and classical trajectory-based methods 4.1. Semiclassical approximation for the calculation of spectra Traditionally, semiclassical spectra and related tunneling splittings are addressed by stationary WKB theory [138] whose multidimensional formulation [56] starts with the N-dimensional wavefunction (q) expressed as   i (q) = exp S(q) , (82) h¯ and the complex-valued function S(q) is expanded in a power series with respect to h, ¯ S(q) = S0 (q) +

h¯ S1 (q) + i

 2 h¯ S2 (q) + · · · . i

(83)

Inserting this ansatz into the stationary Schrödinger equation and keeping only the lowest order terms in h, ¯ one obtains the Hamilton–Jacobi equation for the action function S0 and the transport equation for the amplitude S1 . Notice that since S0 is in general a complex-valued function the configuration space can be conveniently divided into regions where S0 is real- (R), imaginary- (I ), or complex-valued (C) [56]. For the calculation of tunneling splittings Herring’s formula [139]  

(sc) 2 N−1 (L) j (R) (R) j (L)  = h¯ , (84) d q   −  jX jX X=0 has been widely used [56,57,140] where (R/L) are the right- and left-localized wavefunctions and X is the reaction coordinate. In the following we will give a trajectory-based picture of the tunneling splitting. In this respect, it is important to note that according to Eq. (84) a priori the tunneling splitting is not related to a tunneling path, but to the properties of the wave function along a symmetry surface (or line in 2D) . Moreover, as was pointed out by Benderskii et al. [57], if one interprets tunneling in terms of classical trajectories by means of semiclassical theory, such trajectories need to be followed only up to . Especially, these trajectories may never reach the opposite well, but, nevertheless, their contribution can be significant [141]. Let us focus on tunneling in a two-dimensional SMC-potential. Takada and Nakamura [56] showed that there are two extreme cases of tunneling in a SMC potential: (a) tunneling proceeds through the C region or (b) tunneling proceeds through the I region. Case (a) and (b) correspond to small and large coupling, respectively, between the reaction and harmonic coordinate. The situation is sketched in Fig. 24. In case (a) the saddle point falls into the C region. No trajectories can be defined in that region and the action function is complex; especially, it is complex along the symmetry line. In the I region the semiclassical wave function can be expressed as   1 ˜ (X, Q) = exp − SI (X, Q) − S(X, Q) , (85) h¯

244

K. Giese et al. / Physics Reports 430 (2006) 211 – 276

Q

C

I # IRP

C

(a) X Q IRP LSLA

#

instanton

I

C

rT HU

C (b) X Fig. 24. Tunneling mechanisms (schematic) for the SMC-PES. The saddle point () and the IRP (dash–dotted) is indicated. Thick solid and dashed lines correspond to caustics of the allowed and forbidden regions, respectively. (a) Small coupling, and (b) intermediate coupling. The instanton trajectories (dotted) and the trajectories for the locally separable linear approximation (thin solid) are indicated.

˜ SI satisfies the Hamilton–Jacobi equation for the upside-down PES, E − V (r) and where S0 ≡ iSI and S1 ≡ −S. p(r) = ∇SI (r) is the momentum of a classical trajectory on the upside-down PES traversing position r. Here, vectors r denote a position in the 2D space (X, Q). It is assumed that there is a single point rs = (0, Qs ) on the symmetry line  where the exponent is stationary, jSI (X, Q) (86) = 0, jQ rs i.e., there is a particular trajectory with momentum perpendicular to the symmetry line at rs . In the vicinity of rs on , the action SI can be expanded into a second-order Taylor series with respect to Q, SI (X = 0, Q) ≈ SI (0, Qs ) + 0.5 SI (0, Qs )(Q − Qs )2 , where SI is the second derivative with respect to Q. The tunneling splitting can be determined by inserting Eq. (85) into Herring’s formula for the SMC case, evaluating the derivative with respect to X, inserting the truncated Taylor expansion for , and performing the Gaussian integral analytically in the semiclassical limit h¯ → 0 [56],  = B exp{−2SI (0, Qs )/h}, ¯ with the prefactor  B = h¯



SI (0, Qs )

˜

e−2S(0,qs ) 2|px (0, Qs )|,

(87)

(88)

K. Giese et al. / Physics Reports 430 (2006) 211 – 276

245

where px = jSI /jX is the momentum in X-direction. Thus, in this approach, the tunneling splitting is dominated by a particular trajectory. This trajectory is defined by the condition Eq. (86) and the boundary conditions of the field of trajectories corresponding to the I region. Takada and Nakamura [56] derived analytical expressions for the tunneling splitting by making simplifying assumptions: Separability and linearity of the PES in between the hyperbolic umbillic point (denoted HU in Fig. 24) and the take-off point rT . Then, the tunneling splitting according to this locally separable linear approximation is given by Eq. (87), where SI (0, Qs ) is the action along the path denoted LSLA in Fig. 24. The locally separable linear approximation can be applied for ground and excited states as long as the invariant torus exists. An application to multidimensional systems with N > 2 seems not to be feasible because it involves a search for the take-off point in multidimensional space. Under the assumption that the harmonic approximation is valid for the wavefunction in the allowed region one can extend the I region to the minima such that R and C regions vanish and the PES is inverted, −V (X, Q). Consider the right minimum: a field of trajectories that emanate from the right minimum in the infinite past can be defined [59]. Only an extreme tunneling trajectory is considered, that satisfies Eq. (86) at certain point Qs . Such trajectory starts-off from the symmetry line with momentum perpendicular to the symmetry line and approaches the hill along the weakest mode. For a SMC-type PES, Benderskii et al. [59] found a splitting formula similar to Eq. (87), where SI ≡ Sinst is the action along half of the extreme tunneling trajectory. The extreme tunneling trajectory is called the instanton. This approach was also applied to other types of potentials [57,58] as well as to a full-dimensional model of malonaldehyde [142]. The instanton theory was also generalized to treat low vibrationally excited states [59,143]. The crucial point in multidimensional instanton calculations are the determination of the instanton trajectory as well as the calculation of the prefactor. For the latter, an adiabatic separation of different DOFs is often done [144–147].A related approximate approach was introduced by Tautermann et al. [83–85] which combined the search for a Garrett–Truhlar type tunneling path (cf. Eq. (30)) with a prefactor coming from a one-dimensional theory by Garg [148]. Recently, Mil’nikov and Nakamura [149,150] gave a rigorous numerical implementation of instanton theory that is applicable to multidimensional systems. The performance of the new method was demonstrated, e.g., for malonaldehyde [151] and the formic acid dimer [152]. In particular it was shown in Ref. [152] that approximate instanton theories may predict tunneling splittings for vibrationally excited states even in cases where the applicability of instanton theory is not justified. Real-time propagations of the semiclassical wave function have been intensively discussed as an alternative to the stationary WKB approach. In particular the initial value representation of the Herman–Kluk [153] semiclassical propagator [154] attracted considerable interest (for applications in the present context, see, e.g. [155–160]). Using the filter diagonalization method for the spectral analysis of semiclassical correlation functions (see Ref. [161] for a recent review), spectra including tunneling splittings could be obtained for simple SMC-type model PES [156,160]. However, the anharmonicity of ab initio PES seems to render this approach inapplicable as far as high-resolution information, such as tunneling splittings, are concerned. This was demonstrated in Ref. [160] for the HAT in 3,7-dichlorotroplone using a CRS-Hamiltonian model. 4.2. Classical trajectory approach for tunneling 4.2.1. The Makri–Miller model (MM) The dynamics of large systems is often successfully described by classical mechanics. Focussing on HB dynamics and in particular on HAT or PT the immediate question arises how quantum effects can be incorporated into a classical trajectory-based method. Tunneling paths are widely discussed in the context of generalized transition state theories [63,162,163]. Usually adiabatic separability of the intrinsic reaction coordinate and the remaining orthogonal DOF is assumed. In Section 2.3.1 it was already mentioned that the effect of small IRP curvature can be accounted for by means of an effective mass. Avoiding the introduction of the effective mass leads to tunneling paths that lie on the concave side of the IRP and cross the symmetry surface at a position displaced from the saddle point. This effect was already investigated by Marcus and Coltrin and it was termed corner cutting [164]. The same effect can be found, e.g., for the instanton path (cf. Section 4.1). For large IRP curvature the small curvature tunneling approximation breaks down and straight-line paths are favored leading to the large curvature tunneling or sudden approximation. The straight-line paths go from a turning point on one

246

K. Giese et al. / Physics Reports 430 (2006) 211 – 276

side of the IRP to the other side, i.e., they cross a nonadiabatic region, the so-called reaction swath. However, already Marcus and Coltrin [164] pointed out that the “true” tunneling path should satisfy the least action principle of classical mechanics. The formal justification is given by the multidimensional WKB theory (cf. Section 4.1): the solution of the Hamilton-Jacobi equations in the I region is given by classical trajectories propagated on the inverted PES. Still another family of paths was suggested by Garrett and Truhlar and has been discussed already in Section 2.3.1. Finally, tunneling paths (maximum probability paths) have been introduced on the basis of the vibrational ground state density [165]. For a general discussion of various approximations in the context of a two-dimensional model potential, see also Refs. [166,167]. To account for tunneling processes in multidimensional trajectory based simulations, Makri and Miller [168] proposed a simple method that is based on the one-dimensional WKB result for the tunneling splitting  in a symmetric doublewell potential [138],   h 1 ¯ = (89) exp −  , h¯

where  is the classical oscillation frequency of the particle in each of the wells and  is the classical action integral from the left turning point −q0 to the right turning point q0 > 0 (as before q denotes general coordinates),

q0 |p(q)|dq. (90) = −q0

It was noted [168] that according to Eq. (89) one-dimensional tunneling can be interpreted in an intuitive way: a classical particle oscillates back and forth in one well; each time the turning point q0 is encountered, there is a certain (typically small) probability, exp{−/h}, ¯ for tunneling into the other well, with the rate of such events being proportional to the oscillation frequency . For small tunneling probability, the cumulated tunneling probability for the time-interval 0 to t is the sum of the individual probabilities,  (t) = exp{−/h}, (91) ¯ tj  t

where tj are the times at which the turning point q0 is encountered. The function (t) is a staircase function with step-size 2 /. Averaging over an ensemble of trajectories with different phases yields a straight line, the slope of which is proportional to the tunneling splitting,  = 2h¯

d (t). dt

(92)

This equation is the basis for a generalization to the multidimensional case. To this end, it is assumed that tunneling proceeds along straight lines given by a tunneling direction d. A turning point is defined by the condition that the sign of the momentum p projected onto direction d (i.e., p · d) changes from plus to minus. Each time tj when such a classical (cl) turning point qj is encountered, the tunneling probability is computed by [168], j =

max 0



(cl)

d 2[V (qj

(cl)

+ d) − V (qj )],

(93)

where max is the maximum length of the straight-line path for which the square root is real. The cumulated tunneling probability is then given by Eq. (91) with the replacement  → j . The multidimensional (t) is still a staircase but with varying step-size and step-height. However, due to the exponential factor only the smallest j contribute to the sum. The tunneling splitting is obtained from Eq. (92). The ensemble of trajectories is the invariant torus corresponding to a certain eigenstate, i.e., a numerical implementation may use adiabatic switching or normal mode sampling [169]. The quasi-periodicity of the dynamics on the invariant torus ensures (t) to be linear also in the multidimensional case. However, when the initial conditions are generated by an approximate method like normal mode sampling, the sampling error leads to deviations from linearity [169]. For a multidimensional system the choice of the tunneling direction d is nonunique. In the original paper, Makri and Miller proposed to use the shortest straight-line connection between the caustics of the left and right invariant

K. Giese et al. / Physics Reports 430 (2006) 211 – 276

(a)

(b)

247

(c)

Fig. 25. Tunneling directions as proposed by Makri and Miller in Ref. [168]. Tunneling is assumed to occur along straight lines parallel to the given arrows. Thick black lines correspond to caustics of invariant tori. (a) SMC case as well as asymmetric coupling case with large (b) and small (c) displacement.

0.14 0.12 < Γ (t) >

0.10 0.08 0.06 0.04 0.02 0.00

0

5

10

15

20 25 t / fs

30

35

40

45

Fig. 26. Cumulated tunneling probability (t) for an ensemble of 2000 trajectories corresponding to the first-excited OH stretch 12D model of tropolone (cf. Section 3.3). The maximum propagation time was Tmax = 48 fs. The deviation from linearity is due to the approximate nature of the normal mode sampling method.

tori. Their original proposal is illustrated in Fig. 25 for the SMC and asymmetric coupling case. For the SMC case the tunneling direction always coincides with the large amplitude coordinate. For the asymmetric coupling case the tunneling direction not only depends on the coupling strength, but also on the eigenstate of interest. Alternative choices were discussed in the literature [169–173]. It was shown that the Makri–Miller model can be used to estimate tunneling splittings for multidimensional systems. For instance, calculations using 2D model [174] and full-dimensional PES [80] or on-the-fly molecular dynamics [175] were carried out. As an example, we show the cumulated tunneling probability for the first-excited OH stretch of the 12D model of tropolone discussed in Section 3.3. From the slope a tunneling splitting of 32 cm−1 can be determined, in good agreement with the quantum value of about 20 cm−1 (cf. Section 3.3) (Fig. 26). Finally, it was shown that the model can be used to incorporate tunneling into classical molecular dynamics simulations [172]. To this end, a particle that hits a turning point was physically displaced along the straight-line tunneling direction with probability equal to exp{−/h¯ }. 4.2.2. The extended Makri–Miller model (EMM) From the discussion in Section 4.1 it is clear that the straight-line approximation of the Makri–Miller model may break down in a situation sketched in Fig. 24b. Considering an SMC Hamiltonian, ˜  for any set of2 fixed parameters a, c, , and such a situation is met for large enough , since the minima are at (± a/c, ˜ −a/c ˜ ) and the saddle point is at (0, 0). In order to deal with such as situation an extended Makri–Miller (EMM) model has been developed in Ref. [176] that can also approximately account for tunneling through the I region. The following discussion addresses the SMC case, because only in this case a pronounced mode-selectivity (spanning orders of magnitude) can be anticipated. The asymmetric coupling case is complicated by the dependence of the tunneling splitting upon the phase of the wave function along the symmetry line [60]. The squeezed case is not interesting for the present purpose, because straight line paths are well suited, especially the instanton is a straight line in this case (cf. Fig. 5).

248

K. Giese et al. / Physics Reports 430 (2006) 211 – 276

In a one-dimensional system turning points q0 are uniquely defined by the energy E and PES V (q) via condition V (q0 ) = E. In a multidimensional system classically allowed and forbidden regions of configuration space are divided by caustics, that is, a surface (or line in 2D) of points on which trajectories cross from one branch of the momentum function to a neighboring branch. This ensemble is mathematically described by the Lagrange manifold (e.g., invariant torus) on which the branches of the momentum function are defined. For nonseparable potentials one cannot uniquely define the kinetic energy for a particular coordinate, thus it is impossible to obtain the turning points [177]. Therefore, in Ref. [176] the following practical but approximate definition for multidimensional turning points was suggested. Let {dj } be the set of N normal mode vectors corresponding to the right minimum of a double-well system. For a harmonic system, turning points are given by the condition dj · p = 0,

(94)

where p is the momentum vector. In Ref. [176] it was shown that this condition is also reasonable for a modestly anharmonic case. With this prerequisite the algorithm of the extended Makri–Miller model involves the following steps: Initial conditions for an ensemble of trajectories are selected either by adiabatic switching or by normal mode sampling. All the trajectories are initially located either in the right or left well. A trajectory of the ensemble is propagated in the R region (cl) for a certain amount of time T. During this time the trajectory will hit classical turning points qn . These turning points are determined by the condition dj · pt = 0 [Eq. (94)]. (cl) Each time the trajectory hits a classical turning point qn , the so-called parity of motion j of the corresponding direction dj is flipped from +1 to −1. The use of parities of motion allows to introduce real-valued momenta even in the classically forbidden region [178]. More details are given in Appendix C. The trajectory is propagated in the forbidden region according to the equations of motion (C.3–C.4). The condition Eq. (94) is also used to determine nonclassical turning points outside the R-region: When the turning point condition is fulfilled for a trajectory in the forbidden region then another parity k of the corresponding direction is flipped from +1 to −1. Both trajectories in the forbidden region, the one with j = −1 and k = +1 as well as the one with j = −1 and k = −1, will be integrated in this case. When a trajectory, that is propagated in the forbidden region, crosses the symmetry line , the complex (cl) action of that trajectory, Eq. (C.5), is computed along the trajectory from qn to the crossing point on the symmetry (cl) line. There may be several trajectories that emanate from one classical turning point qn . Only that action is kept, that has the smallest imaginary part whose absolute value is denoted Wn . The contributions of all nonclassical trajectories (cl) with smallest action that have emanated from classical turning points qn at time tn are summed up according to    2 (95) exp − Wn . (t) = h¯ t t n

If the initial conditions correspond to an invariant torus then the ensemble average, (t), is a straight line and the tunneling splitting is given by Eq. (92). For the case of a SMC PES the algorithm is illustrated in Fig. 27. Vectors of the 2D space (X, Q) are denoted as r. One may distinguish two important cases: In case (i) a classical trajectory (solid) propagated in the R-region hits a turning point r(cl) “(1)” at the C2 -R caustic. The emanating trajectory (dashed) resides solely in the C2 region and hits the symmetry line. In case (ii) a turning point “(2)” of the C1 -R caustic is touched giving rise to a trajectory in the C1 region. This trajectory hits several nonclassical turning points at the caustics of the C1 and neighboring I regions. For instance, the third non-classical turning point rT of that trajectory is located at the I − C1 caustic. A second parity is flipped at that point giving rise to an additional trajectory that is propagated in the I region. This trajectory then crosses the symmetry line. By this procedure, many trajectories are generated for a single classical turning point r(cl) and several of them will cross the symmetry line. Only the trajectory with smallest imaginary action is kept which is motivated by the least action principle. It is interesting to note that in Ref. [173] it was demonstrated that for the case of HAT in HO− 2 trajectories with one and two parity flips are contributing equally to the tunneling splitting. In the case of small , when the saddle point falls into the C2 region, the parity-based trajectories in that region are expected to have smaller action than those that propagate through the C1 and I region, and vice versa for the case of larger , when the saddle point falls into the I region. The situation in this case is reminiscent of the locally separable linear approximation. Compared to the locally separable linear approximation there are two important differences:

K. Giese et al. / Physics Reports 430 (2006) 211 – 276

249

Q

case (ii) C1 rT

I

(2) case (i) (1) C2

X Fig. 27. Right well of Fig. 24 for an unspecified saddle point. Tunneling trajectories emanating from two classical turning points, (1) and (2), are shown. Cases (i) and (ii) correspond to tunneling through the C and I region, respectively.

(i) the present formulation assumes separability in the C1 region, and (ii) the C1 -I boundary is determined by the turning point condition Eq. (94), i.e., the turning point condition is used for the classical forbidden and allowed region. For the locally separable linear approximation the C1 -I boundary was approximated by a straight-line continuation of the exact C2 -R boundary. The two approximations, have been justified by numerical tests in Refs. [173,176] (see also next section). 4.2.3. Application of the MM and EMM methods The performance of MM and EMM methods has been scrutinized in Ref. [176] on the basis of a generic SMC model potential which allowed for a wide variation of parameters and related tunneling regimes. The findings can be summarized by plotting the tunneling probability Ptunnel (Q) = exp{−2Wn /h} ¯

(96)

vs. position on the symmetry line Q (where X = 0 on ) for the semiclassical methods, EMM and MM, as compared with the density |(0, Q)|2 of the quantum mechanical symmetric ground state along the symmetry line. In Fig. 28 results for an SMC Hamiltonian with parameters g=0.04 and =0.64 and for three values of , i.e. coupling strengths, are shown. The maxima of the EMM and quantum data are normalized to one, i.e., the pre-exponential factor is assumed to vary only slowly along Q compared to the exponential part. The data of the MM is normalized by the normalization constant obtained for the EMM. This is justified, because both methods only yield the exponential part of the wavefunction. In the upper panel ( = 0.16) the saddle point is located at Q = 0.25. Due to corner cutting the maximum of ||2 at Q = 0.16 does not coincide with the saddle point where the energy is lowest along the symmetry line. The scatter plot of the EMM appears to be shifted systematically by about 0.025 to the right, but it agrees quite well with the numerically exact results for Q > 0. For Q < 0 there are two branches of the scatter plot. This can be attributed, respectively, to trajectories that tunnel from the C2 into the I region and to trajectories that tunnel from the C1 into the I region. The scatter plot of the MM shows a single branch that rises to about the maximum of the QM curve. Note that MM and EMM are normalized by the same constant. The maximum of the MM data is only about 4% lower than the maximum of the EMM data, but the MM sampling deviates strongly from the QM results. This is not surprising, because there are no trajectories in the I region. However, for the tunneling splittings only the magnitudes of exp{−2Wn /h} ¯ matter, therefore the magnitude of the MM tunneling splitting agrees with the exact result [176].

250

K. Giese et al. / Physics Reports 430 (2006) 211 – 276

Fig. 28. SMC Hamiltonian with g = 0.04,  = 0.8 and three different values of : The QM density |(0, Q)|2 (dimensionless) along the symmetry line  (Q dimensionless) is compared with scatter plots of the tunneling probability Ptunnel (Q) [cf. Eq. (96)] computed with, respectively the EMM and MM model. The position of the saddle point is indicated by arrows.

In the middle panel of Fig. 28 ( = 0.32) the saddle point is located at Q = 0.5; the maximum of the QM curve is at Q = 0.38. The EMM scatter plot agrees with QM result (systematic shift of 0.025) quite well; the second branch at about Q ≈ 0 almost vanishes. Conversely, the MM results show a strong deviation. Especially, the maximum of the MM is 84% smaller than the EMM maximum. Upon a further increase of the MM scatter values become very small as compared to the EMM values. Exemplary, this is shown in the lower panel of Fig. 28 ( = 0.64). Here, the saddle point is at Q = 1.0 and the QM maximum is at Q = 0.76, respectively. The EMM scatter plot agrees quite well with

K. Giese et al. / Physics Reports 430 (2006) 211 – 276

251

Fig. 29. Geometries of FAD (schematic): (a) minimum (C2h ); and (b) transition state (D2h ) geometry (adapted from Ref. [179]).

the QM results (systematic shift of 0.01). Clearly, in this case only the EMM yields the correct order of magnitude of the tunneling splitting [176,132]. To sum up, since the MM approach does not account for tunneling through the I-region, one observes its breakdown for = 0.32 and 0.64. On the other hand, the approximate EMM method performs reasonably well in such cases [176]. In the following we will focus on a concrete example, that is the double HAT in carboxylic acid dimers. The formic acid dimer (FAD) is the most simple carboxylic acid dimer. The minimum (C2h ) and transition state (D2h ) geometry of gaseous FAD is shown in Fig. 29 [179]. FAD was thoroughly investigated by experimental and theoretical techniques. Recently, the ground state tunneling splitting of FAD (more precisely (DCOOH)2 ) was measured by Madeja and Havenith [29] to be 0.00286 cm−1 . In an early theoretical gas phase study, the structure of FAD was found to keep approximately C2h symmetry along the IRP [180] indicating a concerted transfer as suggested by the transition state geometry. Shida et al. [77] investigated a 3D model of FAD; their tunneling splitting value of 0.004 cm−1 corresponds rather well with the experimental finding. The barrier height was 49.4 kJ/mol and the step wise transfer mechanism was found to contribute about 20% to the tunneling splitting. More recently, Luckhaus obtained a tunneling splitting of 0.0015 cm−1 using a 6D model [89]. Approximate [85,147] as well as a rigorous implementation [152] of instanton theory have been applied to the FAD. The most reliable instanton-based value for the ground state tunneling splitting to date is 0.0038 cm−1 [152]. For elevated temperatures (or excitation energies) the reaction mechanism may change [181,182] or not [183]; a clarification considering also the limit T → 0 was not yet achieved. For the present purpose a concerted HAT for FAD is assumed, i.e., the transition state has D2h symmetry. Thus, the tunneling in FAD is very similar to that of tropolone and related molecules and the main characteristics of the PES can be covered by a SMC-PES. Along the IRP first the monomers approach each other and then the hydrogen transfer takes place [181]. This means, one may vary the IRP curvature by substituting the two hydrogens, H5 and H5 by, e.g., fluorines or phenyl rests. In terms of SMC parameters, the coupling strength can be varied by choosing different substituents. The SMC parameters itself can be determined by the method of Takada and Nakamura [60] as outlined in Appendix A. Various parameters of FAD are reported in Table 3. Results for a number of different levels of quantum chemistry were also reported in Ref. [180]. Among them, barrier heights E deviate by several kJ/mol and no rigorous convergence was achieved. The value in Table 3, E = 47.6 kJ/mol, is assumed to be a reasonable estimate, although a very recent high-level quantum chemistry calculation [DFT/aug-cc-pVTZ with CCSD(T) corrections] points towards a lower barrier of 33.1 kJ/mol [152]. ˚ of the rOO and rCC distances (definition: see The most prominent geometrical feature is the decrease by 0.33 A caption of Table 3) for the saddle point geometry compared to the minimum geometry. Conversely, the rCR distances do not change, and the two monomers approach each other as a whole. There are parameters for two derivatives of FAD in Table 3, too: the fluoro formic acid dimer (FFAD) and the benzoic acid dimer (BAD). In these compounds the two hydrogens, H5 and H5 (cf. Fig. 29) are substituted by either two fluorine atoms (F) or two phenyl rests (C6 H5 ). Both molecules, FFAD and BAD, are planar. The barrier height of FFAD is only slightly increased by 1.5 kJ/mol as compared to FAD, while for BAD it is 7.4 kJ/mol lower. The variation of the geometrical parameters of the carboxylic ring is moderate among the species. Most prominent is the ˚ as compared to FAD (3.91 A) ˚ and BAD (3.89 A). ˚ In all three decreased inter-monomer distance rCC for FFAD (3.80 A) compounds the monomers move as a whole when approaching the saddle point, which can be deduced from the almost unchanged rCR bond length. The dimensionless geometrical parameter  for FFAD (1.19) is only slightly smaller than that for FAD (1.28). For BAD (1.82) it is significantly larger, because the difference between mass-weighted minimum geometry XR and center

252

K. Giese et al. / Physics Reports 430 (2006) 211 – 276

Table 3 Parameters for carboxylic acid dimers [MP2/6 − 31 + G(d)] Molecule

FAD

FFAD

BADa

E (kJ/mol) Ezpe (kJ/mol)

47.6 1.2

49.1 1.5

40.2 1.5

˚ rOO (A) ˚ rOH (A) ˚ rCC (A) ˚ rCR (A)

2.77 (2.44) 1.00 (1.22) 3.91 (3.58) 1.09 (1.09)

2.75 (2.43) 1.00 (1.21) 3.80 (3.50) 1.34 (1.33)

2.69 (2.42) 1.00 (1.21) 3.89 (3.63) 1.49 (1.49)

b

1.28

1.19

1.82

as /2 c (cm−1 ) m /2 c (cm−1 ) ˜ 1 /2 c (cm−1 )

˜

2 /2 c (cm−1 )

3298 203 (497) 1480 471

3350 173 (382) 1181 439

3107 119 (199) 1236 377

˜ 1 /8E g= ˜ 2 / ˜1 = = 2

0.0477 0.318 0.129

0.0371 0.372 0.164

0.0477 0.305 0.169

Parenthesized values correspond to the saddle point geometry. rOO is the O2 .O4 distance, rCC is the C3 .C3 distance, rCR is the C3 –H5 distance, and rOH is the O2 –H1 distance. a B3LYP/6-31 + G(d) geometry. b cf. Eq. (A.1).

geometry XC is primarily the displaced hydrogens, while the difference between saddle point geometry XTS and XC is the displacement of the two monomers. Obviously, the monomers are considerably heavier for BAD leading to larger mass-weighted distances than for FAD and FFAD. ˜ 1 and ˜ 2 , were obtained by the projection of normal mode frequencies onto w1 and w2 , The effective frequencies, ˜ 2 decrease by going respectively. Direction w2 corresponds to an effective inter-monomer vibration. The frequencies from FAD to BAD; this behavior is similar to that of the minimum normal modes of the inter-monomer vibration m ˜ 1 correspond to direction w1 which includes mainly the hydrogen and it is mainly a mass effect. The frequencies ˜ motion. The effective frequency 1 of FFAD and BAD is rather different from the parent compound. No such difference can be anticipated from inspecting the anti-symmetric OH stretch frequencies as . The dimensionless SMC parameters are also given in Table 3 for the three species. The two incorporated DOF are the antisymmetric stretch of the two OH bonds (as ) and the inter-monomer vibration (m ), respectively. They were excluded from the ZPE difference EZPE , which is less than 2 kJ/mol in all cases. The tunneling splitting E obtained for the three methods QM (quantum), MM, and EMM, respectively, is given in Table 4. For FAD, the numerical exact result (3.69 × 10−3 cm−1 ) corresponds well to the experimental value (2.86 × 10−3 cm−1 ), indicating that the effective barrier height is reasonable. The semiclassical methods, MM and EMM, yield the same order of magnitude as the QM method which is satisfactory in view of exponential accuracy. The significant deviation of the experimental BAD value from the theoretical predictions was explained by the fact that the experiment was performed for doped crystals of BAD [85]. Recently, Tautermann et al. [85] determined tunneling splittings of various carboxylic acid dimers using an approximate version of the instanton method. The values are given in the last column of Table 4 and marked as “GT path” (cf. also Section 2.3.1). The quantum chemistry was performed by a hybrid approach using a combination of B3LYP/6-31+G(d) and G2(MP2) theory. The barriers for FAD, FFAD, and BAD are, respectively, 37.0, 36.6, and 32.0 kJ/mol. The present SMC model and the Garrett–Truhlar path results predict that the tunneling splitting of FFAD drops significantly as compared to FAD and BAD. Concerning the SMC model, inspection of Table 4 unveils that the ˜ 1 ; this leads to smaller g, and larger  and as compared to FAD. reason is the decrease of the projected frequency ˜ 2 , which relates to an effective inter-monomer The effect is compensated for in BAD by the substantial decrease of vibration. To summarize, the parameters of FAD (and derivatives) do not lead to predominately tunneling through the I region, and the MM and EMM methods yield a reasonable ground state tunneling splitting. Recently, we have investigated

K. Giese et al. / Physics Reports 430 (2006) 211 – 276

253

Table 4 Tunneling splittings for carboxylic acid dimers as calculated by semiclassical and quantum methods for a SMC model Hamiltonian Molecule

E (cm−1 ) QM

MM

EMM

Exp.

GT patha

FAD

3.69 (−3)

3.20 (−3)

5.45 (−3)

b 2.86 (−3)

2.2 (−3)

FFAD

3.55 (−5)

2.05 (−5)

4.55 (−5)

BAD

1.20 (−3)

0.64 (−3)

1.41 (−3)

5.3 (−4) c 0.22 (0)

2.2 (−3)

a Ref. [85]. b Ref. [29]. c Ref. [184] (seleno-indigo doped BAD crystals).

the mode specific tunneling splitting for the HAT in the hydroperoxyl-anion (HO− 2 ) which is probably the simplest example for multidimensional vibrationally assisted HAT [185]. Here one finds substantial deviations between MM and EMM [173]. 5. Condensed phase Hamiltonian The standard model of condensed phase theory contains a reaction coordinate (s) coupled to a harmonic bath (coordinates {Z }, frequencies {  }) as expressed via the Hamiltonian [38] ⎛  2 ⎞ F (s) 1 ⎝ 2  ⎠. P + 2 Z − H = Ts + V (s) + 2 2

 

(97)

Here F (s) is the coupling between both subsystems. Comparing this expression with the CRS Hamiltonian in Eq. (39) we notice that provided the coordinate dependence of the Hessian is neglected the CRS Hamiltonian actually has the form of the generic system–bath Hamiltonian. While Eq. (97) is only linear in the bath coordinates it can easily be extended to include quadratic couplings as well. In other words the CRS Hamiltonian approach offers the possibility to calculate a system–bath Hamiltonian at least in principle. In practice this approach is restricted to situations where the number of bath DOF is finite such that the Hessian can be calculated. This approach has been taken in Ref. [186] to determine the relaxation rates for the HB dynamics in 3-chlorotropolone (neglecting any solvent) and in particular their microscopic origin in terms of specific anharmonic couplings. In order to derive an alternative more flexible description let us start by considering the scenario of an intramolecular HB in a medium-sized molecule. In most general terms the total system can be separated into an intramolecular and a solvent part. Focussing on the molecule itself the experimental conditions usually will suggest a separation between relevant (i.e. active) and passive DOF. For instance, femtosecond IR spectroscopy ofA–H vibrations highlights processes involving the strongly coupled modes such as the HB vibration on a short time scale only. All other modes, whose coupling to the A–H vibration is weak, will be less important. This does not imply that their effect is negligible. Indeed, together with the solvent DOF they provide a reservoir (or heat bath) for energy release out of the active modes. Furthermore, their fluctuation may cause a modulation of the active DOF’s transition frequencies and by that cause coherence dephasing. The circumstances described so far are also characteristic for other system–bath situations, e.g., electron transfer reactions [38]. The only specification for the situation of HBs so far has been the introduction of two types of reservoirs, i.e. an intra- and an intermolecular one. The motivation for this subdivision also derives from the desire to be able to have a clear assignment of the relaxation processes. Moreover, in particular for the relaxation of higher-frequency modes specific intramolecular reservoir modes might become important and identifying them may suggest to include them into the active system. Notice that for smaller system such as HOD dissolved in D2 O [187,188] it is, of course, possible to include all intramolecular DOF into the relevant system.

254

K. Giese et al. / Physics Reports 430 (2006) 211 – 276

Q

S V (s,Q,Z)

Z

Fig. 30. Separation of the DOFs of an intramolecular HB system in a solvent, into an active (s), a intramolecular reservoir (Q), and an intermolecular (solvent) reservoir (Z).

Thus the total Hamiltonian can be written as [189] H = Hs + Hq + HZ + V (s, Q, Z),

(98)

where s, Q, and Z stand for the active, the intra- and intermolecular reservoir modes, respectively (see. Fig. 30). It should be emphasized that no assumption has been made concerning the type of coordinates so far. If, for instance, normal mode coordinates are used for the molecule, the Hamiltonian will be of the Watson form, Eq. (9). This implies that one has to take care not only of the translational motion but also of the vibration–rotation coupling. These terms are usually taken as part of the bath or the interaction Hamiltonian [187,188]. The system–bath interaction Hamiltonian can be expressed in the standard form [38]  Hs.b ≡ V (s, Q, Z) = K (u) (s)(u) (Q, Z) (99) u

with the operators K (u) and (u) belonging to the system and the bath, respectively. Explicit expressions for these operators can be obtained from a Taylor expansion of the interaction potential V (s, Q, Z). As noted before this approach offers a convenient means for characterizing the interaction in terms of transitions associated with specific DOF. The simplest term is of bilinear form (LL coupling):   (LL) (LL) Hs.b = h si (100) ¯  gi ()Z 

i

(LL)

where we assumed that the solvent bath is harmonic with frequencies  and gi () is the dimensionless coupling strength between the ith system coordinate and the th solvent coordinate. Of course, there is a similar bilinear interaction between si and the intramolecular bath modes {Qj }. While being important for large amplitude system motion, this LL coupling will vanish if the dynamics takes place in the vicinity of a potential minimum where si corresponds to a normal mode vibration. In passing we note that Eq. (100) assumes that the coordinates are dimensionless as well. This can be always achieved by defining an appropriate length scale for a given coordinate [116]. Eq. (100) gives rise to simultaneous one quantum transitions in the system and solvent bath. One of the next terms in the Taylor expansion of Eq. (99) would involve the interaction Hamiltonian of quadratic-linear (QL) form   (QL) (QL) Hs.b = h si2 (101) ¯  gi ()Z 

i

which gives rise to two-quantum transitions in the system and additionally is the leading term responsible for pure dephasing. Again, an interaction ∝ si2 Qj is possible as well, even in the vicinity of a potential minimum. The lowest order term which involves transitions in all three subsystems is given by   (LLL) (LLL) Hs.b = h si ()Qj Z . (102) ¯  gij i



j

K. Giese et al. / Physics Reports 430 (2006) 211 – 276

255

By virtue of this term solvent-assisted transitions become possible. Such processes may be important if, for instance, the deexcitation of the system coordinate is not resonant with the excitation of some intramolecular mode. Of course, other terms of the Taylor expansion are possible and the actual choice will depend on the system at hand and should be guided by experimental observations (an example will be given in Section 6.3). Since the relevant system Hamiltonian will usually be low-dimensional it can be diagonalized, i.e. Hs |a = Ea |a.

(103)

In principle the eigenstates |a will depend on all bath coordinates. In order to make the perturbation of the system states due to the interaction with the bath small, it has been proposed to include the bath-averaged interaction Hamiltonian into the definition of the relevant system Hamiltonian, i.e. [186,188] Hs → Hs + V (s, Q, Z)B .

(104)

This way, only the fluctuations of the system–bath interaction are entering the model Hs.b → Hs.b − V (s, Q, Z)B .

(105)

But there is also a bath perspective: for anharmonic systems a|s|a is likely to be different for different states which changes the system–bath interaction. The interdependence of system and bath is often neglected or even does not occur by virtue of the model Hamiltonian, e.g., assuming a harmonic bath the bath-averaged couplings, Eqs. (100)–(102) would vanish. There are, however, situations where such a treatment is of relevance, e.g. if the system–bath interaction is treated not by a Taylor expansion but within a classical molecular dynamics calculation of the quantum system embedded in a classical bath. Here, it has been proposed to account for the above mentioned interdependence self-consistently within an iteration procedure [188]. The strategy outlined so far can be used to describe nonreactive and reactive dynamics. The separation of the total Hamiltonian, Eq. (98), is tailored for application of perturbation theory with respect to the system–bath interaction. However, often the influence of, e.g., charge fluctuations, of the surrounding solvent or protein is substantial and a nonperturbative treatment is necessary. This holds in particular for PT reactions which come along with a pronounced charge translocation. A nonperturbative approach which does not resort to such a partitioning of the Hamiltonian comes at the expense of neglecting quantum effects in the nuclear motion. In this case the forces driving classical motion can be generated “on the fly”. For small systems they can even be calculated using quantum chemical methods [181,183,190,191]. This has the advantage that the bond-breaking and making processes are correctly accounted for at the given level of quantum chemistry. For large systems such as proteins, however, the majority of the DOF can be described using classical force fields within a hybrid quantum mechanics/molecular mechanics scheme, e.g., by using Car–Parinello molecular dynamics [192–194]. As an alternative to this density functional type of approach, empirical valence bond theory [195] enjoys great popularity as well [196–199]. Here, it is assumed that the ground state PES of the reactive PT system can be obtained from the diagonalization of a nuclear coordinate dependent Hamiltonian matrix whose matrix elements are chosen such as to reproduce, e.g., the energetics of suitable reference reactions. Notice that these hybrid approaches are not restricted to the limit of classical nuclei. Under the assumption that the motion of the proton is much faster than that of the other nuclei, PES for the proton motion can be defined for each instantaneous configuration of the heavy nuclei along a classical trajectory. This corresponds to the second Born–Oppenheimer separation discussed in Section (2.1) [200–204]. Besides the proton motion also a HB coordinate has been incorporated within this strategy in Ref. [205]. 6. Dissipative quantum dynamics of hydrogen bonds 6.1. Lineshape theories and nonlinear response functions The models of anharmonic coupling introduced in Section 2.1 have also been used to discuss the broadening mechanism of IR spectra in the condensed phase. For the simple two-mode model, Eq. (2), two approaches can be distinguished: Using Kubo’s lineshape theory, Rösch and Ratner [206] put forward a model where the dephasing of the AH fundamental transition is a consequence of the direct coupling of the associated dipole moment to the fluctuating field of the surroundings. On the other hand, the dephasing has been modeled by assuming an indirect coupling via

256

K. Giese et al. / Physics Reports 430 (2006) 211 – 276

the damping of the low-frequency mode described either quantum mechanically [207,208] or via a stochastic process [209]. For an increased HB strength additional modes such as the A–H bending vibration have to be taken into account. The interplay between the dephasing coupling and the arising Fermi resonance involving the A–H stretching fundamental and its bending overtone transition has been extensively studied by Henri–Rousseau and coworkers who supplemented the above mentioned two mode quantum model by a Fermi resonance coupling to a bending overtone transition while taking into account the direct dephasing of the bending and stretching modes [46,47]. This way they extended the earlier relaxation-free approach by Witkowski and Wójcik [210] as well as the stochastic description given by Bratos and coworker for the static modulation limit [211,212]. An overview of these developments can be found in Refs. [45,48]. By now it is well-established that the processes contributing to the lineshape can be unravelled by means of timeresolved IR spectroscopy [23]. A unified description of linear and nonlinear IR spectroscopy can be given within the response function formalism in Liouville space [213] (see also Refs. [214,215]). Given the molecule–field interaction Hamiltonian, Eq. (70), the linear polarization can be expressed as [213]

∞ (1) P (t) = dt1 R (1) (t1 )E(t − t1 ) (106) 0

with the response function being defined as R (1) (t1 ) =

i ˆ 1 )[, Weq ]}. tr{G(t h¯

(107)

ˆ ˆ Here, G(t) is the Liouville space Green’s function, i.e. G(t) • =U (t) • U + (t), and Weq the equilibrium statistical operator of the total system. Eq. (107) is essentially a dipole–dipole correlation function (cf. Eq. (77)) and contains all information on the linear IR absorption spectrum. In third order with respect to E(t) the polarization is given by [213]

∞ P (3) (t) = dt3 dt2 dt1 R (3) (t3 , t2 , t1 )E(t − t3 )E(t − t3 − t2 )E(t − t3 − t2 − t1 ) (108) 0

with the response function now being defined as a function of the three time intervals between the interactions  3 i (3) ˆ 3 )[, G(t ˆ 2 )[, G(t ˆ 1 )[, Weq ]]]}. R (t3 , t2 , t1 ) = tr{G(t (109) h¯ This response function contains all information on third-order nonlinear IR experiments such as pump–probe or photonecho methods. Eq. (109) can be read in the following way: the first interaction with the laser field via  promotes the diagonal equilibrium (population) density matrix into an off-diagonal one (coherence). The latter evolves during the interval t1 up to the second interaction which turns the coherences into a nonequilibrium population subsequently propagating during t2 . After the third interaction the system evolves in a coherence during t3 . R (1) (t1 ) and R (3) (t3 , t2 , t1 ) can be simplified if the discussion is restricted to a single transition, say AH = 0 → 1. In this case and after performing a cumulant expansion to second order one obtains for the dipole–dipole correlation function [213] tr{(t)(0)Weq } ≈ |10 |2 exp{−i10 t − g(t)}

(110)

with

t

g(t) = 0

dt 

t

dt  C(t  ),

(111)

0

where we introduced the correlation function of the transition frequency fluctuations 10 (t) 0 C(t) = tr{10 (t)10 (0)Weq }.

(112)

0 is the equilibrium density operator for the state Here, the trace is defined with respect to the remaining DOF and Weq AH = 0. This correlation function can be calculated analytically, e.g., for the shifted harmonic Brownian oscillator model [216]. However, it should be emphasized that this expression is adapted from electronic spectroscopy where the implied Condon approximation is usually fulfilled. This must not be the case for HB, where the coordinate dependence of the dipole moment can be of importance (for a recent study of non-Condon effects see Ref. [217]).

K. Giese et al. / Physics Reports 430 (2006) 211 – 276

257

Thus for the understanding of linear and nonlinear IR spectra energy gap correlation functions, Eq. (112), are required. They contain information about dephasing processes, i.e. pure and energy relaxation type dephasing [38]. In a seminal paper, Oxtoby and coworkers determined dephasing rates from classical molecular dynamics correlation functions and studied the interplay between anharmonicity of the PES and dephasing [218]. Notice that gap fluctuations may also play a role for the determination of relaxation rates [219]. An alternative description is provided by making use of the fact that it is usually only a small subset of relevant DOF in the sense of Section 5 which will interact directly with the external field while the majority of DOF form a bath. Tracing out those bath coordinates and assuming a weak system–bath coupling in the Markovian limit together with an initial factorization of the equilibrium statistical operator into a bath and a relevant system part one obtains, e.g. by using the projection operator formalism, an equation of motion for the reduced density operator (t) = tr B (W ) involving the following Green’s function [38,213]: ˆ s (t) = (t) exp{−iLˆ s t − Rt}. ˆ G

(113)

Here, Lˆ s is the Liouville superoperator for the relevant system defined as 1 Lˆ s • = [Hs , •] h¯

(114)

and Rˆ is the relaxation superoperator which will be discussed in the following section. In order to introduce this approach into Eq. (109) one has to make the additional assumption that the trace over the bath DOF can be performed for each time interval separately [213]. This implies a phase randomization between the different propagation intervals ˆ ˆ s (t). An early application of this approach in the context of linear IR spectra of which allows us to replace G(t) →G HBs has been given in Ref. [220]. Nonlinear IR spectra of a HAT system have been modeled along these lines, for instance, in Ref. [221]. 6.2. Quantum master equation Besides providing a means to calculate the nonlinear IR response functions the Green’s function (113) can also be used to study the dynamics of the reduced density matrix, e.g., triggered by arbitrary laser driving. To this end one needs to solve the following QME [38]: j ˆ = −i[Lˆ s + Lˆ field (t)] − R, jt

(115)

where the Liouville superoperator for the coherent evolution has been introduced in Eq. (114) and Lˆ field (t) • = ˆ [Hfield (t), •]/h. ¯ The relaxation and dephasing due to the system–bath interaction is contained in the superoperator R. Using the explicit form of the system–bath interaction, Eq. (99), it can be expressed in the energy representation given by {|a} (cf. Eq. (103)) as [38]   Rab,cd = ac be,ed (de ) + bd ae,ec (ce ) − ca,bd (db ) − db,ac (ca ). (116) e

e

Here we introduced the damping matrix  (u) (u ) ab,cd () = Kab Kcd Cuu ()

(117)

uu

being defined in terms of the Fourier transform of the bath correlation function as

∞ 1   dteit (u) (t)(u ) (0)B , Cuu () = 2 Re h¯ 0

(118)

where •B denotes the equilibrium expectation value of the bath. The energy representation of the QME is sometimes also called Redfield equation with Rabcd being the Redfield tensor [38,222,223]. The various approximations contained in Eq. (115) have already been mentioned in the previous section. (Notice that Eq. (116) neglects in addition the energy

258

K. Giese et al. / Physics Reports 430 (2006) 211 – 276

renormalization of the spectrum.) There are several approaches which are aimed to go beyond the present limitations, e.g., by including initial correlations or non-Markovian dynamics [224–227] (for an overview, see also Ref. [38]). The actual form of the damping matrix depends on the system–bath interaction operator at hand. For a bilinear (LL) coupling, Eq. (100), one obtains  (LL) (LL) ab,cd () = a|si |bc|sj |d Cij () (119) ij

with (LL)

Cij

(LL)

() = 2 (1 + n())[Jij

(LL)

() − Jij

(−)],

where we defined the spectral density  (LL) (LL) (LL) Jij () = gi ()gj ()( −  ).

(120)

(121)



Within the QME approach the spectral density is the central quantity determining the strength of system–bath interaction at a certain frequency. Eq. (121) has been obtained under the assumption of a harmonic bath which is likely to be of limited use for real solvents. However, within linear response theory the actual fully anharmonic solvent can be mapped onto an effective harmonic solvent. This essentially assumes that there is a uniform distribution of coupling strengths over a large number of solvent modes. The effective spectral density for the coordinate si then reads [228] (see also Ref. [229])   ∞ 2 h ¯ (LL) (eff) Jii () → J () = 2 tanh dt cos(t)C (cl) (t). (122) 2kB T

h¯ 0 Here C (cl) (t) is the classical correlation function of the force, −jHs.b /jsi acting on some system coordinate si whose value is fixed. Note that the use of (classical) force–force correlation functions has a long history in this respect [230–233]. The vibrational relaxation rate between a pair of states can be expressed as [38] kab = 2ab,ba (ab )

(123)

from which the T1 vibrational life-time say of state |a follows as  1 kab . = T1 (a)

(124)

b =a

In other words, this approach which is sometimes also called Landau–Teller method, allows to calculate vibrational life-times microscopically from classical molecular dynamics simulations. In practice, the treatment is not restricted to the LL coupling case, that is, Eqs. (117) and (118) are more general (see below). The use of these expressions   often comes along with the replacement (u) (t)(u ) (0)B → (u) (t)(u ) (0)class [187,188]. Notice, however, that this introduces some ambiguity since classical and quantum correlation functions have different properties, e.g., the classical function will not lead to detailed balance [38]. This is usually accounted for by introducing certain (problem specific) correction factors [234,235]. Alternatively to a classical description of the spectral density for vibrational energy relaxation one can use model (LL) and a new spectral functions. Suppose that the coupling strength can be expressed in terms of the single parameter gi (LL) (LL) (LL) (LL) (LL) () is introduced as being equal for all system coordinates, i.e. Jij = gi gj j . A common form density j of the spectral density is [38] 2 j (LL) () = () exp(−/(LL) ), c (LL)

(125)

is some cutoff parameter. As shown in Ref. [236] Eq. (125) can rather accurately reproduce calculated where c spectral densities for HB systems in simple systems.

K. Giese et al. / Physics Reports 430 (2006) 211 – 276

259

In terms of the PES expansion of Section 5, accounting for pure dephasing requires to incorporate at least the quadratic-linear (QL) coupling which gives  (QL) (QL) a|si2 |bc|sj2 |dCij (). (126) ab,cd () = ij

Within the QME, (115), the total dephasing rate between two levels is then given by [38] ab = Rab,ab =

 1 (u) (u ) (QL) (kac + kbc ) − Kaa Kbb Cuu ( = 0). 2 c 

(127)

uu

Here, pure dephasing is accounted for by the last term for which one can use, e.g., the approximate form [116,237,238] (QL)

lim Cij

→0

(II) (II) 4 kB T

() ≈ gi gj

lim (j (QL) () − j (QL) (−)) →0 h¯ (QL) (QL) ≈ gi gj pd ,

(128)

where pd is a (temperature dependent) parameter. Higher order system–bath couplings such as Eq. (102) can be systematically included yielding expressions for the damping matrix in terms of spectral densities [38,116] and thus additional contributions to relaxation and dephasing rates. Of course, this approach will generate even more spectral densities which requires further approximations and lead to additional (usually unknown) parameters. One possibility to obtain at least the intramolecular part of the Redfield tensor has been given in Ref. [186] for chlorotropolone on the basis of an ab initio CRS Hamiltonian. So far we have only discussed those contributions to the Redfield tensor, Eq. (116), which have a straightforward interpretation in terms of energy and phase relaxation rates. This restriction is possible in the Bloch limit where populations and coherences within the density matrix are completely decoupled such that   Rab,cd cd → (1 − ab )Rab,ab ab + ab Raa,cc cc . (129) c

cd

The validity of the Bloch limit can be appreciated by looking at the secular approximation to the equations of motion, (115). Here, in the spirit of a rotating wave approximation, one retains only those terms on the right-hand side which obey |ab − c d  | = 0, i.e.   Rab,cd cd → Rab,c d  c d  . (130) cd

c d 

The conditions for the secular approximation are apparently fulfilled by the Bloch limit terms. However, the right-hand side of Eq. (130) also mixes different coherence matrix elements, i.e. a bath-induced coherence transfer, ab → c d  , is possible [236,239,240]. Beyond the secular approximation the bath-induced conversion of coherences into populations, ab → cc , and vice versa becomes possible as well [241,242]. 6.3. Vibrational energy cascading in an intramolecular hydrogen bond In Section 3.2.1 we have considered the intramolecular HB dynamics of PMME triggered by ultrafast laser-driving. It turned out that without accounting for the influence of the surrounding solvent the experimental observation of rather rapid population relaxation of the OH/OD stretching vibration could not be explained. This issue was reinvestigated using a system–bath model in Refs. [116,236] and supplemented by further experimental data and simulations in Refs. [112,113]. In the following we summarize the conclusions which have been reached on the basis of a QME study (Section 6.2) combined with an ab initio PES including up to four-mode correlations in a normal mode representation of the relevant system (Section 2.1) as well as two-color IR pump–probe experiments. In Fig. 31 we show the IR absorption spectrum of PMME-H dissolved in CCl4 together with the spectra of the different pulses which were chosen such as to probe in the OH -bending fundamental region and to excite in the OH stretching, OH -bending as well as in the C–O stretching region. An additional probe pulse around 2200 cm −1 served

260

K. Giese et al. / Physics Reports 430 (2006) 211 – 276

Fig. 31. Stationary IR absorption spectrum of PMME-H in CCl4 solution. In the lower part of the figure the pulse spectra of the pump and probe pulses are shown. (For more details, see Ref. [112]).

Fig. 32. Two-color IR signals as a function of the delay time of probing the OH = 0 → 1 transition after excitation in the fundamental transition regions of the OH , OH , and C–O stretching vibration, cf. Fig. 31. (For more details, see Ref. [112]).

to monitor the excited state absorption of the OH -stretching band. From the latter the decay time of about 220 fs of the OH = 1 state has been inferred. In Fig. 32 the transient signals probed in the OH = 0 → 1 range are shown for the various excitation conditions. For direct excitation of this transition a fast component (800 fs) is observed which can be attributed to the decay of the bending fundamental state. The same fast component is seen upon excitation of the OH = 0 → 1 stretching fundamental transition, whereas for excitation of the C–O transition only a slow component is found, which is present in all data and due to the overall vibrational cooling process. The initial ultrafast dynamics can be well described within a 5D dissipative quantum model. The five normal modes spanning the PES are depicted in Fig. 33 together with a schematic level scheme showing the fundamental, overtone, and combination transitions of the 4 fast coordinates, i.e. OH , OH , and two modes having out-of-plane bending character, OH1 and OH2 . Each of these states is dressed by an anharmonic potential energy curve corresponding to the motion along the low-frequency HB mode HB . The coupling between these 5 DOF and the remaining intra– and intermolecular DOF has been treated by including into the system–bath Hamiltonian a LL term, Eq. (100), which is responsible for energy dissipation out of HB into the solvent as well as a third order term, Eq. (102), which serves to transfer energy from the modes OH1 and OH2 into intramolecular bath modes. The latter process is assisted by the solvent and therefore could not be observed in the gas phase simulations discussed in Section 3.2. An overall picture of the relaxation process can be obtained from Fig. 34 where we show the population dynamics of the diabatic states defined for the fast coordinates as a function of time and energy. The cascading nature of the dynamics

K. Giese et al. / Physics Reports 430 (2006) 211 – 276

261

Fig. 33. Normal mode displacement vectors of the 5D dissipative model of the PMME-H relaxation dynamics developed in Ref. [112]. In the right part of the figure possible transitions between the four fast mode states, (OH , OH , OH1 , OH2 ), after excitation of the OH = 0 → 1 transition are shown.

is clearly visible and the associated time scales for the OH and OH fundamental transitions are well reproduced. Notice that the modeling includes essentially a single parameter, i.e. the strength of the third-order coupling. The relaxation time for the low-frequency HB mode was determined from a classical simulation to be about 1.7 ps [236]. The importance of the third-order coupling can be appreciated by simply switching it off. Respective results are shown in Fig. 35. In this case the only means for energy release into the environment is the LL coupling via the low-frequency mode. Apparently, in the present model this causes a much slower decay of the initially excited state manifold. Finally, we note that the low-frequency modulations of the signal after OH -stretching excitation in Fig. 32 are due to quasi-coherent wave packet motion of the HB mode (cf. discussion in Section 3.2.1). The nature of this wave packet motion has been investigated, e.g., for a 3D dissipative model of PMME-D in Ref. [236] where it was shown that it is due to coherences in the diabatic ground state potential which are excited by the laser field, but may also be triggered by bath-induced coherence transfer. The example of PMME shows that the combination of time-resolved IR spectroscopy and dissipative quantum theory enables one to obtain a detailed understanding of the dynamics underlying complex stationary IR spectra as given in Fig. 31. It should be noted that similar observations (ultrafast dynamics, coherent wave packet motion) have also been reported for the intermolecular HB in the acetic acid dimer [243–245]. 6.4. HOD in D2 O and H2 O The dynamics of water as the most important solvent has been scrutinized in some detail during the last years [23,246]. In particular the microscopic details of the longstanding problem of unusual proton and hydroxide ion mobility became accessible to computer simulations [247,248]. Being interested in the fundamental building blocks of these systems,

262

K. Giese et al. / Physics Reports 430 (2006) 211 – 276

Fig. 34. Population dynamics of the fast mode diabatic states (shown at their respective energies) after ultrafast (100 fs) excitation of the OH =0 → 1 region. The dissipative model includes a LL coupling to the solvent as well as third-order solvent assisted coupling between the modes OH1 and OH2 and two intramolecular bath modes of comparable frequency. Notice that the population of each state is calculated by summing up all vibrational levels of the low-frequency HB mode for that state. The ground state population is not shown.

Fig. 35. Same as Fig. 34 but without the third-order coupling.

protonated water clusters have triggered quite some attention [249] with the Zundel cation, H5 O+ 2 , emerging as the prototype system for a strong HB whose IR spectrum is still providing a challenge for interpretation [250–252]. These clusters can also be viewed as a means to obtain a deeper insight into the dynamics of proton translocation along “proton wires” in biological water networks such as bacteriorhodopsin, the prototypical light-driven proton pump [193,253] . Returning to bulk water, there is a host of ultrafast spectroscopic investigations which usually focus on HOD in D2 O or H2 O because it allows for a spectral separation between the primarily excited vibrations and the modes of the surrounding. The central questions are concerned with the understanding of the relaxation and dephasing processes which lead, for instance, to the rather broad absorption band. While early time-resolved studies seemed to be limited by their time-resolution [254], more recent IR pump–probe experiments point to a vibrational relaxation time of the excited OH-stretching mode of about 0.7 ps [255,256] to 1.3 ps [257]. Considerable mechanistic information comes

K. Giese et al. / Physics Reports 430 (2006) 211 – 276

263

from IR-pump plus Raman-probe experiments which besides an OH-relaxation time of 1 ps allowed to assign the primary relaxation channel to involve the intramolecular HOD bending and the solvent DOD bending vibrations [258]. First theoretical investigations of vibrational relaxation have been presented by Rey and Hynes using a Landau–Teller type approach which combines a quantum description of the solute with a classical solvent [187]. Their calculated relaxation time of 8 ps is too long and indeed the more recent extensive work by Skinner and coworkers gave an improved value of about 2.3 ps [235]. Both agree, however, on the relaxation pathway starting with the HOD bending overtone. A rather detailed comparison has been given in Ref. [259]. It should be noted that these theoretical approaches assume “intact” intramolecular coordinates for HOD, i.e. the possibility for HB breaking by virtue of energy release into the HB mode is not included. It has, however, been observed in fully classical simulations of OH-frequency gap correlation functions [260] and at least in Ref. [256] it has been argued from the experimental point of view, that the uptake of energy by the HB mode contributes to the OH stretching relaxation. For a more detailed discussion of this partly unresolved issue, see Refs. [23,259]. The OH-vibration of HOD in D2 O has also been investigated using IR photon echo spectroscopy. Two-pulse photon echoes probing the evolution during t1 and t3 while t2 is zero in Eq. (109) gave evidence for a rather fast pure dephasing time of about 90 fs within a Bloch model [261]. This has been assigned within Oxtoby’s model [218] (see Section 6.2) to be a consequence of the large anharmonicity of the OH-mode which allows for dephasing already by a linear coupling model, i.e., Hs.b = sOH f (t), where f (t) is the fluctuating force. Notice that such a term would vanish for a 2 f (t) would lead to pure dephasing (see Section harmonic oscillator coordinate and only the quadratic coupling ∝ sOH 6.1). A more detailed theoretical investigation in Ref. [262] points to the relevance of electric field fluctuations on length scales which originate from the local environment during the first 200 fs and from more collective long-ranged motions at longer times. Spectral diffusion, that is, the randomization of an initially prepared distribution by the environmental fluctuations, has been probed using 3-pulse photon echo spectroscopy, i.e., the echo peak shift is monitored as a function of the delay time t2 . It is well established that this gives direct information on the gap correlation function, Eq. (112). For the present system it was found that the decay of C(t) has picosecond components indicating structural correlations for times well beyond the vibrational relaxation time [263]. Subsequently, even signatures of low-frequency underdamped HB motion have been observed [264] in accord with the simulations of gap-correlation functions in Refs. [260,262,265]. 6.5. HAT and PT reactions In the foregoing two sections we have considered nonreactive vibrational dynamics of HBs. On the other hand, the kinetics and dynamics of reactive HAT and PT in the condensed phase is at the heart of many chemical and biochemical processes (for a general review see, e.g., Ref. [33]). In principle one can distinguish adiabatic and nonadiabatic transfer processes depending on whether or not transitions between the quantum states of the proton are important [38]. Within the QME approach of Section 6.2 there is no need for this distinction since the full information on the quantum states and respective transitions is part of the model and naturally accounted for by solving the equations of motion [266–271]. In the context of the hybrid quantum–classical simulations discussed in Section 5, however, the distinction between adiabatic and nonadiabatic is essential. Adiabatic PT including the quantum character of the proton wave function has been reported, for instance, for model HB complexes in polar solution [272], excess proton motion in water [273], and enzyme-catalyzed hydrolysis [202]. In Ref. [202] it was also shown that the quantum-classical approach is a limiting case of a TDSCF ansatz for the total wave function, i.e. it is essentially a mean-field approach. Nonadiabatic transitions between the adiabatic PES can be accounted for along a classical trajectory of the heavy nuclei by means of the surface hopping procedure [274,275]. Combined with a multiconfiguration description of the proton wave function this approach has been termed MC-MDQT (multiconfigurational molecular dynamics with quantum transitions) [203]. MC-MDQT has been applied to various processes ranging from excess PT in water wires [276] to hydride transfer in enzymes [277–279]. 6.6. Multidimensional IR spectroscopy Inspired by multidimensional NMR spectroscopy, considerable efforts have been directed towards carrying some of the respective ideas into IR spectroscopy. The advantages are obvious: additional information is obtained on congested

264

K. Giese et al. / Physics Reports 430 (2006) 211 – 276

Fig. 36. Homodyne-detected signal, Eq. (131), of a dissipative double minimum system with LL and QL coupling. Left and right column correspond to the Bloch limit and the full Redfield tensor simulation, respectively. In the upper panels t1 = 0 whereas for the lower panels we have t1 = 1000 fs. For more details see Ref. [221].

spectra with the benefit—as compared to NMR spectroscopy—of having a high time-resolution [216,280–282]. Within the context of HBs the fluctuations, e.g. of the OH/OD-stretching transition in HOD diluted in D2 O/H2 O have been a target for 2D-IR studies [283,284]. But also processes in other systems such as the HB breaking in methanol-OD oligomers [285,286] have been investigated. Rather promising are the prospects for unravelling HAT or PT reaction dynamics as shown in theoretical calculations for malonaldehyde [81] and generic dissipative model systems [221,287]. Tanimura and coworker, for instance, have demonstrated a novel technique for obtaining chemical reaction rates directly from 2D-IR spectra [287]. 2D-IR spectroscopy is rather sensitive to anharmonic couplings of the PES. This feature has been used in Ref. [221] to investigate the dependence of the dissipative dynamics of a reaction coordinate on the type of system–bath coupling. In particular it was demonstrated that the cross-terms appearing for simultaneous LL and QL coupling, i.e. Rabcd ∝ a|si |bc|si2 |d, can lead to a breakdown of the Bloch approximation to the QME. Specifically, this term introduces an efficient bath-induced coherence transfer where the frequency mismatch in the Redfield tensor is of the order of the ground state tunneling splitting. In other words, although the latter can be small and the related dynamics well in the picosecond range, there is an influence on the ultrafast dynamics due to coherence transfer. As an example we show in Fig. 36 the Fourier-transformed homodyne-detected signal Shom ( 3 , t2 , t1 ) =

∞ 0

2 dt3 ei 3 t3 R (3) (t3 , t2 , t1 )

(131)

for different coherence times t1 . The model system has two pairs of states below the barrier, (|0, |1) and (|2, |3). In −1 the figure one sees peaks at 3 = 43 = 1048, 21 = 1630, √ and 30 = 1936 cm . Initially the system has been prepared in a nonequilibrium superposition state (|2 + |3)/ 2. Panels (a) and (b) correspond to the Bloch limit for t1 = 0 and 1000 fs, respectively. The decay of the different peaks is in the subpicosecond range due to population relaxation (t2 -axis). In addition, coherence dephasing during t1 leads to an overall reduction of the amplitude. A similar behavior is seen when the full QME is solved, Fig. 36c and d. However, there is an important difference for the resonance at 3 = 21 , for which the dephasing and relaxation times are apparently longer as compared to the Bloch case. In Ref. [221] it was shown that this is a consequence of coherence transfer of the type 21 → 20 → 21 .

K. Giese et al. / Physics Reports 430 (2006) 211 – 276

265

7. Sundry topics 7.1. Geometric isotope effects The anharmonicity of the PES together with the wave nature of hydrogen motion leads to a measurable quantum effect already for the vibrational ground state. Consider the case of an asymmetric double minimum potential. The expectation value of the bond length RA.H quant for an A.H · · · B bond will be larger than the respective classical value, RA.H . Upon deuteration the lower ZPE leads to a reduction of the deviation between quantum and classical results such that RA.H quant > RA.D quant > RA.H . This is called primary geometric isotope effect. For symmetric single minimum potentials of strong HBs, A · · · H · · · A, on the other hand, the expectation value of the H or D position is always midway between the donor and the acceptor. So far we have not yet considered the coupling between the A–H motion and the hydrogen bond A · · · B motion. In fact, this may cause a simultaneous change of the expectation value of the A · · · B distance upon deuteration. This is called secondary geometric isotope effect. It was found in the 1930s in molecular crystals by Ubbelohde (Ubbelohde effect) [288] (for an early review, see also Ref. [289]). Generally speaking, for the asymmetric A.H · · · B case the reduction of the A–H distance upon deuteration will cause a weakening of the HB, i.e. an increase of the A · · · B distance. The converse will occur for the single minimum A · · · H · · · A case. Upon increasing the A · · · A distance, e.g. for different crystals, the magnitude of the secondary geometric isotope effect will change. Ichikawa compiled a correlation diagram based on available crystal structures which showed that a HB expansion occurs for A · · · A distances ˚ The maximum expansion is about 0.03 A ˚ [290,291]. In fact, much of this behavior has in between 2.43 and 2.65 A. been rationalized by using a simple empirical chain model including the effect of quantum fluctuation of the H/D atoms [292]. In passing we note that solid state [293] as well as low-temperature liquid state [294] NMR offer an alternative to X-ray or neutron scattering to study geometric isotope effects. An interesting case occurs when the ZPEs are just below or above the barrier such that there is a substantial difference between the H and D wave functions. While the expectation value of the H/D position will be the same, the wave function will have a single maximum in the former case while it has two local maxima in the latter case. Consequently, the HB will contract in the former case and expand in the latter one. Recently, this effect has been studied in the context of an empirical valence bond model by Limbach et al. [295]. The limiting case where the ZPE for H is just above and that for D is just below the barrier seems to be realized in the H3 O− 2 ion which has been studied using ab initio path integral molecular dynamics [296]. The calculation of the geometric isotope effect requires knowledge of the respective vibrational ground state wave function. This becomes particularly straightforward in cases where this function is rather localized either in a single minimum potential [299] or in a double minimum potential where the tunneling coupling can be neglected [297,298]. Here one can use a PES spanned by the local normal mode coordinates as given in Eq. (10). In Fig. 37 we show an example of a four-dimensional model for the double HB in porphycene. From the two-dimensional cut which

Fig. 37. Geometric isotope effect in porphycene: Two normal mode displacement vectors for which the PES cut in the right panel is shown together with the ground state vibrational probability density (for more details, see Refs. [297,298]).

266

K. Giese et al. / Physics Reports 430 (2006) 211 – 276

covers the symmetric NH-stretching as well as a low-frequency collective ring deformation mode one notices that the maximum of the vibrational probability density is not at the coordinate origin (classical limit). Since the doubly deuterated species has a lower ZPE the probability density will cover a less anharmonic part of the PES and its maximum will be closer to the classical one. Thus, the simultaneous change of NH-distances and HB geometry mediated by the anharmonic coupling leads to the primary and secondary geometric isotope effects. Finally, we note that multicomponent quantum chemistry [39] offers in principle a rather well-suited alternative to the construction of a Born–Oppenheimer PES. Here one has the advantage that a geometry optimization of the classical coordinates yields automatically information on the geometric isotope effects. However, this does not yet include the effect of quantum mechanical correlations between the proton and the heavy atom motion, which needs to be taken into account in a subsequent heavy atom wave function calculation on the electron–proton PES. Unfortunately, for a molecule of the size of porphycene, for instance, this approach is only feasible upon making further approximations concerning the quantum chemical model [297]. 7.2. Excited state PT Inter- and intramolecular PT after electronic excitation has attracted considerable attention, e.g., in the context of industrial photostabilizers [21,300] and biological photoprotection [301]. In contrast to the IR driven dynamics in the ground state, the forces experienced by the ground state wave packet upon promotion to the excited state are considerable, leading to rather rapid PT reactions. It is also important to note that the PT will be triggered by the laser interaction at a specific instant of time, opening the way to a better understanding of such reactions in various environments. For instance, using femtosecond pump–probe spectroscopy it was found that the photostabilizer Tinuvin completes a reaction cycle—excitation, PT, internal conversion, PT—in about 1 ps, where the excited state PT itself takes about 150 fs [21]. This suggest an essentially barrierless process which is, however, part of a multidimensional wave packet motion involving several low-frequency DOFs. Oscillatory spectra signaling coherently excited nuclear wave packets where first found in Ref. [302] (see, also Ref. [303]). An even better resolution of the coherent wave packet motion was achieved in Ref. [304] where a benzothiazole compound has been investigated which has a much longer S1 lifetime (for a review see also Ref. [305]). The excitation of four vibrational modes could be identified and the 60 fs excited state PT was termed to occur ballistically, i.e. without much wave packet dispersion. A theoretical study in terms of normal modes projected on the minimum energy path was given in Ref. [306]. Condensed phase excited state intermolecular PT is of particular importance in the context of photoacid chemistry. Upon electronic excitation photoacids, such as certain aromatic dye molecules increase their acidity substantially and in the presence of a base PT occurs. From the kinetics point of view, this problem has traditionally been handled as a bimolecular reaction whose rate is governed by diffusional motion and PT in the encounter complex (for a review see Ref. [307]). Recently, it was demonstrated that the anharmonic couplings between the reaction coordinate for PT and the vibrations of the proton donor and acceptor can be used to follow the reaction on its intrinsic time scale. To this end, the optical excitation pulse is combined with a mid-IR pulse which observes the transient changes of specific marker modes (for a review, see Ref. [308]). This strategy had been introduced into excited state PT research in Ref. [309] about 20 years ago. It allowed to unravel, for instance, the ultrafast PT steps in the complex between the photoacid 8-hydroxy-1,3,6-trisulfonate-pyrene and acetate in deuterated water [310], including the possibility of proton hopping along water wires in so-called loose complexes [311]. Although much theoretical work has been focussed on aspects of PT in acid–base chemistry [25], a microscopic modeling of the transient spectra has not yet been accomplished. There are also a number of double PT systems which have been studied, putting emphasis on the issue of stepwise versus concerted transfer. For instance, for BP(OH)2 ([2, 2 -bipyridyl]-3, 3 -diol) it was found that step-wise and concerted transfer occur at the same time, both with a 100 fs time scale [312]. Here, the concerted pathway showed a pronounced dependence on the excitation wave length. Wave packet motion of skeleton modes associated with the two reaction pathways have been investigated in Ref. [305]. Excited state double PT has also been investigated for the 7-azaindole dimer as a model mimicking DNA base pairs. The issue of the transfer mechanism appeared to be controversial for a long time [313,314]. The most recent time-resolved resonance-enhanced multiphoton ionization (REMPI) study of the gas phase dimer gives evidence for a concerted mechanism [315]. Other model base pairs, like the 2-aminopyridine dimer, have been studied with regard to the photostability. Here, conical intersections reached along the PT coordinate are testifying the breakdown of the Born–Oppenheimer approximation and are the key for efficient deactivation after electronic excitation [316].

K. Giese et al. / Physics Reports 430 (2006) 211 – 276

267

7.3. Laser control of hydrogen transfer Coherent laser control is emerging as a powerful tool to influence and study nonreactive and reactive molecular dynamics [317–320]. An aspect crucial for the realization of control is the coherence of the evolving wave packet. With the demonstration of coherent wave packet dynamics of HBs [108] in solution, the question can be raised whether electronic ground state PT or HAT can be triggered by adapted IR pulses [321]. The case where the dynamics is restricted to the two lowest vibrational states in a double well potential has been extensively investigated (for a review, see Ref. [322]). Taking into account higher excited states, the pump–dump scheme [323] is perhaps the most straightforward control strategy. For an asymmetric PES the initial state which may be localized in the left minimum is first excited to a state which is energetically above the barrier (pump). A second pulse then de-excited (dumps) the system to a state localized in the right minimum. This strategy has been applied to a generic two-dimensional model in Ref. [324] and to a two-dimensional model of thioacetylacetone in Ref. [269]. While parameterized pulse sequences have the advantage that the optimized pulse has a straightforward interpretation, optimal control theory which optimizes the laser field such as to achieve, e.g., a maximum overlap between the wave packet and some target state at a given time turned out to be much more powerful in predicting, e.g., new reaction pathways. For instance, in Ref. [269] it was demonstrated that upon increasing the penalty for high pulse intensities during the optimization, the HAT dynamics changes from an above-the-barrier pump–dump like scheme to that of a through-the-barrier tunneling. The latter mechanism was termed “Hydrogen Subway” in Ref. [325]. This mechanism which does not involve highly excited vibrational states was later shown to particularly useful in condensed phase situation where phase and energy relaxation competes with coherent laser driving [270,271,326,327]. (For a general account on the influence of pure dephasing on HAT control, see also Ref. [328].) Given the success of low-dimensional model studies one might ask whether the concepts can be applied to multidimensional systems. In this case making a guess for an analytical pulse sequence would be a real challenge and much hope lies on optimal control theory. However, the implementation of the obvious combination of optimal control theory and MCTDH has not yet been reported. In Fig. 38 we show an example indicating the complications one faces in multidimensional HAT control [122]. For illustration we consider the 7D model of deuterated salicylaldimine discussed in Section 3.2.2 [121]. The two dimensions corresponding to the planar HAT are used to define diabatic states (cf. discussion in Section 3.2.2). Each diabatic state is dressed by a 5D harmonic PES describing the modes shown in Fig. 13. The coordinate dependent diabatic state couplings are included but their effect on the PES is not shown in the figure. The initial state is the ground state wave packet in the potential curve corresponding to the diabatic state localized in the keto well (cf. Fig. 14 d) (thick line). For the generation of the initial wave packet the state couplings are switched-off. This could be a choice of a product target state in a control scheme. Due to the anharmonic couplings

Fig. 38. Left: diabatic (with respect to the D-motion) potential energy curves along mode 14 for a 7D model of deuterated salicylaldimine (cf. Fig. 13). The thick solid curves represent the ground state ( = 1) and the initial state ( = 10) which is populated at time zero in the dynamics simulation shown in the right panel. Right: Population dynamics of the ten lowest diabatic states after initial population of the HAT product like state  = 10 (curves are shifted by ( − 1) × 0.25). For details, see also Ref. [122].

268

K. Giese et al. / Physics Reports 430 (2006) 211 – 276

between the harmonic modes and between the diabatic states the decay of the initial state back to the overall ground state is rather rapid as seen from the diabatic state populations in Fig. 38. However, there appears to remain some small population in the initial state during the first picosecond. Moreover, a detailed analysis shows that the energy is almost equally redistributed over all five harmonic modes. In other words, at least for this system, the dynamics is truly multidimensional and the choice of a reaction path or even of a sufficiently stable target state for laser driving is not obvious.

8. Summary The theoretical description of hydrogen bond dynamics continues to provide a challenge since it comprises multidimensionality and pronounced quantum effects. This setting would be interesting on its own because it can be argued that methods developed in its context are rather likely to be transferable to other systems with less quantum behavior. But, the intriguing aspect of hydrogen bond research is its broad scope in terms of possible application ranging from the properties of water via biological light-triggered proton pumps to enzyme catalysis. The last years have witnessed an enormous increase of the microscopic insight into these processes due to experimental and theoretical progress. In this review we have focussed on recent theoretical efforts to understand the stationary and time-resolved infrared spectroscopy of hydrogen-bonded systems starting from a microscopic multidimensional Hamiltonian. Here, the strategies range from anharmonic normal mode potentials for nonreactive vibrational dynamics up to reaction surfaces for reactive processes. Considerable space has been devoted to the presentation of the CRS Hamiltonian approach, because it is not only a reasonable compromise between accuracy and the desire for a truly multidimensional ab initio potential, but it also provides the interface to commonly used generic system–bath Hamiltonians for condensed phase reactions. The availability of a multidimensional potential is only useful if that the respective dynamical equations can be solved. Here, the CRS Hamiltonian turns out to be ideally suited for the combination with the MCTDH method for propagating wave packets. This statement has been underlined by three applications in Section 3. First, we have shown how the anharmonic coupling between the OH stretching vibration and some low-frequency hydrogen-bond mode in a medium strong single minimum system can lead to coherent wave packet motions of the low-frequency mode after laser excitation of the OH-mode. While this example did not show any ultrafast intramolecular vibrational energy redistribution (IVR), increased anharmonicity as present in a double minimum system was shown to lead to subpicosecond energy redistribution even in this gas phase situation. In the final example we examined the issue of mode-specific tunneling in a symmetric double minimum hydrogen bond, which provides just a specific realization of the general theme of the influence of promoting and reorganization motions in proton or hydrogen atom transfer reactions. Stepping from the gas to the condensed phase the question can be raised whether classical molecular dynamics with its favorable scaling can be supplemented by quantum effects at moderate costs. We have looked at this issue in terms of the ability of trajectory-based methods to account for quantum tunneling. Thereby, we outlined an extension to the well-known Makri–Miller method which is based on trajectories in the classically forbidden region. Although we showed that this new method is capable of improving the description, for instance, of tunneling splittings, it must be stated that this comes along with substantial numerical effort. We have also discussed the more traditional approach to condensed phase dynamics in terms of a quantum master equation. This approach has been illustrated for the description of energy cascading in an intramolecular hydrogen bond, where its combination with a nontrivial system–bath Hamiltonian provided a microscopic picture for the time-scales observed in ultrafast pump–probe spectroscopy. We further presented some results for the two-dimensional infrared spectroscopy of a generic proton transfer system which showed how this rather promising technique can give detailed information on the anharmonic couplings. Nonlinear time-resolved infrared spectroscopy is a valuable means to unravel even complex dynamics which is otherwise buried in broad stationary spectra. On the other hand, the availability of the appropriate laser source will trigger the desire to not only analyze hydrogen bond dynamics, but to control it by virtue of specifically shaped pulses. Some results along these lines have been discussed at the end of this review together with the challenges posed in this respect by the multidimensional nature of hydrogen bond dynamics.

K. Giese et al. / Physics Reports 430 (2006) 211 – 276

269

Acknowledgment This work has been financially supported by the Deutsche Forschungsgemeinschaft, partly through the Sfb 450. The authors thank N. Došli´c, J. Dreyer, T. Elsaesser, K. Heyne, H.-H. Limbach, J. Manz, H.-D. Meyer, E. T. J. Nibbering, M. F. Shibl, J. Stenger, K. Takatsuka, Y. Tanimura, and H. Ushiyama for many stimulating discussions during the last years. Appendix A. Determination of the SMC Hamiltonian parameters In Ref. [60] a simple method was proposed to determine the SMC parameters entering Eq. (18). It makes use of the reaction plane idea which is introduced in Section 2.3.1. The saddle point of the SMC-PES is at (0, ) where  = /2 is the ratio of the geometry displacements, ≡

|XC − XTS | . |XR − XC |

(A.1)

˜ 1 . The frequencies ˜ 1 and ˜ 2 correspond to the reaction plane vectors ˜ 2 / The frequency parameter  is given by  = w1 and w2 , respectively, and are determined by weighting the contribution of each normal mode with a component into the w1 and w2 direction, respectively,  ˜ 2k =

(Yj · wk )2 2j , k = 1, 2, (A.2) j

where Yj with j = 1, . . . , N − 6 are the N = 3 · Nat dimensional normal mode vectors of the (right or left) minimum ˜ 1 /8E , where E and j are the corresponding normal mode frequencies. The SMC parameter g is given by g = is the ZPE corrected barrier height. This procedure has been applied in Ref. [60] to tropolone and in Section 4.2.3 to various carboxylic acid dimers. Appendix B. Calculation of multidimensional eigenstates In the following we briefly outline how low lying eigenstates of a multidimensional double minimum PES can be determined approximately by a combination of improved MCTDH relaxation together with solution of a generalized (0,R) (0,L) eigenvalue problem. For this purpose the states n (n ) corresponding to the harmonic approximation to the right (left) minimum serve as initial states for MCTDH improved relaxations (cf. Section 3.1) which gives the relaxed (rel,R) (rel,L) and n . These states form an optimized basis for the diagonalization of the Hamiltonian. The left states n and right states are not orthogonal, thus eigenstates and eigenvectors of the Hamiltonian can be obtained by solving the generalized eigenvalue problem,    ,  ,  (Hn,n (B.1)  − εl Sn,n ) · cn ,l = 0, n  =R,L



(rel,)

 , where  = R or L, Hn,n  = n 

(rel, )

|H |n



(rel,)

,  is the Hamiltonian matrix, Sn,n  = n

matrix, and cn ,l are the expansion coefficients for the lth eigenstate,    ) l = cn,l (rel, , n

(rel, )

|n

 is the overlap

(B.2)

n =R,L

with eigenenergy εl . The present procedure requires less single particle functions than the attempt to compute the gerade and ungerade states directly by improved relaxation. Appropriate symmetry adapted gerade and ungerade superpositions can be formed according to (rel,±) ∝ (rel,R) ± (rel,L) . n n n An application is given in Section 3.3 (cf. Table 2).

(B.3)

270

K. Giese et al. / Physics Reports 430 (2006) 211 – 276

Appendix C. Parities of motion The assignment of parities of motion, j = ±1, to each DOF has been suggested by Takatsuka an coworkers [178]. It insures that real-valued momenta can always be defined as √ p¯ j = pj / j . (C.1) For j = +1 and −1 the motion is classically allowed and forbidden, respectively. Of course, the transformation Eq. (C.1) is noncanonical if j = −1. Inserted into the usual Hamiltonian it yields  j ¯ q; ) = p¯ 2 + V (q), (C.2) H¯ (p, 2 j j

¯ q) space can be generated by modified where  = (1 , . . . , N ) is the vector of parities of motion. Trajectories in the (p, Hamilton’s equation of motion, jH¯ jV p˙¯ j = − =− , jqj jqj jH¯ q˙j = = j p¯ j . jp¯ j

(C.3) (C.4)

The first equation, Eq. (C.3), is the unchanged Newton’s equation of motion. The second equation, Eq. (C.4), determines that for j =−1 the velocity q˙j and momentum p¯ j have opposite directions. Notice that Eqs. (C.3)–(C.4) are equivalent to the formulation in Ref. [178] and have a canonically invariant form. This implies that in principle the method of characteristics can be used to construct an action function based on a field of trajectories [329]. However, for the present purpose we defined the action along a trajectory on a -sheet (for constant energy) by  t √ S(p¯ 0 , q0 , ; t) = j p¯ j ()q˙j () d, (C.5) j

0

where (p¯ 0 , q0 ) are the initial conditions of the trajectory. References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18] [19] [20] [21] [22]

P. Schuster, G. Zundel, C. Sandorfi (Eds.), The Hydrogen Bond Theory, North-Holland, Amsterdam, 1976. T. Bountis (Ed.), Proton Transfer in Hydrogen-Bonded Systems, Plenum Press, New York, 1992. D.A. Smith (Ed.), Modeling the Hydrogen Bond, ACS Symposium Series, vol. 569, ACS, Washington, 1994. G.A. Jeffrey, An Introduction to Hydrogen Bonding, Oxford, New York, 1997. S. Scheiner, Hydrogen Bonding, Oxford, New York, 1997. D. Hadži (Ed.), Theoretical Treatments of Hydrogen Bonding, Wiley, Chichester, 1997. H.-H. Limbach, J. Manz, Ber. Bunsenges. Phys. Chem. 102 (1998) 289. G.R. Desiraju, T. Steiner, The Weak Hydrogen Bond, Oxford University Press, Oxford, 1999. P. Schuster, P. Wolschann, Chemical Monthly 130 (1999) 947. T. Elsaesser, H.J. Bakker (Eds.), Ultrafast Hydrogen Bonding Dynamics and Proton Transfer Processes in the Condensed Phase, Kluwer Academic Publishers, Dordrecht, 2002. A. Kohen, H.-H. Limbach (Eds.), Isotope Effects in Chemistry and Biology, Taylor and Francis, Boca Raton, 2005. R.L. Schowen (Ed.), Handbook of Hydrogen Transfer, Wiley-VCH, Weinheim, 2006. W.M. Latimer, W.H. Rodebush, J. Am. Chem. Soc. 42 (1920) 1419. M.L. Huggins, Angew. Chem. Int. Ed. 10 (1971) 147. J.D. Bernal, R.H. Fowler, J. Chem. Phys. 1 (1933) 515. M.L. Huggins, J. Phys. Chem. 40 (1936) 723. L. Pauling, The Nature of the Chemical Bond, Cornell University Press, Ithaka, 1939. R.M. Badger, S.H. Bauer, J. Chem. Phys. 5 (1936) 839. L. Pauling, R.B. Corey, H.R. Branson, Proc. Nat. Acad. Sci. USA 37 (1951) 205. J.D. Watson, F.H.C. Crick, Nature 171 (1953) 737. T. Elsaesser, in: J. Manz, L. Wöste (Eds.), Femtosecond Chemistry, vol. 2, Verlag Chemie, Weinheim, 1994, p. 563. A. Douhal, F. Lahmani, A.H. Zewail, Chem. Phys. 207 (1996) 477.

K. Giese et al. / Physics Reports 430 (2006) 211 – 276

271

[23] E.T.J. Nibbering, T. Elsaesser, Chem. Rev. 104 (2004) 1887. [24] M. Benoit, D. Marx, Chem. Phys. Chem. 6 (2005) 1738. [25] P.M. Kiefer, J.T. Hynes, in: T. Elsaesser, H.J. Bakker (Eds.), Ultrafast Hydrogen Bonding Dynamics and Proton Transfer Processes in the Condensed Phase, Kluwer Academic Publishers, Dordrecht, 2002, p. 73. [26] S. Bratos, J.-C. Leicknam, G. Gallot, H. Ratajczak, in: T. Elsaesser, H.J. Bakker (Eds.), Ultrafast Hydrogen Bond Dynamics and Proton Transfer Processes in the Condensed Phase, Kluwer Academic, Dordrecht, 2002, p. 5. [27] R.P. Bell, The Tunnel Effect in Chemistry, Chapman & Hall, London, 1980. [28] T. Baba, T. Tanaka, I. Morino, K.M.T. Yamada, K. Tanaka, J. Chem. Phys. 110 (1999) 4131. [29] F. Madeja, M. Havenith, J. Chem. Phys. 117 (2002) 7162. [30] H. Sekiya, Y. Nagashima, Y. Nishimura, J. Chem. Phys. 92 (1990) 5761. [31] S. Hammes-Schiffer, Chem. Phys. Chem. 3 (2002) 33. [32] L. Rodriguez-Santiago, M. Sodupe, A. Oliva, J. Bertran, J. Am. Chem. Soc. 121 (1999) 8882. [33] M.V. Basilevsky, M.V. Vener, Russ. Chem. Rev. 72 (2003) 3. [34] Y. Tomioka, M. Ito, N. Mikami, J. Phys. Chem. 87 (1983) 4401. [35] C. Manca, C. Tanner, S. Coussan, A. Bach, S. Leutwyler, J. Chem. Phys. 121 (2004) 2578. [36] P.K. Agarwal, S.R. Billeter, P.T.R. Rajagopalan, S. Benkovic, S. Hammes-Schiffer, Proc. Nat. Acad. Sci. USA 99 (2002) 2794. [37] R.S. Sikorski, L. Wang, K.A. Markham, P.T.R. Rajagopalan, S.J. Benkovic, A. Kohen, J. Am. Chem. Soc. 126 (2004) 4778. [38] V. May, O. Kühn, Charge and Energy Transfer Dynamics in Molecular Systems, 2nd ed., Wiley-VCH, Weinheim, 2004. [39] M.V. Pak, C. Swalina, S.P. Webb, S. Hammes-Schiffer, Chem. Phys. 304 (2004) 227. [40] B.I. Stepanov, Nature 157 (1946) 808. [41] N. Sheppard, in: D. Hadži (Ed.), Hydrogen Bonding, Pergamon Press, New York, 1957, p. 85. [42] S. Bratos, D. Hadži, J. Chem. Phys. 27 (1957) 991. [43] A. Witkowski, J. Chem. Phys. 47 (1967) 3645. [44] Y. Marechal, A. Witkowski, J. Chem. Phys. 48 (1968) 3697. [45] O. Henri-Rousseau, P. Blaise, Adv. Chem. Phys. 103 (1998) 1. [46] O. Henri-Rousseau, D. Chamma, Chem. Phys. 229 (1998) 37. [47] D. Chamma, O. Henri-Rousseau, Chem. Phys. 229 (1998) 51. [48] O. Henri-Rousseau, P. Blaise, D. Chamma, Adv. Chem. Phys. 121 (2002) 241. [49] J. Dreyer, J. Chem. Phys. 122 (2005) 184306. [50] J.K.G. Watson, Mol. Phys. 15 (1968) 479. [51] S. Carter, J.M. Bowman, J. Chem. Phys. 108 (1998) 4397. [52] N. Došli´c, O. Kühn, Z. Phys. Chem. 217 (2003) 1507. [53] S. Carter, S.J. Culik, J.M. Bowman, J. Chem. Phys. 107 (1997) 10458. [54] X. Huang, B.J. Braams, S. Carter, J.M. Bowman, J. Am. Chem. Soc. 126 (2004) 5042. [55] J.M. Bowman, X. Zhang, A. Brown, J. Chem. Phys. 119 (2003) 646. [56] S. Takada, H. Nakamura, J. Chem. Phys. 100 (1994) 98. [57] V.A. Benderskii, S.Y. Grebenshchikov, G.V. Mil’nikov, E.V. Vetoshkin, Chem. Phys. 188 (1994) 19. [58] V.A. Benderskii, S.Y. Grebenshchikov, G.V. Mil’nikov, Chem. Phys. 198 (1995) 281. [59] V.A. Benderskii, S.Y. Grebenshchikov, G.V. Mil’nikov, Chem. Phys. 194 (1995) 1. [60] S. Takada, H. Nakamura, J. Chem. Phys. 102 (1995) 3977. [61] M.J. Wojcik, H. Nakamura, S. Iwata, W. Tatara, J. Chem. Phys. 112 (2000) 6322. [62] E.B. Wilson, J.C. Decius, P.C. Cross, Molecular Vibrations, Dover, New York, 1955. [63] D. Heidrich (Ed.), The Reaction Path in Chemistry: Current Approaches and Perspectives, Kluwer Academic, Dordrecht, 1995. [64] W.H. Miller, N.C. Handy, J.E. Adams, J. Chem. Phys. 72 (1980) 99. [65] Y.-P. Liu, G.C. Lynch, T.N. Truong, D.-H. Lu, D.G. Truhlar, B.C. Garrett, J. Am. Chem. Soc. 115 (1993) 2408. [66] W.H. Miller, B.A. Ruf, Y.-T. Chang, J. Chem. Phys. 89 (1988) 6298. [67] B. Fehrensen, D. Luckhaus, M. Quack, Z. Phys. Chem. 209 (1999) 1. [68] S. Carter, N.C. Handy, J. Chem. Phys. 113 (2000) 987. [69] D.P. Tew, N.C. Handy, S. Carter, Mol. Phys. 99 (2001) 393. [70] R. Meyer, T.-K. Ha, Mol. Phys. 101 (2003) 3263. [71] D.P. Tew, N.C. Handy, S. Carter, Phys. Chem. Chem. Phys. 3 (2001) 1958. [72] S. Schweiger, B. Hartke, G. Rauhut, Phys. Chem. Chem. Phys. 6 (2004) 3341. [73] J. Manz, J. Römelt, Chem. Phys. Lett. 81 (1981) 179. [74] T. Carrington, W.H. Miller, J. Chem. Phys. 81 (1984) 3942. [75] T. Carrington, W.H. Miller, J. Chem. Phys. 84 (1986) 4364. [76] N. Shida, P.F. Barbara, J.E. Almlöf, J. Chem. Phys. 91 (1989) 4061. [77] N. Shida, P.F. Barbara, J.E. Almlöf, J. Chem. Phys. 94 (1991) 3633. [78] I. Matanovi´c, N. Došli´c, J. Phys. Chem. A 109 (2005) 4185. [79] B.A. Ruf, W.H. Miller, J. Chem. Soc. Faraday Trans. 2 84 (1988) 1523. [80] K. Yagi, T. Taketsugu, K. Hirao, J. Chem. Phys. 115 (2001) 10647. [81] T. Hayashi, S. Mukamel, J. Phys. Chem. A 107 (2003) 9113. [82] B.C. Garrett, D.G. Truhlar, J. Chem. Phys. 79 (1983) 4931.

272 [83] [84] [85] [86] [87] [88] [89] [90] [91] [92] [93] [94] [95] [96] [97] [98] [99] [100] [101] [102] [103] [104] [105] [106] [107] [108] [109] [110] [111] [112] [113] [114] [115] [116] [117] [118] [119] [120] [121] [122] [123] [124] [125] [126] [127] [128] [129] [130] [131] [132] [133] [134] [135] [136] [137] [138] [139] [140] [141]

K. Giese et al. / Physics Reports 430 (2006) 211 – 276 C.S. Tautermann, A.F. Voegele, T. Loerting, K.R. Liedl, J. Chem. Phys. 117 (2002) 1962. C.S. Tautermann, A.F. Voegele, T. Loerting, K.R. Liedl, J. Chem. Phys. 117 (2002) 1967. C.S. Tautermann, A.F. Voegele, K.R. Liedl, J. Chem. Phys. 120 (2004) 631. K. Giese, O. Kühn, J. Chem. Phys. 123 (2005) 054315. A. Nauts, X. Chapusiat, Chem. Phys. Lett. 136 (1987) 164. F. Gatti, Y. Justum, M. Menou, A. Nauts, X. Chapusiat, J. Mol. Spectr. 181 (1997) 403. D. Luckhaus, J. Phys. Chem. A 110 (2005) 3151. P.A.M. Dirac, Proc. Cambridge Philos. Soc. 26 (1930) 376. R.B. Gerber, V. Buch, M.A. Ratner, J. Chem. Phys. 77 (1982) 3022. N. Makri, W.H. Miller, J. Chem. Phys. 87 (1987) 5781. O. Vendrell, H.-D. Meyer, J. Chem. Phys. 122 (2005) 104505. A.D. Hammerich, R. Kosloff, M.A. Ratner, Chem. Phys. Lett. 171 (1990) 97. H.-D. Meyer, U. Manthe, L.S. Cederbaum, Chem. Phys. Lett. 165 (1990) 73. M.H. Beck, A. Jäckle, G.A. Worth, H.-D. Meyer, Phys. Rep. 324 (2000) 1. H. Meyer, G.A. Worth, Theor. Chem. Acc. 109 (2003) 251. J.C. Light, I.P. Hamilton, J. Chem. Phys. 82 (1985) 1400. D.T. Colbert, W.H. Miller, J. Chem. Phys. 96 (1992) 1982. G. Worth, M. Beck, A. Jäckle, H.-D. Meyer, The MCTDH package, Version 8.3, University of Heidelberg, Heidelberg, 2005. G.A. Worth, H.-D. Meyer, L.S. Cederbaum, J. Chem. Phys. 109 (1998) 3518. N. Makri, Chem. Phys. Lett. 169 (1990) 541. G.A. Worth, H.-D. Meyer, L.S. Cederbaum, J. Chem. Phys. 105 (1996) 4412. H. Wang, M. Thoss, J. Chem. Phys. 119 (2003) 1289. G.A. Worth, I. Burghardt, H.-D. Meyer, Chem. Phys. Lett. 368 (2003) 502. A. Jäckle, H.-D. Meyer, J. Chem. Phys. 104 (1996) 7974. J. Stenger, D. Madsen, J. Dreyer, E.T.J. Nibbering, P. Hamm, T. Elsaesser, in: T. Elsaesser, S. Mukamel, M. Murnane, N. Scherer (Eds.), Ultrafast Phenomena XII, Springer Series in Chemical Physics, Springer, New York, 2000, p. 542. J. Stenger, D. Madsen, J. Dreyer, E.T.J. Nibbering, P. Hamm, T. Elsaesser, J. Phys. Chem. A 105 (2001) 2929. D. Madsen, J. Stenger, J. Dreyer, P. Hamm, E.T.J. Nibbering, T. Elsaesser, Bull. Chem. Soc. Japan 75 (2002) 909. J. Stenger, D. Madsen, J. Dreyer, P. Hamm, E.T.J. Nibbering, T. Elsaesser, Chem. Phys. Lett. 354 (2002) 256. D. Madsen, J. Stenger, J. Dreyer, E.T.J. Nibbering, P. Hamm, T. Elsaesser, Chem. Phys. Lett. 341 (2001) 56. K. Heyne, E.T.J. Nibbering, T. Elsaesser, M. Petkovi´c, O. Kühn, J. Phys. Chem. A 108 (2004) 6083. K. Heyne, E.T.J. Nibbering, T. Elsaesser, M. Petkovi´c, O. Kühn, in: T. Kobayashi, T. Okada, T. Kobayashi, K.N. Nelson, S. DeSilvestri (Eds.), Ultrafast Phenomena XIV, Springer, 2004, p. 389. G.K. Paramonov, H. Naundorf, O. Kühn, Eur. J. Phys. D 14 (2001) 205. H. Naundorf, G.A. Worth, H.-D. Meyer, O. Kühn, J. Phys. Chem. A 106 (2002) 719. O. Kühn, J. Phys. Chem. A 106 (2002) 7671. H. Naundorf, O. Kühn, in: A. Douhal, J. Santamaria (Eds.), Femtochemistry and Femtobiology, World Scientific, Singapore, 2002, p. 438. M. Fores, M. Duran, M. Sola, Chem. Phys. 260 (2000) 53. M. Fores, M. Duran, M. Sola, Chem. Phys. 234 (1998) 1. A. Simperler, W. Mikenda, Chem. Monthly 128 (1997) 969. M. Petkovi´c, O. Kühn, J. Phys. Chem. A 107 (2003) 8458. M. Petkovi´c, O. Kühn, Chem. Phys. 304 (2004) 91. M. Petkovi´c, O. Kühn, in: M.M. Martin, J.T. Hynes (Eds.), Ultrafast Molecular Events in Chemistry and Biology, Elsevier, Amsterdam, 2004, p. 181. K. Heyne, J. Stenger, J. Dreyer, E. Nibbering, T. Elsaesser, unpublished. R.L. Redington, R.L. Sams, J. Phys. Chem. A 106 (2002) 7494. A.C.P. Alves, J.M. Hollas, Mol. Phys. 23 (1972) 927. A.C.P. Alves, J.M. Hollas, Mol. Phys. 25 (1973) 1305. K. Tanaka, H. Honjo, T. Tanaka, H. Kohguchi, Y. Ohshima, Y. Endo, J. Chem. Phys. 110 (1999) 1969. R.L. Redington, T.E. Redington, J.M. Montgomery, J. Chem. Phys. 113 (2000) 2304. R.L. Redington, J. Chem. Phys. 113 (2000) 2319. R.K. Frost, F. Hagemeister, D. Schleppenbach, G. Laurence, T.S. Zwier, J. Phys. Chem. B 100 (1996) 16835. K. Giese, Multidimensional tunneling in hydrogen transfer reactions, Ph.D. Thesis, Freie Universität Berlin, Berlin, 2005. M.D. Coutinho-Neto, A. Viel, U. Manthe, J. Chem. Phys. 121 (2004) 9207. E.J. Heller, J. Phys. Chem. 99 (1995) 2625. K. Giese, H. Ushiyama, K. Takatsuka, O. Kühn, J. Chem. Phys. 122 (2005) 124307. M. Gruebele, Adv. Chem. Phys. 114 (2000) 193. K. Giese, D. Lahav, O. Kühn, J. Theor. Comp. Chem. 3 (2004) 567. M. Razavy, Quantum Theory of Tunneling, World Scientific, Hoboken NJ, 2003. C. Herring, Rev. Mod. Phys. 34 (1962) 631. M. Wilkinson, Physica D 21 (1986) 341. V.A. Benderskii, S.Y. Grebenshchikov, E.V. Vetoshkin, G.V. Mil’nikov, D.E. Makarov, Chem. Phys. 98 (1994) 3309.

K. Giese et al. / Physics Reports 430 (2006) 211 – 276 [142] [143] [144] [145] [146] [147] [148] [149] [150] [151] [152] [153] [154] [155] [156] [157] [158] [159] [160] [161] [162] [163]

273

V.A. Benderskii, E. Vetoshkin, I. Irgibaeva, H. Trommsdorff, Chem. Phys. 262 (2000) 393. E. Vetoshkin, V.A. Benderskii, S.Y. Grebenshchikov, L. von Laue, H.P. Trommsdorff, Chem. Phys. 219 (1997) 119. Z. Smedarchina, W. Siebrand, M.Z. Zgierski, J. Chem. Phys. 103 (1995) 5326. Z. Smedarchina, W. Siebrand, M.Z. Zgierski, J. Chem. Phys. 104 (1996) 1203. Z. Smedarchina, A. Fernandez-Ramos, W. Siebrand, J. Comput. Chem. 22 (2001) 787. Z. Smedarchina, A. Fernandez-Ramos, W. Siebrand, Chem. Phys. Lett. 395 (2004) 339. A. Garg, Am. J. Phys. 68 (2000) 430. G. Mil’nikov, H. Nakamura, J. Chem. Phys. 115 (2001) 6881. G.V. Mil’nikov, H. Nakamura, J. Chem. Phys. 122 (2005) 124311. G.V. Mil’nikov, K. Yagi, T. Taketsugu, H. Nakamura, K. Hirao, J. Chem. Phys. 119 (2003) 10. G.V. Mil’nikov, O. Kühn, H. Nakamura, J. Chem. Phys. 123 (2005) 074308. M.F. Herman, E. Kluk, Chem. Phys. 91 (1984) 27. W.H. Miller, J. Phys. Chem. A 105 (2001) 2942. S. Garashchuk, D. Tannor, Chem. Phys. Lett. 262 (1996) 477. V.A. Mandelshtam, M. Ovchinikov, J. Chem. Phys. 108 (1998) 9206. V. Guallar, V.S. Batista, W.H. Miller, J. Chem. Phys. 110 (1999) 9922. V. Guallar, V.S. Batista, W.H. Miller, J. Chem. Phys. 113 (2000) 9510. T. Yamamoto, W.H. Miller, J. Chem. Phys. 118 (2003) 2135. K. Giese, O. Kühn, J. Chem. Phys. 120 (2004) 4207. V.A. Mandelshtam, Progr. Nucl. Magn. Res. 38 (2001) 159. D.G. Truhlar, A.D. Isaacson, B.C. Garrett, in: M. Baer (Ed.), Theory of Chemical Reaction Dynamics, vol. IV, CRC, Boca Raton, 1985, p. 65. T.C. Allison, D.G. Truhlar, in: D.L. Thompson (Ed.), Modern Methods for Multidimensional Dynamics Computations in Chemistry, World Scientific, Singapore, 1998, p. 618. [164] R.A. Marcus, M.E. Coltrin, J. Chem. Phys. 67 (1977) 2609. [165] N. Shida, J. Almlöf, P.F. Barbara, J. Phys. Chem. 95 (1991) 10457. [166] V.A. Benderskii, D.E. Makarov, P.G. Grinevich, Chem. Phys. 170 (1993) 275. [167] V.A. Benderskii, D.E. Makarov, C.A. Wright, Adv. Chem. Phys. 88 (1994) 1. [168] N. Makri, W.H. Miller, J. Chem. Phys. 91 (1989) 4026. [169] Y. Guo, S. Li, D.L. Thompson, J. Chem. Phys. 107 (1997) 2853. [170] T.D. Sewell, Y. Guo, D.L. Thompson, J. Chem. Phys. 103 (1995) 8557. [171] Y. Guo, D.L. Thompson, J. Phys. Chem. A 106 (2002) 8374. [172] V. Guallar, B.F. Gherman, W.H. Miller, S.J. Lippard, R.A. Friesner, J. Am. Chem. Soc. 124 (2002) 3377. [173] K. Giese, O. Kühn, J. Chem. Theor. Comp. 2 (2006) 717. [174] E. Bosch, M. Moreno, J.M. Lluch, Chem. Phys. 159 (1992) 99. [175] M. Ben-Nun, T.J. Martinez, J. Phys. Chem. A 103 (1999) 6055. [176] K. Giese, H. Ushiyama, O. Kühn, Chem. Phys. Lett. 371 (2003) 681. [177] M. Razavy, A. Pimpale, Phys. Rep. 168 (1988) 305. [178] K. Takatsuka, H. Ushiyama, A. Inoue-Ushiyama, Phys. Rep. 322 (1999) 347. [179] J.-H. Lim, E.K. Lee, Y. Kim, J. Phys. Chem. A 101 (1997) 2233. [180] Y. Kim, J. Am. Chem. Soc. 118 (1996) 1522. [181] S. Miura, M.E. Tuckerman, M.L. Klein, J. Chem. Phys. 109 (1998) 5290. [182] H. Ushiyama, K. Takatsuka, J. Chem. Phys. 115 (2001) 5903. [183] K. Wolf, A. Simperler, W. Mikenda, Chem. Monthly 130 (1999) 1031. [184] C. Rambaud, H.P. Trommsdorff, Chem. Phys. Lett. 306 (1999) 124. [185] W.-T. Chan, I.P. Hamilton, Chem. Phys. Lett. 292 (1998) 57. [186] R. Xu, Y.J. Yan, O. Kühn, Eur. Phys. J. D 19 (2002) 293. [187] R. Rey, J.T. Hynes, J. Chem. Phys. 104 (1996) 2356. [188] C.P. Lawrence, J.L. Skinner, J. Chem. Phys. 117 (2002) 5827. [189] D. Borgis, J.T. Hynes, Chem. Phys. 170 (1993) 315. [190] M. Tuckerman, K. Laasonen, M. Sprik, M. Parrinello, J. Chem. Phys. 103 (1995) 150. [191] M.V. Vener, O. Kühn, J. Sauer, J. Chem. Phys. 114 (2001) 240. [192] A. Laio, J. VandeVondele, U. Röthlisberger, J. Chem. Phys. 116 (2002) 6941. [193] R. Rousseau, V. Kleinschmidt, U.W. Schmitt, D. Marx, Phys. Chem. Chem. Phys. 6 (2004) 1848. [194] R. Rousseau, V. Kleinschmidt, U.W. Schmitt, D. Marx, Angew. Chem. Int. Ed. 43 (2004) 4804. [195] J. Aaqvist, A. Warshel, Chem. Rev. 93 (1993) 2523. [196] J. Lobaugh, G.A. Voth, J. Chem. Phys. 104 (1996) 2056. [197] U.W. Schmitt, G.A. Voth, J. Phys. Chem. B 102 (1998) 5547. [198] R. Vuilleumier, D. Borgis, in: B.J. Berne, G. Ciccotti, D.F. Coker (Eds.), Classical and Quantum Dynamics in Condensed Phase Simulations, World Scientific, Singapore, 1998, p. 721. [199] S. Braun-Sand, M.H.M. Olsson, J. Mavri, A. Warshel, in: J.T. Hynes, H.-H. Limbach (Eds.), Handbook of Hydrogen Transfer, vol. 1: Physical and Chemical Aspects of Hydrogen Transfer, Wiley-VCH, Weinheim, 2006. [200] H.J.C. Berendsen, J. Mavri, J. Phys. Chem. 97 (1993) 13464.

274

K. Giese et al. / Physics Reports 430 (2006) 211 – 276

[201] A. Staib, D. Borgis, J.T. Hynes, J. Chem. Phys. 102 (1995) 2487. [202] P. Bala, P. Grochowski, B. Lesyng, J.A. McCammon, J. Phys. Chem. 100 (1996) 2535. [203] S. Hammes-Schiffer, J. Chem. Phys. 105 (1996) 2236. [204] S.R. Billeter, W.F. van Gunsteren, Comp. Phys. Commun. 107 (1997) 61. [205] S.Y. Kim, S. Hammes-Schiffer, J. Chem. Phys. 119 (2003) 4389. [206] N. Rösch, M.A. Ratner, J. Chem. Phys. 61 (1974) 3344. [207] B. Boulil, O. Henri-Rousseau, P. Blaise, Chem. Phys. 126 (1988) 263. [208] B. Boulil, J.-L. Dejardin, N.E. Ghandour, O. Henri-Rousseau, J. Mol. Struct. (Theochem) 314 (1994) 83. [209] G.N. Robertson, J. Yarwood, Chem. Phys. 32 (1978) 267. [210] A. Witkowski, M. Wójcik, Chem. Phys. 1 (1973) 9. [211] S. Bratos, J. Chem. Phys. 63 (1975) 3499. [212] S. Bratos, H. Ratajcak, J. Chem. Phys. 76 (1982) 77. [213] S. Mukamel, Principles of Nonlinear Optical Spectroscopy, Oxford, New York, 1995. [214] S. Bratos, J.C. Leicknam, J. Chem. Phys. 101 (1994) 4536. [215] S. Bratos, A. Laubereau, in: D. Hadži (Ed.), Theoretical Treatments of Hydrogen Bonding, Wiley, Chichester, 1997, p. 187. [216] Y. Tanimura, S. Mukamel, J. Chem. Phys. 99 (1993) 9496. [217] J.R. Schmidt, S.A. Corcelli, J.L. Skinner, J. Chem. Phys. 123 (2005) 044513. [218] D.W. Oxtoby, D. Levesque, J.-J. Weis, J. Chem. Phys. 68 (1978) 5528. [219] H.J. Bakker, J. Chem. Phys. 121 (2004) 10088. [220] C. Scheurer, P. Saalfrank, J. Chem. Phys. 104 (1996) 2869. [221] O. Kühn, Y. Tanimura, J. Chem. Phys. 119 (2003) 2155. [222] A.G. Redfield, Adv. Magn. Reson. 1 (1965) 1. [223] W.T. Pollard, A.K. Felts, R.A. Friesner, Adv. Chem. Phys. 93 (1996) 77. [224] N. Makri, J. Phys. Chem. A 102 (1998) 4144. [225] C. Meier, D.J. Tannor, J. Chem. Phys. 111 (1999) 3365. [226] R. Xu, Y. Yan, J. Chem. Phys. 116 (2002) 9196. [227] Q. Shi, E. Geva, J. Chem. Phys. 119 (2003) 12063. [228] N. Makri, J. Phys. Chem. B 103 (1999) 2823. [229] D. Xu, K. Schulten, Chem. Phys. 182 (1994) 91. [230] R. Zwanzig, J. Chem. Phys. 34 (1961) 1931. [231] R. Zwanzig, Nonequilibrium Statistical Mechanics, Oxford University Press, New York, 2001. [232] D.W. Oxtoby, Adv. Chem. Phys. 47 (1981) 487. [233] J.C. Owrutsky, D. Raftery, R.M. Hochstrasser, Annu. Rev. Phys. Chem. 45 (1994) 519. [234] S.A. Egorov, K.F. Everitt, J.L. Skinner, J. Phys. Chem. A 103 (1999) 9494. [235] C.P. Lawrence, J.L. Skinner, J. Chem. Phys. 119 (2003) 1623. [236] O. Kühn, H. Naundorf, Phys. Chem. Chem. Phys. 5 (2003) 79. [237] D. Reichman, R.J. Silbey, A. Suarez, J. Chem. Phys. 105 (1996) 10500. [238] Y. Yan, F. Shuang, R. Xu, J. Cheng, X.-Q. Li, C. Yang, H. Zhang, J. Chem. Phys. 113 (2000) 2068. [239] Y. Ohtsuki, Y. Fujimura, J. Chem. Phys. 91 (1989) 3903. [240] J.M. Jean, G.R. Fleming, J. Chem. Phys. 103 (1995) 2092. [241] A.M. Walsh, R.D. Coalson, Chem. Phys. Lett. 198 (1992) 293. [242] O. Kühn, V. Sundström, T. Pullerits, Chem. Phys. 275 (2002) 15. [243] K. Heyne, N. Huse, E.T.J. Nibbering, T. Elsaesser, Chem. Phys. Lett. 369 (2003) 591. [244] N. Huse, K. Heyne, J. Dreyer, E.T.J. Nibbering, T. Elsaesser, Phys. Rev. Lett. 91 (2003) 197401. [245] K. Heyne, N. Huse, J. Dreyer, E.T.J. Nibbering, T. Elsaesser, S. Mukamel, J. Chem. Phys. 121 (2004) 902. [246] B. Bagchi, Chem. Rev. 105 (2005) 3197. [247] D. Marx, M.E. Tuckerman, J. Hutter, M. Parrinello, Nature 397 (1999) 601. [248] M.E. Tuckerman, D. Marx, M. Parinello, Nature 417 (2002) 925. [249] R. Ludwig, Chem. Phys. Chem. 5 (2004) 1495. [250] K.R. Asmis, N.L. Pivonka, G. Santambrogio, M. Brümmer, C. Kaposta, D.M. Neumark, L. Wöste, Science 299 (2003) 1375. [251] N.I. Hammer, E.G. Diken, J.R. Roscioli, M.A. Johnson, E.M. Myshakin, K.D. Jordan, A.B. McCoy, X. Huang, J.M. Bowman, S. Carter, J. Chem. Phys. 122 (2005) 244301. [252] J. Sauer, J. Döbler, Chem. Phys. Chem. 6 (2005) 1706. [253] W. Kühlbrandt, Nature 406 (2000) 569. [254] H. Graener, G. Seifert, A. Laubereau, Phys. Rev. Lett. 66 (1991) 2092. [255] S. Woutersen, U. Emmerichs, H.J. Bakker, Nature 278 (1997) 658. [256] H.-K. Nienhuys, S. Woutersen, R.A. van Santen, H.J. Bakker, J. Chem. Phys. 111 (1999) 1494. [257] G.M. Gale, G. Gallot, F. Hache, N. Lascoux, S. Bratos, J.-C. Leicknam, Phys. Rev. Lett. 82 (1999) 1068. [258] J.C. Deak, S.T. Rhea, L.K. Iwaki, D.D. Dlott, J. Phys. Chem. A 104 (2000) 4866. [259] R. Rey, K.B. Moller, J.T. Hynes, Chem. Rev. 104 (2004) 1915. [260] R. Rey, K.B. Moller, J.T. Hynes, J. Phys. Chem. A 106 (2002) 11993. [261] J. Stenger, D. Madsen, P. Hamm, E.T.J. Nibbering, T. Elsaesser, Phys. Rev. Lett. 87 (2001) 027401.

K. Giese et al. / Physics Reports 430 (2006) 211 – 276 [262] [263] [264] [265] [266] [267] [268] [269] [270] [271] [272] [273] [274] [275] [276] [277] [278] [279] [280] [281] [282] [283] [284] [285] [286] [287] [288] [289] [290] [291] [292] [293] [294] [295] [296] [297] [298] [299] [300] [301] [302] [303] [304] [305] [306] [307] [308] [309] [310] [311] [312] [313] [314] [315] [316] [317] [318] [319] [320] [321]

275

J.D. Eaves, A. Tokmakoff, P.L. Geissler, J. Phys. Chem. A 109 (2005) 9424. J. Stenger, D. Madsen, P. Hamm, E.T.J. Nibbering, T. Elsaesser, J. Phys. Chem. A 106 (2002) 2341. C.J. Fecko, J.D. Eaves, J.J. Loparo, A. Tokmakoff, P.L. Geissler, Science 301 (2003) 1698. C.P. Lawrence, J.L. Skinner, J. Chem. Phys. 118 (2003) 264. R. Meyer, R.R. Ernst, J. Chem. Phys. 93 (1990) 5518. O. Brackhagen, O. Kühn, J. Manz, V. May, R. Meyer, J. Chem. Phys. 100 (1994) 9007. O. Brackhagen, C. Scheurer, R. Meyer, H.-H. Limbach, Ber. Bunsenges. Phys. Chem. 102 (1998) 303. N. Došli´c, K. Sundermann, L. González, O. Mo, J. Giraud-Girard, O. Kühn, Phys. Chem. Chem. Phys. 1 (1999) 1249. N. Došli´c, O. Kühn, Chem. Phys. 255 (2000) 247. H. Naundorf, K. Sundermann, O. Kühn, Chem. Phys. 240 (1999) 163. D. Borgis, G. Tarjus, H. Azzouz, J. Phys. Chem. 96 (1992) 3188. U.W. Schmitt, G.A. Voth, J. Chem. Phys. 111 (1999) 9361. S. Hammes-Schiffer, J.C. Tully, J. Chem. Phys. 101 (1994) 4657. J.-Y. Fang, S. Hammes-Schiffer, J. Chem. Phys. 110 (1999) 11166. H. Decornez, K. Drukker, S. Hammes-Schiffer, J. Phys. Chem. A 103 (1999) 2891. S.P. Webb, P.K. Agarwal, S. Hammes-Schiffer, J. Phys. Chem. B 104 (2000) 8884. S.R. Billeter, S.P. Webb, T. Iordanov, P.K. Agarwal, S. Hammes-Schiffer, J. Chem. Phys. 114 (2001) 6925. P.K. Agarwal, S.R. Billeter, S. Hammes-Schiffer, J. Phys. Chem. B 106 (2002) 3283. S. Mukamel, Annu. Rev. Phys. Chem. 51 (2000) 691. M. Cho, Phys. Chem. Comm. 5 (2002) 40. M. Khalil, N. Demirdöven, A. Tokmakoff, J. Phys. Chem. A 107 (2003) 5258. J.B. Asbury, T. Steinel, K. Kwak, S.A. Corcelli, C.P. Lawrence, J.L. Skinner, M.D. Fayer, J. Chem. Phys. 121 (2004) 12431. T. la Cour Jansen, T. Hayashi, W. Zhuang, S. Mukamel, J. Chem. Phys. 123 (2005) 114504. J.B. Asbury, T. Steinel, C. Stromberg, K.J. Gaffney, I.R. Piletic, A. Goun, M.D. Fayer, Chem. Phys. Lett. 374 (2003) 362. J.B. Asbury, T. Steinel, C. Stromberg, K.J. Gaffney, I.R. Piletic, M.D. Fayer, J. Chem. Phys. 119 (2003) 12981. A. Ishizaki, Y. Tanimura, J. Chem. Phys. 123 (1) (2005) 014503. J.M. Robertson, A.R. Ubbelohde, Proc. Roy. Soc. A 170 (1939) 222. A.R. Ubbelohde, K.J. Gallagher, Acta Crystallogr. 8 (1955) 71. M. Ichikawa, Act. Crystallogr. B 34 (1978) 2074. M. Ichikawa, J. Mol. Struct. 552 (2000) 63. S. Tanaka, Phys. Rev. B 42 (1990) 10488. H. Benedict, H.-H. Limbach, M. Wehlan, W.-P. Fehlhammer, N.S. Golubev, R. Janoschek, J. Am. Chem. Soc. 120 (1998) 2939. I.G. Shenderovich, H.-H. Limbach, S.N. Smirnov, P.M. Tolstoy, G.S. Denisov, N.S. Golubev, Phys. Chem. Chem. Phys. 4 (2002) 5488. H.-H. Limbach, M. Pietrzak, H. Benedict, P.M. Tolstoy, N.S. Golubev, G.S. Denisov, J. Mol. Struct. 706 (2004) 115. M. Tachikawa, M. Shiga, J. Am. Chem. Soc. 127 (2005) 11908. M.F. Shibl, M. Tachikawa, O. Kühn, Phys. Chem. Chem. Phys. 7 (2005) 1368. M.F. Shibl, M. Pietrzak, H.-H. Limbach, O. Kühn, 2006, submitted. J.E. Almlöf, Chem. Phys. Lett. 17 (1972) 49. T. Elsaesser, in: T. Elsaesser, H.J. Bakker (Eds.), Ultrafast Hydrogen Bonding Dynamics and Proton Transfer Processes in the Condensed Phase, Kluwer Academic Publishers, Dordrecht, 2002, p. 119. C. Crespo-Hernandez, P.M.H.B. Cohen, B. Kohler, Chem. Rev. 104 (2004) 1977. C. Chudoba, E. Riedle, M. Pfeiffer, T. Elsaesser, Chem. Phys. Lett. 263 (1996) 622. A. Douhal, Science 276 (1997) 221. S. Lochbrunner, A.J. Wurzer, E. Riedle, J. Chem. Phys. 112 (2000) 10699. S. Lochbrunner, E. Riedle, in Recent Research Developments in Chemical Physics, vol. 4, Transworld Research Network, Trivandrum, Kerala, 2003, p. 31. R. deVivie Riedle, V. DeWaele, L. Kurtz, E. Riedle, J. Phys. Chem. A 107 (2003) 10591. N. Agmon, J. Phys. Chem. A 109 (2005) 13. E.T.J. Nibbering, H. Fidder, E. Pines, Annu. Rev. Phys. Chem. 56 (2005) 337. T. Elsaesser, W. Kaiser, Chem. Phys. Lett. 128 (1986) 231. M. Rini, B.-Z. Magnes, E. Pines, E.T.J. Nibbering, Science 301 (2003) 349. O.F. Mohammed, D. Pines, J. Dreyer, E. Pines, E.T.J. Nibbering, Science 310 (2005) 83. M. Glasbeek, Isr. J. Chem. 39 (1999) 301. A. Douhal, S.K. Kim, A.H. Zewail, Nature 378 (1995) 260. S. Takeuchi, T. Tahara, J. Phys. Chem. A 102 (1998) 7740. K. Sakota, C. Okabe, N. Nishi, H. Sekiya, J. Phys. Chem. A 109 (2005) 5245. T. Schultz, E. Samoylova, W. Radloff, I.V. Hertel, A.L. Sobolewski, W. Domcke, Nature 306 (2005) 1765. S. Rice, M. Zhao, Optimal Control of Molecular Dynamics, Wiley, New York, 2001. M. Shapiro, P. Brumer, Principles of the Quantum Control of Molecular Processes, Wiley, Hoboken, 2003. T. Brixner, G. Gerber, Chem. Phys. Chem. 4 (2003) 418. I. Walmsley, H. Rabitz, Phys. Today August (2003) 43. O. Kühn, L. González, in: J.T. Hynes, H.-H. Limbach (Eds.), Handbook of Hydrogen Transfer, vol. 1: Physical and Chemical Aspects of Hydrogen Transfer, Wiley-VCH, Weinheim, 2006.

276 [322] [323] [324] [325] [326] [327] [328] [329]

K. Giese et al. / Physics Reports 430 (2006) 211 – 276 M. Grifoni, P. Hänggi, Phys. Rep. 304 (1998) 229. G.K. Paramonov, M.V. Korolkov, J. Manz, Adv. Chem. Phys. 101 (1997) 327. N. Došli´c, O. Kühn, J. Manz, Ber. Bunsenges. Phys. Chem. 102 (1998) 292. N. Došli´c, O. Kühn, J. Manz, K. Sundermann, J. Phys. Chem. A 102 (1998) 292. O. Kühn, Eur. Phys. J. D 6 (1999) 49. O. Kühn, Y. Zhao, F. Shuang, Y.-J. Yan, J. Chem. Phys. 112 (2000) 6104. E. Geva, J. Chem. Phys. 116 (2002) 1629. V.P. Maslov, M.V. Fedoriuk, Semiclassical Approximations in Quantum Mechanics, D. Reidel, Dordrecht, 1981.