Ground and asymmetric CO-stretch excited state

Ground and asymmetric CO-stretch excited state tunneling splittings in the formic acid dimer I. Matanović, N. Došlić, and O. Kühn Citation: J. Chem. Phys. 127, 014309 (2007); doi: 10.1063/1.2748048 View online: http://dx.doi.org/10.1063/1.2748048 View Table of Contents: http://jcp.aip.org/resource/1/JCPSA6/v127/i1 Published by the AIP Publishing LLC.

Additional information on J. Chem. Phys. Journal Homepage: http://jcp.aip.org/ Journal Information: http://jcp.aip.org/about/about_the_journal Top downloads: http://jcp.aip.org/features/most_downloaded Information for Authors: http://jcp.aip.org/authors

Downloaded 03 Oct 2013 to 202.116.1.148. This article is copyrighted as indicated in the abstract. Reuse of AIP content is subject to the terms at: http://jcp.aip.org/about/rights_and_permissions

THE JOURNAL OF CHEMICAL PHYSICS 127, 014309 共2007兲

Ground and asymmetric CO-stretch excited state tunneling splittings in the formic acid dimer I. Matanović and N. Došlić Department of Physical Chemistry, Rudjer Bošković Institute, 10000 Zagreb, Croatia

O. Kühna兲 Institut für Chemie und Biochemie, Freie Universität Berlin, D-14195 Berlin, Germany

共Received 26 March 2007; accepted 16 May 2007; published online 5 July 2007兲 There has been some controversy concerning the assignment of measured tunneling splittings for the formic acid dimer in the vibrational ground state and the asymmetric CO-stretching excited state. The discussion is intimately related to the question whether the fundamental excitation of the CO-vibration promotes or hinders tunneling. Here we will address this issue on the basis of a five-dimensional reaction space Hamiltonian which includes three large amplitude coordinates as well as two harmonic modes whose linear superposition reproduces the asymmetric CO-vibrational mode. Within density functional theory using the B3LYP functional together with a 6-311+ + G共3df , 3pd兲 basis set we obtain a ground state tunneling splitting which is about 2.4 larger than the one for the CO-stretching excited state. © 2007 American Institute of Physics. 关DOI: 10.1063/1.2748048兴 I. INTRODUCTION

Carboxylic acid dimers are considered to be prototypical for systems with intermolecular double hydrogen bonds. In particular, infrared absorption line splittings due to concerted double proton tunneling provide a serious challenge for theory. In this work we will focus on the formic acid dimer 共FAD兲 for which high-resolution gas-phase spectroscopy in the region of the asymmetric CO-stretching vibration is available. Studying deuterated FAD 共DCOOH兲2 in the 1245 cm−1 region Madeja and Havenith assigned a ground state tunneling splitting of ⌬0 = 0.002 86 cm−1 together with a vibrationally excited state splitting of ⌬1 = 0.009 99 cm−1.1 In other words, it was concluded that vibrational excitation of the CO-stretching mode promotes tunneling. It is important to note that the basis for this assignment has been the simultaneous observation of a- and b-type transitions. However, the experimental limitations as far as the intensity of the a-type transitions is concerned allowed for an alternative assignment as well. Here, the ground state splitting is ⌬0 = 0.0125 cm−1, while the CO-excited state splitting would be smaller 共⌬1 = 0.0031 cm−1兲, i.e., this mode would not promote tunneling. Very recent investigations of 共HCOOH兲2 by the same group in the 1225 cm−1 range allowed for an unambiguous assignment which points to a hindering of the tunneling by the asymmetric CO-stretching vibration.2 Specifically, it was found that ⌬0 = 0.0158 cm−1 and ⌬1 = 0.010 cm−1. This, of course, makes the alternative assignment for 共DCOOH兲2 rather likely as well. The first theoretical estimate of the ground state tunneling splitting of 共HCOOH兲2 given by Chang et al. was 0.3 cm−1.3 Using a three-dimensional reaction surface model Shida et al.4 later predicted a value of 0.004 cm−1. In a recent a兲

Electronic mail: [email protected]

0021-9606/2007/127共1兲/014309/7/$23.00

model study on the basis of an ab initio two-dimensional symmetric coupling Hamiltonian a value of 0.0037 cm−1 had been obtained.5 Further, a reduced dimensionality model focusing on the effect of tunneling in the OH-stretch region gave a ground state tunneling splitting of 0.3 cm−1 as well as 300 cm−1 splitting due to concerted tunneling upon excitation of the symmetric OH-stretching fundamental transition.6 The original assignment of Madeja and Havenith in the asymmetric CO-stretching region1 was recently supported by Luckhaus who obtained a splitting of 0.0013 cm−1 within a five-dimensional model of 共HCOOH兲2 using density functional theory 共DFT兲/B3LYP with a 6-31+ G共d兲 basis set. This model, however, did not include the CO-stretching vibration.7 It is important to mention this restriction since the original assignment of Ref. 1 had been questioned by Smedarchina et al. in Refs. 8 and 9 which put forward a ground state splitting of 0.0147/ 0.0149 cm−1 and a CO-excited state splitting of 0.004/ 0.0031 cm−1 for 共DCOOH兲2 / 共HCOOH兲2. The level of quantum chemistry has been a combination of MC-QCISD/3 and MCG3 calculations. In fact, Ref. 8 presents the only simulation so far, which includes the COstretching mode, such as to enable prediction of its splitting in the fundamental excitation range. One should note that this assignment is based on approximate instanton theory. Approximate versions of instanton theory are widely used, for instance, Tautermann et al. predicted the ground state tunneling splitting for the FAD to be 0.0022 cm−1 共Ref. 10兲 at the DFT/B3LYP, 6-31+ G共d兲 level. However, especially for vibrationally excited states, the general validity of approximate instanton theories has to be questioned. In Ref. 11 it has been shown for the FAD that, in contrast to approximate instanton theory, no stable semiclassical instanton solution for frequencies higher than 1000 cm−1 could be found if a numerically exact form of instanton theory is used.12 This work reported a ground state tunneling split-

127, 014309-1

© 2007 American Institute of Physics


014309-2

J. Chem. Phys. 127, 014309 共2007兲

Matanović, Došlić, and Kühn

ting of 0.0038 cm−1 共Ref. 11兲 using a DFT/B3LYP, 6-311+ + G共3df , 3pd兲 instanton with CCSD共T兲 energy corrections. The wide range of tunneling splittings reported so far points to the strong sensitivity with respect to the method for solving the Schrödinger equation as well as the employed level of quantum chemistry. In view of the discussion of the semiclassical calculations8,9,11 as well as the new experimental results2 we will present a quantum mechanical model which includes the relevant coordinates participating in the tunneling process, that is, in particular, a coordinate which is of asymmetric CO-stretching character. Multidimensional potential energy surfaces 共PESs兲 including large amplitude motions in a few coordinates are often described by means of a reaction path13 or reaction surface14–16 共for a review, see also Ref. 5兲. Here, an exact treatment of the large amplitude coordinates is combined with a harmonic expansion with respect to the remaining orthogonal degrees of freedom 共DOFs兲. The first application to a double proton transfer system has been demonstrated for FAD by Shida et al.4 This type of reaction surface formulation has the disadvantage of relying on, e.g., internal coordinates leading to couplings in the kinetic energy operator whose treatment often requires additional approximations. To cope with this problem it has been suggested to use an allCartesian form of the reaction surface Hamiltonian. In the original formulation by Ruf and Miller it was based on atomic coordinates,17 see also the applications in Refs. 18–20. While this has the drawback of a nonphysical separation of coordinates leading to large couplings, a formulation in terms of collective Cartesian coordinates has been proposed in Ref. 21. It has been stimulated by earlier work of Takada and Nakamura22 and the idea for defining suitable large amplitude coordinates has also been used in the context of a Shepard interpolation of the PES for malonaldehyde23 and for calculating nonlinear IR spectra of the same molecule.24 It should be noted that the use of Cartesian coordinates necessitates further approximations as far as the vibration-rotation coupling is concerned which is usually neglected. In this contribution we will address the assignment of the tunneling splittings in the ground and CO vibrationally excited states on the basis of an all-Cartesian reaction space 共CRS兲 Hamiltonian, a term which had been coined by Yagi et al.23 for the case of three large amplitude coordinates 共LACs兲. The choice of the LAC is outlined in Sec. II A; the relevant orthogonal harmonic DOFs are introduced in Sec. II B. Results of the numerical diagonalization of the CRS Hamiltonian are presented in Sec. III and the results of the paper are discussed in Sec. IV. II. REACTION SPACE HAMILTONIAN A. Large amplitude coordinates

It has been suggested by Takada and Nakamura22 that a reduced dimensional PES for symmetric proton transfer reactions can be constructed on the basis of the known geometries of the stationary points on the PES corresponding to the right-hand minimum 共XR兲, the transition state 共XTS兲, and

FIG. 1. Displacements along the vectors defining the five-dimensional reaction space model. In harmonic approximation the directions of the LAC correspond to the following D2h transition state normal modes: d1: O–H stretching mode 共b3g兲 with a frequency of 1265i cm−1, d2: O¯O stretching mode 共ag symmetry兲 with a frequency of 505 cm−1, and d3: CO2 rocking 共b3g symmetry兲 with a frequency of 232 cm−1. The displacements along the vectors qa 共b2u兲 and qs 共b1u兲 represent two components of the asymmetric CO-stretching mode 共bu兲 of the C2h minimum structure.

the left-hand minimum 共XL兲, where mass-weighted coordinates are used. Specifically, for Nat atoms one can define the following two 3Nat dimensional vectors d1 and d2 as d1 =

XR − XL , 兩XR − XL兩

共1兲

d2 =

XC − XTS , 兩XC − XTS兩

共2兲

where XC represents the center geometry XC = 21 共XR + XL兲.

共3兲

Notice that given a molecular symmetry transformation, i.e., permutation of atoms and rotation, with the property TXR = XL ,

共4兲

the basis vectors will transform as Td1 = − d1,

Td2 = d2 .

共5兲

The following calculations have been performed for 共HCOOH兲2 using the DFT/B3LYP level of theory as implemented in GAUSSIAN 03.25 Unless stated otherwise a 6-311 + + G共3df , 3pd兲 basis set has been used. For comparison also results for a 6-31+ G共d兲 basis set will be reported. The atomic displacements along the directions d1 and d2 are shown in Fig. 1. This figure highlights that d1 corresponds to the concerted double proton transfer, whereas d2 describes the symmetric rearrangement of the molecular scaffold associated with the transfer. The two vectors are constructed in such a way as to reproduce geometries and energetics of the stationary points on the PES, which therefore are contained in the spanned reaction plane. Apart from these extrema, however, the relevance of any point on the reaction plane is not obvious and has to be justified according to the problem at hand. In Ref. 21 it had been suggested that for describing tunneling a reasonable demand is that the intrinsic reaction path 共IRP兲 in full-dimensional space does not deviate appreciably from the reaction plane. In the upper panel of Fig. 2 we show the energy profile along the full-dimensional IRP together with the deviation between this IRP and its projection onto the two-dimensional 共2D兲 reaction plane. In passing we note that all geometries along the IRP have been rotated with respect to the XTS in order to satisfy


014309-3

J. Chem. Phys. 127, 014309 共2007兲

Tunneling splitting in the formic acid dimer

as explained in the following section, it is obvious that the more anharmonicity one can capture with suitable large amplitude coordinates, the better will be the performance of the harmonic approximation for the remaining DOF. Figure 2 suggests to construct a third LAC from the two symmetrically related geometries XPL and XPR for which ␴共2兲共s兲 has its maximum. We define d3 according to d3 =

XPL − XPR , 兩XPL − XPR兩

共7兲

which is subsequently Gram-Schmidt orthogonalized to 共d1 , d2兲 and which transforms as 共8兲

Td3 = − d3 .

Following Yagi et al.23 the subspace 共d1 , d2 , d3兲 will be called reaction space. The atomic displacements along this direction are shown in the middle panel of Fig. 1. One notices that d3 involves a kind of wagging motion of the formic acid monomers accompanying double proton transfer. The rms deviation ␴共3兲共s兲 of the IRP with respect to the reaction space 共d1 , d2 , d3兲 is shown in Fig. 2. As anticipated the IRP lies approximately in the three-dimensional 共3D兲 reaction space; the maximum rms difference is as small as 7.3⫻ 10−3 Å at s = ± 2.1a0 amu1/2. The energy profile along the projected IRP is almost indistinguishable from that of the full IRP with the largest deviation being about 25 cm−1 共not shown兲. At this point it can be safely assumed that remaining deviations are part of the harmonic valley around the IRP. Their importance will be scrutinized in the following section. FIG. 2. Upper panel: Energy profile along the unconstrained IRP 共solid squares兲. The barrier height is 2280 cm−1. Energy difference between the IRP geometries and their projections on a reaction plane spanned by vectors d1 and d2 共crosses兲. Lower panel: Crosses—root mean squared 共rms兲 difference ␴n=2共s兲 between the IRP geometries 共X共s兲兲 and its projection on the reaction plane spanned by vectors d1 and d2. Solid dots—same for the reaction space 共n = 3兲 spanned by d1, d2, and d3.

the Eckart conditions.26 While this guarantees the uniqueness of the choice of XR and XL, it does not ensure rotational invariance of an arbitrary point on the reaction plane 共for a more detailed discussion see also the Appendix of Ref. 21兲. From Fig. 2 it is clear that the deviation of the IRP from the reaction plane is appreciable. This can be further quantified by inspecting the root mean squared 共rms兲 difference between geometries along the IRP parametrized by s, X共s兲, and its projection on the reaction plane. For the case of an n-dimensional subspace we define the rms deviation as

冏冑冏冏冉

␴共n兲共s兲 =

1

Nat

M−1/2

n

冊

冏

1 − 兺 didTi 共X共s兲 − XTS兲 , i=1

共6兲

where M is the diagonal mass matrix. For the present case of n = 2 this function is plotted in the lower panel of Fig. 2. By definition, the rms difference vanishes at the geometries of the transition state 共s = 0兲 and the minimum 共s = 3.3a0 amu1/2兲. It takes its largest value of 5.1⫻ 10−2 Å at s = ± 1.1a0 amu1/2, where the energy difference exceeds 1600 cm−1. Although the deviations from the IRP will approximately be described by displaced harmonic oscillators

B. Harmonic expansion and reduced dimensional reaction space Hamiltonian

A full-dimensional potential can be obtained by harmonic expansion with respect to the 3D reaction space. For this purpose it is advantageous to define the actual normal mode displacements with respect to the saddle point geometry which has the highest symmetry, D2h. Denoting the mass-weighted force constant matrix at the transition state as K共f兲共XTS兲 one readily obtains the set of 3Nat − 9 = 21 normal modes 共three reaction coordinates as well as translations and rotations projected out兲兵Qi其 which by construction will be symmetric or antisymmetric with respect to T. Up to quadratic order in Q one arrives at the following expression for the PES:5 V共d1,d2,d3,Q兲 = V共d1,d2,d3,Q = 0兲 − F共p兲共d1,d2,d3兲Q + 21 QK共p兲共d1,d2,d3兲Q,

共9兲

where F共p兲 and K共p兲 are the force and Hessian, respectively, projected onto the subspace of the normal mode vectors Qi. The remaining deviations of the IRP from the reaction space 共d1 , d2 , d3兲 can now be analyzed in terms of linear displacements of the oscillator coordinates ⌬Qi. In order to quantify these displacements it has been suggested21 to relate them to the width of the respective vibrational ground states. For all modes we find that the maximum displacement is smaller than 15% of the variance of the ground state Gaussian wave function, i.e., the numerical value of the wave function has dropped to about 99% of its maximum only such


014309-4

J. Chem. Phys. 127, 014309 共2007兲


TABLE I. Matrix which transforms the 共d1 , d2 , d3 , qa兲 and 共d1 , d2 , d3 , qs兲 . vectors to reduced 4D normal modes of the minima, Y共4D兲 j

共f兲 FIG. 3. Absolute values of the overlaps 兩Y共f兲 11 · Qi兩, where Y11 denotes the CO asymmetric stretching mode with symmetric 共b1u, solid兲 and antisymmetric 共b2u, dashed兲 modes.

that the effect will be barely noticeable. Therefore, as far as the proper description of the IRP itself is concerned, there is no need to include any oscillator DOF. With this background we are in the position to choose additional oscillator coordinates targeted to describe the asymmetric CO-stretching vibration. First, we notice that the assignment of this vibration refers to one of the minima, i.e., to the C2h configuration. Denoting the set of normal modes of the full Hessian at the minima, K共f兲共XMIN兲, as 兵Y 共f兲 k 其 one 共f兲 finds the asymmetric CO-stretching vibration as mode Y 11 −1 with the harmonic frequency of 1260 cm . In Fig. 3 we 共f兲 show the overlap between the vectors Qi and Y11 , i.e., 共f兲 Y11 · Qi. Appreciable overlaps exist with the symmetric modes 共b1u, in parentheses harmonic frequencies兲: Q9 共1397 cm−1兲, Q13 共1235 cm−1兲, and Q16 共804 cm−1兲 and the antisymmetric modes 共b2u兲: Q3 共1725 cm−1兲, Q6 共1566 cm−1兲, and Q7 共1403 cm−1兲. In the following we will make use of the freedom that any linear combination of modes Qi is equivalent to define two new modes qa and qs as a linear combination of the antisymmetric and symmetric modes, respectively. The coefficients in the linear combination are chosen in such a way 共f兲 that the overlaps 兩Y11 · Qi兩 are reproduced when calculating 共f兲兩Y11 · qa/s兩. Specifically, we have used qa = 冑2 / 5Q3 +冑1 / 5Q6 + 冑2 / 5Q7 and qs = 冑9 / 11Q9 − 冑1 / 11Q13 + 冑1 / 11Q16. Needless to say that apart from reproducing the overlaps these combinations are rather arbitrary. By construction these two new modes transform as Tqs = qs,

Tqa = − qa .

共10兲

The displacements along the vectors qa and qs are shown in the right part of Fig. 1. In order to appreciate the significance of the set of coordinates 共d1 , d2 , d3 , qa , qs兲 it is convenient to introduce a set of normal modes relating to a certain minimum in a reduced space as well. If the dimensionality is ˜n this set of reduced ˜ D兲其. The decomposition normal modes will be labeled as 兵Y 共n k of these reduced normal modes is first analyzed for the separate subspaces 共d1 , d2 , d3 , qa兲 and 共d1 , d2 , d3 , qs兲 in Table I. It can be seen that for both cases, the reduced mode Y 共4D兲 2

Mode

␻ 共cm−1兲

d1

d2

d3

qa

Y共4D兲 1 Y共4D兲 2 Y共4D兲 3 Y共4D兲 4

2709 1527 299 190

−0.54 ⬇0 0.63 −0.56

0.19 ⬇0 −0.55 0.81

0.82 ⬇0 0.54 0.18

⬇0 0.998 ⬇0 ⬇0

Mode

␻ 共cm−1兲

d1

d2

d3

qs

Y共4D兲 1 Y共4D兲 2 Y共4D兲 3 Y共4D兲 4

2709 1297 299 190

−0.54 ⬇0 0.63 −0.56

0.19 ⬇0 −0.55 0.81

0.82 ⬇0 0.54 0.18

⬇0 0.998 ⬇0 ⬇0

exclusively involves only mode qa or qs, while the directions 共4D兲共4D兲 Y共4D兲 are linear combination of the reaction 1 , Y3 , and Y4 space vectors d1, d2, and d3. This pattern reflects the fact that in harmonic approximation the coupling between the qa or qs mode and the d1 , d2 , d3 coordinates is forbidden by symmetry 共b1u 丢 b3g ⫽ ag or b2u 丢 b3g ⫽ ag兲. In order to establish the connection to the experimental assignment in full-dimensional space it is further necessary to inspect the relation between the reduced and full normal modes. Following Ref. 21 we define the Nat ⫻ ˜n matrix B that transforms between the reduced and the full-dimensional space, i.e., B = 共d1d2d3qs ¯ 兲,

共11兲

which allows us to calculate the overlap ˜ D兲共n T p jk = 共Y共f兲 . j 兲 BYk

共12兲

These overlaps are compiled for the two four-dimensional 共4D兲 cases in Table II. Reduced mode Y 共4D兲 has the largest 1 overlap with the OH stretching mode 关full normal mode Y 共f兲 4 兴 with a frequency of 2709 cm−1. The reduced modes Y 共4D兲 3 and Y 共4D兲 have the largest overlap with the low frequency 4 共f兲共f兲 skeletal vibrations 关full normal modes Y 21 and Y 23 兴 which are the O¯O stretching mode with a frequency of 207 cm−1 and a CO2 rocking mode with the frequency of 190 cm−1 at the minima. It can be seen from Table II that the reduced 共4D兲 −1 −1 modes Y 共4D兲 2 共1527 cm 兲 and Y 2 共1297 cm 兲 cannot repro共f兲兴. They duce the C–O stretching mode of the minima 关Y 11 共f兲兴 include both the asymmetric C–O 关full normal mode Y 11 and asymmetric C v O stretching modes 关full normal mode Y 共f兲 5 兴. This is expected as both vectors qa and qa are constructed from the eigenvectors of a projected force constant matrix at the transition state for which the difference between the C–O 共donor兲 and C v O 共acceptor兲 bonds does not exist. To reproduce the asymmetric C–O stretching mode of the minima both of the modes have to be included in the calculation that leads to a five-dimensional 共5D兲 problem. In Tables III and IV we show the decomposition of the 5D subspace 共d1 , d2 , d3 , qa , qs兲 into its normal modes defined at the minimum as well as the important overlaps between reduced 5D normal modes and full normal modes. Table III reveals the same separation between the large amplitude and harmonic coordinates. Further, the mixing between qa and qs


014309-5

J. Chem. Phys. 127, 014309 共2007兲


TABLE II. Overlaps p jk between the reduced normal mode vectors Y共4D兲 j from Table I with the full normal modes Yk共f兲 of the minima. Harmonic frequencies of the full normal modes of the minima are calculated on a B3LYP/ 6-311+ + G共3df , 3pd兲 level of theory.

TABLE IV. Overlaps p jk between the reduced normal modes Y共5D兲 and Y共5D兲 2 3 from Table III with the full normal modes Yk共f兲 of the minima. It can be seen represents the asymmetric CO stretching mode. that mode Y共5D兲 3

␻ 共cm−1兲

Full mode

␻ 共cm−1兲

p jk

Y共5D兲 2

1616

Y共f兲 5 Y共f兲 8 Y共f兲 10 Y共f兲 11 Y共f兲 17

1770 1452 1403 1259 725

−0.74 0.23 −0.54 0.11 0.23

Y共5D兲 3

1245

Y共f兲 10 Y共f兲 11 Y共f兲 17

1403 1259 725

0.11 0.97 −0.18

Mode Mode

␻ 共cm−1兲

Mode

␻ 共cm−1兲

p jk

Y共4D兲 1

2709

Y共f兲 2 Y共f兲 4 Y共f兲 6 Y共f兲 9

3071 3050 1698 1403

0.31 −0.78 −0.46 −0.25

Y共4D兲 2

1527


1770 1452 1403 1259

0.64 −0.24 0.41 −0.59

Y共4D兲 2

1297


1770 1403 1259 725

−0.39 −0.33 0.75 0.41

Y共4D兲 3

299

Y共f兲 21 Y共f兲 23

207 176

0.83 0.53

Y共4D兲 4

190

Y共f兲 21 Y共f兲 23

207 176

−0.54 0.84

is obvious. Turning to Table IV we notice that the reduced 共f兲 mode Y 共5D兲 has a 97% overlap with the target mode Y 11 3 which allows us to conclude that the present 5D model comprises the subspace relevant for explanation of the experimental findings. Putting Tables III and IV together we can also state that the asymmetric C–O stretching mode can be written as a linear combination of the qa and qs modes, i.e., Y 共5D兲 = −0.86qs − 0.51qa. 3 The present model can be cast into the form of the following 5D reaction space Hamiltonian 共q = 共qa , qs兲兲: 3

H=−

2

1 ⳵2 1 ⳵2 − 兺兺 + V共d1,d2,d3,q = 0兲 2 i=1 ⳵d2i 2 k=1 ⳵q2k

− F共p兲共d1,d2,d3兲q + 21 qK共p兲共d1,d2,d3兲q.

共13兲

Focusing on this relevant 5D subspace we set all other harmonic coordinates equal to zero. The stationary Schrödinger equation for this Hamiltonian has been solved numerically, for the 5D case but for reference also for 4D and 3D 共no harmonic modes兲 cases. The potential energy surface has been obtained from 1500 共3000 when using symmetry兲 single point energies, gradients, and Hessians on a 共d1 , d2 , d3兲 grid. For this TABLE III. Matrix which transforms the directions 共d1 , d2 , d3 , qs , qa兲 to a 5D space of reduced normal modes.

purpose density functional theory using B3LYP/ 6-311 + + G共3df , 3pd兲 has been employed. The modified Shepard method27 was used for the interpolation of the PES in all directions. The grid parameters for the 3D, 4D, and 5D calculations are given in Table V. Eigenvalues and the eigenvectors of the 3D and 4D problems were obtained using the Fourier grid Hamiltonian method.28 For the diagonalization we used the Lanczos method as implemented in an ARPACK suits of programs.29 Eigenvalues and the eigenvectors of the 5D Hamiltonian 关Eq. 共13兲兴 were obtained by expressing the 5D eigenvectors in the basis of the lower dimensional wave functions ␾共3D兲共d1 , d2 , d3兲 and ␾共2D兲共qa , qs兲: i j

共d1,d2,d3,qa,qs兲 ⌿共5D兲 k 共d1,d2,d3兲␾共2D兲共qa,qs兲, = 兺 cijk␾共3D兲 i j

共14兲

ij

共d1 , d2 , d3兲关i = 0共±兲 , 1共±兲 , . . . 兴 and ␾共2D兲共qa , qs兲 where ␾共3D兲 i j 关j = 共Na , Ns兲; Na = 0 , 1 , . . ., Ns = 0 , 1 , . . .兴 represent the eigenvectors of a 3D and a 2D problem with the potentials V共d1 , d2 , d3兲 and V共qa , qs兲共all other coordinates set to zero兲. For the summation we found that convergence with respect to the target transitions can be achieved by including the lowest 12 basis functions from each set. Below we will only report on transition energies; IR intensities of the transitions were not calculated since it was not attempted to obtain an appropriate dipole moment surface. TABLE V. Grid parameters for the eigenvalue calculation on 3D 共d1 , d2 , d3兲, 4D 共d1 , d2 , d3 , qs兲, and 5D 共d1 , d2 , d3 , qs , qa兲 PESs. For every DOF, the grid extends from xmin to xmax with n points. The convergence of the grid parameters has been checked by means of various test calculations in different dimensions.

Mode

␻ 共cm−1兲

d1

d2

d3

qs

qa

DOF

n

xmin 共a0 amu1/2兲

xmax 共a0 amu1/2兲

Y共5D兲 1 Y共5D兲 2 Y共5D兲 3 Y共5D兲 4 Y共5D兲 5

2709 1616 1245 299 190

0.54 ⬇0 ⬇0 0.63 −0.56

−0.19 ⬇0 ⬇0 −0.55 −0.81

−0.82 ⬇0 ⬇0 0.54 −0.18

⬇0 0.51 −0.86 ⬇0 ⬇0

⬇0 −0.86 −0.51 ⬇0 ⬇0

d1 d2 d3 qa qs

35 31 23 21 21

−3.5 −0.6 −1.4 −0.38 −0.42

3.5 5.0 1.4 0.38 0.42


014309-6

J. Chem. Phys. 127, 014309 共2007兲


TABLE VI. Ground state and excited state tunneling splittings calculated on various potentials for which the grid parameters are given in Table V. Model

Transition

3D 共d1 , d2 , d3兲 ground state

共3D兲共3D兲 → ⌿0共−兲 ⌿0共+兲

4D 共d1 , d2 , d3 , qa兲 ground state qa mode 4D 共d1 , d2 , d3 , qs兲 ground state qs mode 5D 共d1 , d2 , d3 , qs , qa兲 ground state qs + qa mode

␻ij 共cm−1兲

⌬i / cm−1

0.197

共4D兲共4D兲 → ⌿0共−兲 ⌿0共+兲共4D兲共4D兲 ⌿0共+兲 → ⌿1共−兲共4D兲共4D兲 ⌿0共+兲 → ⌿1共+兲

1571.601 34 1571.713 62

共4D兲共4D兲 ⌿0共+兲 → ⌿0共−兲共4D兲共4D兲 ⌿0共+兲 → ⌿1共−兲共4D兲共4D兲 ⌿0共+兲 → ⌿1共+兲

1373.441 7 1373.640 39

共5D兲共5D兲 ⌿0共+兲 → ⌿0共−兲共5D兲共5D兲 ⌿0共+兲 → ⌿1共+兲共5D兲共5D兲 ⌿0共+兲 → ⌿1共−兲

1284.688 51 1284.754 01

0.162 0.112

0.197 0.199

of 1245 cm−1 obtained in harmonic approximation 共see Table III兲 indicating some reaction space mediated anharmonic coupling. The assignment of the CO-stretching fundamental transition has been based on visual inspection of the 2D projection of the 5D wave function. The splitting is obtained as ⌬1 = 0.066 cm−1 which is 2.4 times smaller value in the ground state 共0.155 cm−1兲. As far as the nature of these excited states is concerned it is interesting to inspect the expansion coefficients for the 5D wave function, as defined in Eq. 共14兲. It turns out that the excited state wave functions can approximately be written as 共3D兲共5D兲共3D兲共2D兲 = 关0.87␾0共±兲 + 0.10␾1共±兲兴␾01 ⌿1共±兲共3D兲共3D兲共2D兲 + 关0.47␾0共⫿兲 + 0.08␾1共⫿兲兴␾10 .

0.155 0.0655

III. GROUND AND VIBRATIONALLY EXCITED STATE TUNNELING SPLITTINGS

The results of the matrix diagonalization for the different models are summarized in Table VI. The tunneling splitting of the ground state, i.e., the vibrational transition frequency from the symmetric 共⫹兲 to the antisymmetric 共⫺兲 ground state, is determined as the difference between the two lowest eigenvalues of the Hamiltonian 关Eq. 共13兲兴. Comparing the ground state splittings for the different models we notice that the effect of the dimensionality is not very pronounced. For the 3D case it is 0.197 cm−1 which is almost unchanged upon inclusion of the qs in the 4D model 共0.197 cm−1兲. The value is somewhat reduced for the case of the antisymmetric mode qa 共0.162 cm−1兲 which hints to an anharmonic coupling to the large amplitude coordinates which cannot be seen from the harmonic values in Table II. The splitting is further reduced for the 5D model which gives ⌬0 = 0.155 cm−1. In passing we note that the actual ground state 共zero-point兲 energy is 2976 cm−1, i.e., above the barrier 共2280 cm−1兲. In other words, similar to other systems with comparable barrier heights 共see, e.g., Refs. 30 and 31兲 the splitting at least within the present 5D model is likely to be a consequence of dynamical tunneling. Consideration of tunneling splitting upon CO-stretching excitation is only meaningful in the 5D case since only here both qa and qs modes are mixed giving rise to the correct symmetry of the asymmetric CO-stretching vibration around the PES minima. This is also supported by the frequencies of the respective fundamental transitions. In fact, the 1285 cm−1 obtained for the 5D model are in reasonable agreement with the experimental value of 1245 cm−1 共Ref. 1兲 given the circumstance that we did not attempt to introduce any scaling accounting for the deficiencies of DFT. A further source of error comes from the all-Cartesian formulation which cannot account for vibration-rotation coupling. We also note that the 5D transition frequency of 1285 cm−1 is larger than the value

共15兲

In other words, this is the anticipated combination between the qa and qs fundamental excitations, but it contains some contribution from the reaction space wave function 共3D兲 ␾1共±兲共d1 , d2 , d3兲 as well. From the perspective of the asymmetric CO-stretching mode this mixing is responsible for the shift of the fundamental frequency as mentioned above. From the reaction space tunneling point of view, this mixing causes the change in the tunneling splitting. Inspecting the 共3D兲共d1 , d2 , d3兲 one finds that it nature of the wave function ␾1共±兲 involves excitation along both d1 and d2 directions. In the spirit of an anharmonic coupling model as put forward in Ref. 8 the interaction Hamiltonian would be proportional at least to d1d2qaqs which is possible since b3g 丢 ag 丢 b2u 丢 b1u = ag. In passing we note that in Ref. 8 only a third order coupling had been considered which if translated to the present model would roughly correspond to a term ⬀d1qaqs. IV. DISCUSSION

A five-dimensional all-Cartesian reaction space Hamiltonian has been presented for the concerted double proton transfer in the formic acid dimer. It comprises three large amplitude coordinates coupled to two harmonic degrees of freedom. While the selection of large amplitude coordinates has been based on the comparison with intrinsic reaction path calculations, the harmonic degrees of freedom were chosen such as to describe the asymmetric CO-stretching vibration. The calculated tunneling splittings are ⌬0 = 0.155 cm−1 for the ground state and ⌬1 = 0.0655 cm−1 for the COvibrationally excited state. In other words the excitation of the asymmetric CO vibration hinders tunneling. This conclusion is in accord with the results of approximate instanton theory put forward by Smedarchina et al.8 which had been the only account giving both splittings so far. In their approach the adiabatic treatment of a third order anharmonic coupling of their tunneling coordinate 共essentially the present d1兲 to two transition state modes 共roughly the present qa and qs兲 led to an increased effective mass for the tunneling particle as well as to a higher effective barrier, both effects reducing ⌬1. The physical argument of that paper has been that “the tunneling process displaces the C–O and C v O oscillators” thus exerting a “drag on the tunneling” such that


014309-7

J. Chem. Phys. 127, 014309 共2007兲


excitation of these modes should slow down tunneling. Within the present 5D quantum approach we showed how the decrease of the tunneling splitting emerges from the mixing of 3D reaction space and 2D CO-oscillator wave functions. Finally, we have to comment on the absolute numbers obtained for the tunneling splittings in the present 5D model. Adopting the new experimental assignment for 共HCOOH兲2 共ground state: 0.0158 cm−1, excited state: 0.010 cm−1兲,2 we find the ratio of the splittings to be 1.58 which is reasonably close to the present theoretical value of 2.4. Only the absolute values of the splittings are calculated to be about an order of magnitude larger. This fact we attribute mostly to the failure of density functional theory to accurately predict the reaction barrier 共see also discussion in Ref. 8兲. Indeed in Ref. 11 it has been shown for 共DCOOH兲2 that DFT/B3LYP with a 6-311+ + G共3df , 3pd兲 basis set gives a ground state tunneling splitting within instanton theory of 0.17 cm−1, a value which is rather close to the present 0.155 cm−1. The barrier in this case has been 2280 cm−1. Including CCSD共T兲 corrections to the energy the barrier increases to 2837 cm−1 and the instanton tunneling splitting dropped to 0.0038 cm−1 which, in view of the new assignment,2 is slightly below the experimental value. In order to further highlight the issue of the dependence on the level of quantum chemistry we have repeated the calculation using DFT/B3LYP with a 6-31+ G共d兲 basis set that gives ⌬0 = 0.0032 cm−1 and ⌬1 = 0.0012 cm−1, which can be attributed mostly to the higher barrier of 2933 cm−1. First, we notice that ⌬0 is close to the value obtain for another 5D model but with the same basis set 共0.0013 cm−1兲,7 suggesting that both models describe a similar region of the PES as far as the vibrational ground state is concerned. However, these values are an order of magnitude below the experimental ones.2 Finally, ⌬0 is also close to the instanton result for 共DCOOH兲2 which was 0.0077 cm−1.11 However, the latter calculation also suggests that there might be substantial error cancellation at work. Indeed, incorporating the CCSD共T兲 correction to the energy with the same basis set gave a tunneling splitting as small as 6 ⫻ 10−5 cm−1 共barrier increases to 4515 cm−1兲. Thus in view of the trend of the instanton results with respect to the quantum chemical method the present use of a 6-311+ + G共3df , 3pd兲 basis set appears to better justified even though it does not give the correct absolute numbers. In any case, the more important point is that both basis sets predict the ordering ⌬1 ⬍ ⌬0. Summing up, we have presented a five-dimensional quantum model which qualitatively reproduced the experimental observation that excitation of the asymmetric COstretching fundamental reduces the tunneling splitting in the formic acid dimer.

ACKNOWLEDGMENTS

This work has been supported by the MZOŠ of Croatia under Project No. 098-0352851-2921. Two of the authors 共N.D. and I.M.兲 gratefully acknowledge support from the Humboldt Foundation. The authors thank Dr. D. Luckhaus for his insightful comments on the manuscript and Professor M. Havenith for sharing the results of Ref. 2 prior to publication. F. Madeja and M. Havenith, J. Chem. Phys. 117, 7162 共2002兲. M. Ortlieb and M. Havenith, J. Phys. Chem. A 共to be published兲. 3 Y.-T. Chang, Y. Yamaguchi, W. H. Miller, and H. F. Schäfer III, J. Am. Chem. Soc. 109, 7245 共1987兲. 4 N. Shida, P. F. Barbara, and J. Almlöf, J. Chem. Phys. 94, 3633 共1991兲. 5 K. Giese, M. Petković, H. Naundorf, and O. Kühn, Phys. Rep. 430, 211 共2006兲. 6 O. Kühn, M. V. Vener, and J. M. Bowman, Chem. Phys. Lett. 349, 562 共2001兲. 7 D. Luckhaus, J. Phys. Chem. A 110, 3151 共2005兲. 8 Z. Smedarchina, A. Fernández-Ramos, and W. Siebrand, Chem. Phys. Lett. 395, 339 共2004兲. 9 Z. Smedarchina, A. Fernández-Ramos, and W. Siebrand, J. Chem. Phys. 122, 134309 共2005兲. 10 C. S. Tautermann, A. F. Voegele, and K. R. Liedl, J. Chem. Phys. 120, 631 共2004兲. 11 G. V. Mil’nikov, O. Kühn, and H. Nakamura, J. Chem. Phys. 123, 074308 共2005兲. 12 G. V. Mil’nikov and H. Nakamura, J. Chem. Phys. 122, 124311 共2005兲. 13 W. H. Miller, N. C. Handy, and J. E. Adams, J. Chem. Phys. 72, 99 共1980兲. 14 T. Carrington and W. H. Miller, J. Chem. Phys. 81, 3942 共1984兲. 15 T. Carrington and W. H. Miller, J. Chem. Phys. 84, 4364 共1986兲. 16 N. Shida, P. F. Barbara, and J. E. Almlöf, J. Chem. Phys. 91, 4061 共1989兲. 17 B. A. Ruf and W. H. Miller, J. Chem. Soc., Faraday Trans. 2 84, 1523 共1988兲. 18 H. Naundorf, J. A. Organero, A. Douhal, and O. Kühn, J. Chem. Phys. 110, 11286 共1999兲. 19 M. Petković and O. Kühn, J. Phys. Chem. A 107, 8458 共2004兲. 20 K. Giese, D. Lahav, and O. Kühn, J. Theor. Comput. Chem. 3, 567 共2004兲. 21 K. Giese and O. Kühn, J. Chem. Phys. 123, 054315 共2005兲. 22 S. Takada and H. Nakamura, J. Chem. Phys. 102, 3977 共1995兲. 23 K. Yagi, T. Taketsugu, and K. Hirao, J. Chem. Phys. 115, 10647 共2001兲. 24 T. Hayashi and S. Mukamel, J. Phys. Chem. A 107, 9113 共2003兲. 25 M. J. Frisch, G. W. Trucks, H. B. Schlegel et al., GAUSSIAN 03, Revision B.04, Gaussian Inc., Wallingford, CT, 2004. 26 W. Kohn, E. S. J. Robles, C. F. Logan, and P. Chen, J. Phys. Chem. 97, 4936 共1993兲. 27 T. Taketsugu, N. Watanabe, and K. Hirao, J. Chem. Phys. 111, 3410 共1999兲. 28 J. Stare and G. G. Balint-Kurti, J. Phys. Chem. A 107, 7204 共2003兲. 29 D. Sorensen, Tutorial: Implicitly Restarted Arnoldi/Lanczos Methods for Large Scale Eigenvalue Calculations 共Rice University, Houston, TX, 1995兲. 30 K. Giese, H. Ushiyama, K. Takatsuka, and O. Kühn, J. Chem. Phys. 122, 124307 共2005兲. 31 I. Matanović and N. Došlić, J. Phys. Chem. A 109, 4185 共2005兲. 1 2
