Geometric Isotope Effect on Low Barrier Hydrogen-Bonding Systems of Acetic Acid Dimer, Formic Acid Dimer, and Their Anionic Clusters by Using the Multi-Component Molecular Orbital Method

Masato KANEKO, Taro UDAGAWA and Masanori TACHIKAWA


1 Introduction

In some hydrogen bonding systems, there is small displacement of geometry due to deuterization of shared hydrogen atoms, which is referred to the geometric isotope effect (GIE). This effect is well known to be present in hydrogen-bonded ferroelectric crystal systems. For example, it is believed that large isotope effect in the phase transition temperature of hydrogen-bonded ferroelectric crystals such as potassium dihydrogen phosphate (KDP) is due to a small structural modification of the crystal by deuterization [1]. On the other hand, the GIE is not well known in hydrogen-bonded molecular clusters such as low barrier hydrogen-bonded (LBHB) [2] systems with carboxylic groups, which often play some fundamental roles in many disciplines across physical, chemical, and biological features [3]. Although recent advances of NMR in liquids [4] and gas-phase infrared spectroscopy [5 - 7] can identify spectral features of the hydrogen-bonded proton in acetic or formic acid dimers, it is not obvious how a subtle change in geometry by isotope substitution could be deduced.
The carboxylic acid dimer consists of two monomers held together with two hydrogen bondings. Early static ab initio theoretical work was concerned only with the equilibrium geometry and electronic structure [8, 9]. Later, state-of-the-art ab initio quantum chemical calculations are widely performed for acetic and formic acid dimers [10 - 12]. However, it should be noted that nuclear quantum effect is not taken into account in these static calculations.
Recently, some studies including the nuclear quantum effect have been carried out. Shida [13] et al. analyzed the nuclear motion for the double proton transfer of formic acid dimer by reaction surface Hamiltonian using the minimum energy path on the potential energy surface (PES). Kim [10] found that kinetic isotope effect is considerately larger than that observed for single proton transfer from the rate constant calculated with variational transition state theory. Dreyer [14] suggested that anharmonic coupling to low-frequency modes and Fermi resonance coupling with fingerprint modes are essential mechanisms in acetic acid dimer from IR spectrum study with DFT method. On the other hand, Miura [15] et al. have performed the ab initio DFT path integral molecular dynamics simulations for formic acid dimer and found that the large distribution at higher energy regions than the geometrical saddle point due to proton tunneling. Although there are many theoretical researches of acetic acid dimer and formic acid dimer as mentioned above [7 - 15], there are few reports on geometric isotope effect (GIE) of these dimers to our knowledge.
To describe the GIE within the framework of Born-Oppenheimer (BO) [16] approximation, it is necessary to get the PES corresponding to the eigenvalues with the electronic Hamiltonian for all possible nuclear configurations. This treatment is, however, possible to only a few-atoms containing system. On the other hand, we have proposed the multi-component molecular orbital (MC_MO) methods [17 - 21], which can directly include the nuclear quantum effect. We have already analyzed GIEs and the kinetic isotope effects in hydrogen-bonded systems and hydrogen/proton transfer reaction [22 - 24]. Also, several groups developed the similar treatment with the MC_MO method [25 - 28].
In this study, thus, we would focus on the GIE of LBHB systems such as acetic acid dimer, formic acid dimer, and their anionic clusters with the aid of MC_MO method. We have also estimated the interaction energy including the nuclear quantum effect using MC_MO method. We note here that the acetic acid dimer is often picked up as the model of active site of serine protease [29] and HIV-1 aspartyl protease (PR) [30], which catalyzes the synthesis of proteins for the growth of HIV [31]. We believe that our research should be a first step toward theoretical investigation of real and whole system of HIV-1 PR enzyme.

2 Theory

Here, we briefly introduce the MC_MO method, and more detailed information is published elsewhere [17 - 21, 32]. In the MC_MO treatment the concept of molecular orbital is extended to light nuclei such as proton, deuteron, and so on. Since our MC_MO method can directly include the anharmonicity of the potential, we can analyze geometric isotope effect efficiently by using MC_MO method.
The total Hamiltonian used in multi-component theories is given as

where the i and j indices refer to the electrons, p and q to the quantum nuclei, m and n to the classical nuclei, Ne and Np the number of the electrons and quantum nuclei, and M is the number of the classical nuclei. ZA represents the charge of Ath nuclei, and Mp is the mass of the pth quantum nuclei. The first to third terms are the conventional electronic Hamiltonian, and the fourth term represents classical nucleus-nucleus repulsion. The fifth to seventh terms are quantum nuclear terms, which are arisen by extending to multi-component treatment.
In the Hartree-Fock level of multi-component method (MC_HF), the total wavefunction Y0 is written as a simple product of electronic (F0e) and nuclear (F0p) wavefunctions,

The Hartree-Fock equations are derived by the conventional variational principle as

where feHF and fpHF are fock operators for an electron and a nucleus, fi and fp are spatial MOs of an electron and a nucleus, he and hp are one-electronic and one-nuclear operators, Ji, Jp and Ki, Kp are Coulomb and exchange operators for the electrons and the nuclei, respectively. To solve these Fock equations, we introduce the basis set expansion, and both optimum electronic MO fi and quantum nuclear MO fp are obtained by iteratively solving the electronic and quantum nuclear Roothaan equations.
We have also developed the multi-component hybrid density functional theory (MC_(HF+DFT)) [21]. The Kohn-Sham operators used in the MC_(HF+DFT) method are given as,

where VXC and VC are the exchange-correlation and the correlation functionals, (e-e), (e-p), and (p-p) represent the case of electron-electron, electron-nucleus, and nucleus-nucleus, respectively. It has been already reported that the evaluation of dynamical electron-nuclear correlation should be important for the quantitative discussion [19, 26 - 28]. In the present study, however, we only evaluate electron-electron and neglected and in above equations because of the absence of the suitable functionals for (e-p) and (p-p) correlations [18].

3 Computational detail

We have calculated the acetic acid dimer (AD), formic acid dimer (FD), and their anionic clusters (AD- and FD-) using MC_HF, MC_(HF+DFT), conventional HF, and conventional hybrid type DFT methods with 6-31G** electronic basis set. We have chosen BHandHLYP hybrid type functional [33] for both MC_(HF+DFT) and conventional hybrid type DFT calculations.
In the MC_MO calculation several nuclei can be treated as quantum wave function as well as electrons under the field of fixed-nuclear point charges of heavy atoms. In this study only hydrogen-bonded protons/deuterons are treated quantum-mechanically, since we focus on the H/D isotope effect induced by replacing hydrogen-bonded proton(s) to deuteron(s). We here define the notations of isotopomers: (a) HH, non-deuteronic substituted dimer; (b) HD, singly-deuteronic substituted dimer; (c) DD, doubly-deuteronic substituted dimer. The positions of the nuclear point charges are determined by ordinary optimization procedures using analytical gradients [34]. We have employed single s-type Gaussian type function (GTF), exp{-a(r - R)2} , for each protonic and deuteronic basis functions, and optimized two variational parameters (a and R), simultaneously. Although the nuclear wavefunction should be expanded by suitable numberof nuclear basis sets, we have already demonstrated that single s-type GTF is usable to qualitatively discuss the H/D isotope effect in the present study. The centers of electronic GTFs were fixed on each nucleus or quantum nuclear GTF. All calculations were carried out with the modified version of GAUSSIAN 03 program package [35].

4 Result and discussion

4. 1 Optimized geometry of AD, FD, AD-, and FD- by conventional calculation

Tables 1, - 4 list the optimized geometric parameters of AD, FD, AD-, and FD-, respectively. The notations used in these Tables are depicted in Figure 1. For the case of neutral dimers, the difference between ROH and RH...O is larger than 0.6A in both HF and DFT calculations as shown in Tables 1, 2. This means that central hydrogen is clearly one-sided in the optimized neutral dimers. On the other hand, the difference between ROH and RH...O for anionic compounds with DFT method is less than 0.2A as shown in Tables 3, 4, which means central hydrogen settled near the mid-point of two oxygen atoms, while central hydrogen is clearly one-sided in HF-optimized geometry. This is due to the underestimation of the barrier height of the hydrogen-bonding by using the current density functional. We discuss these results in more detail afterward, comparing with the results with MC_BHandHLYP method.

Figure 1. Optimized structure of acetic acid dimer and formic acid dimer

4. 2 Multi-Component calculations

4. 2. 1 Optimized exponent a values and geometry of neutral dimers

Tables 1, 2 also list the optimized geometric parameters and optimized orbital exponent a values for protonic or deuteronic basis functions using MC_HF and MC_BHandHLYP methods for AD and FD, respectively. Because of the difference of mass, the optimized a value for proton is always smaller than that of deuteron in all results of our multi-component calculations. We have carried out some calculation with additional p or d functions, since one way to improve the nuclear wavefunction is to use multiple Gaussians as shown in some previous papers [25 - 28, 36]. However, we have found that there were little difference between the results with single s-type and multiple Gaussians. This is because we have used 1s GTF for the nuclear wavefunction in the present study.

Table 1. Optimized exponent a value and geometric parameters of AD

Table 2. Optimized exponent a value and geometric parameters of FD

First, we focus on the shortening of covalent bond length (ROL) by deuteronic substitution. We note here that the MC_MO method provide the average bond length (<Rave>), but not the equilibrium bond length (Req). Thus, we have obtained the relation of Req < <ROD> < <ROH> for covalent bond length, due to the anharmonicity of the potential directly taken into account in our MC_MO method. We also mention here that the dR, the difference between ROL and RL...O, becomes smaller when the nuclear quantum effect are taken into account in all species. This change indicates that the effective barrier height becomes lower by including protonic motion directly.
The hydrogen-bonded distance (RL...O) and the distance between two oxygen atoms (RO...O) in HH species are shorter than those in HD and DD species. Especially, our results of AD in Table 1 reproduce the overall tendency of the NMR experimental ones. This so-called secondary geometric isotope effect (GIE) can be interpreted with the idea of the Ubbelohde effect [37] in hydrogen-bonded crystal. The explanation of the effect is generally based on the assumption that the hydrogen-bonded proton mediates the attractive interaction of acceptor atoms more than the hydrogen-bonded deuteron. Since protonic wavefunction is more delocalized than deuteronic one, the hydrogen-bonded proton strongly attracts acceptor atoms more than the deuteron. This idea can be used for the geometrical change of RL...O and RO...O in our calculation, because we have actually observed that the exponent a value of optimum protonic wavefunction is smaller than that of deuteronic one, as mentioned above. Strictly speaking, our calculated GIE is still overestimated by comparison of the corresponding NMR results in Table 1, though the overall qualitative tendency is reproduced. For more quantitative reproduction, dynamical electron-nuclear correlation would be required to include the low-frequency vibrational motions [36, 38].

4. 2. 2 Optimized exponent a value and geometry for anionic compounds

The optimized geometrical parameters and exponent a values for the anionic species, AD- and FD-, are shown in Tables 3, 4. The tendency of geometrical changes induced by H/D isotope effect in anionic species with MC_HF results is similar to that in neutral dimers. In the case of MC_BHandHLYP method, the difference between ROL and RL...O is less than 0.01A, that is, the central hydrogen settles in a mid-point of heavy atoms for both protonic and deuteronic species.

Table 3. Optimized exponent a value and geometric parameters of AD-

Table 4. Optimized exponent a value and geometric parameters of FD-

First, we discuss the optimized a values for anionic compounds. The a value is around 18~19 and 27~28 for anionic H and D species, while 21~22 and 31~32 for neutral HH and DD species, respectively. These values for anionic compounds seem to be significantly smaller than those for neutral compounds and those in the previous work [17 - 21]. Next, in order to interpret the DFT results, we compare the MC_BHandHLYP results with conventional BHandHLYP ones as shown in Tables 1, - - 4. As discussed above, the difference between ROL and RL...O with conventional BHandHLYP is obviously smaller than that with conventional HF method, due to the underestimation of the barrier height [39] with the current exchange-correlation functionals. Consequently, central proton can exceed the underestimated barrier in LBHB by including nuclear quantum effect, and the hydrogen atom is situated near the central position of two oxygen atoms by MC_BHandHLYP method. This geometry, however, is not in agreement with the corresponding experimental results even qualitatively. Experimental results indicate that the shape of potential energy curve is double well with low barrier. Since the lowering of the effective barrier height by including the nuclear quantum effect is the reasonable change and the underestimation of barrier height with conventional DFT functional is a well known problem, we consider this inadequate geometrical feature is due to the exchange-correlation functional. In order to theoretically analyze these anionic dimers with the MC_DFT method, new exchange-correlation functionals, which can adequately describe the reaction barrier, would be required.

4. 3 Interaction energy

We also estimate the interaction energy (DE) of these compounds to analyze the nuclear quantum effect on the interaction energy. Table 5 lists the calculated interaction energy DE with HF, BHandHLYP, MC_HF, and MC_BHandHLYP methods. In the case of neutral species, we defined the DE as

where Edimer_N and Emonomer_N are the energies of neutral dimer and neutral monomer, respectively. For anionic compounds, we defined the DE as

where Edimer_A and Emonomer_A are the energies of anionic dimer and anionic monomer, respectively.

Table 5. Calculated DEa obtained with MC_HF and MC_(HF+DFT) Calculation.

aUnits are in kcal/mol

For all compounds, the absolute value of DE decreases as the hydrogen is substituted by the deuterium, which means the interaction between the two monomers is weakening by deuteronic substitution. The interaction energy differences between H and D species (DDEH-D = DEH - DED) of both anionic species are 1.6 and 1.4 kcal/mol with MC_HF for MC_BHandHLYP calculations, respectively. It should be noted that the interaction energy differences between H and D species DDEH-D are comparable to that with MC_HF, although the GIE is very small in anionic compounds with MC_BHandHLYP as shown in Tables 3, 4.
To analyze this result from a viewpoint of the distribution of nuclear wave functions, we have calculated the energies of AD- with following two conditions: (a) using optimum protonic exponent values (aH) with optimized geometry of H species for monomer and dimer (EHmonomer_N(aH), EHdimer_A(aH)); (b) using optimum deuteronic exponent values (aD) with optimized geometry of D species for monomer and dimer (EHmonomer_N(aD), EHdimer_A(aD)). Here, we classify the contribution to DDEH-D into that of the distribution of wave function and of geometrical change. Additionally, we define DE(aH) and DE(aD) as below,

The energy differences between DEH(aH) and DEH(aD), defined as DDE, are 1.7 and 1.4 kcal/mol for MC_HF and MC_BHandHLYP calculations, respectively. This DDE contains no contribution of geometrical change, but that of the distribution of nuclear wave functions to the isotope effect on interaction energy. We note here that the DDE value is almost the same as DDEH-D, which means that the contribution from difference of the distribution of wavefunction accounts for almost all of the DDEH-D. This result corresponds to our anticipation, because the interaction energy of H species is larger than that of D species even the hydrogen-bonded distance of D species is shorter than that of H species. This result suggests to us that quantum-mechanical treatment is required for quantitative analysis of interaction energy because the optimum exponent a value is easily varied by the effect of surroundings even in the same system, for example, the a values are different in stationary point and transition state.

5 Conclusion

We have analyzed the geometric isotope effect (GIE) of acetic acid dimer (AD), formic acid dimer (FD), and their anionic clusters by using MC_MO method with HF level (MC_HF method) and hybrid DFT level (MC_BHandHLYP) of calculation, which can directly take account of the nuclear quantum effect as well as electrons.
  1. The optimized geometries of the neutral dimers of AD and FD with both MC_HF and MC_BHandHLYP methods reproduce the overall qualitative tendency of the corresponding experimental GIE.
  2. On the other hand, the optimized geometries of anionic compounds of AD- and FD- with MC_BHandHLYP are in poor agreement with the experiment ones, while those with MC_HF method reproduce the overall qualitative tendency. This is due to the poor description of the nuclear wavefunction and barrier height of the low barrier hydrogen-bonding system with conventional exchange-correlation functional such as BHandHLYP.
  3. Our study clearly demonstrates that the interaction energy of hydrogen-bonded systems is sensitive to the distribution of nuclear wavefunction. DE from the change of distribution of nuclear wavefunction DE (1.7 kcal/mol for MC_HF and 1.4 kcal/mol for MC_BHandHLYP) covers the most part of whole DDEH-D (1.6 kcal/mol for MC_HF and 1.4 kcal/mol for MC_BHandHLYP), which includes the effect of geometrical relaxation.
For further analysis of LBHB systems with MC_(HF+DFT) method, we may need to develop the new density functional that can adequaely represent the barrier height of hydrogen-bonded system. To overcome a difficulty of description of the barrier height, several methods have already been proposed such as long-range correction scheme [40]. When we analyze hydrogen-bonded systems with such higher-accuracy method that can include e-p correlation in the near future, multiple Gaussians will be necessary for the quantitative description of the nuclear wavefunction.
In order to analyze the active site of HIV-1 PR in greater detail at the next stage, we need to treat not only for the model of the active space, but also including larger or whole structure of HIV-1 PR. Our MC_MO method has already been extended to the fragment molecular orbital (FMO) method [41], which can treat biological macromolecules such as proteins. By combining our MC_(HF+DFT) method with FMO-MC_MO method [42] and long-range correction scheme, the quantitative analysis of the larger model of the active site with including the nuclear quantum nature will be brought to realization.

This work is supported by Grant-in-Aid for Scientific Research and for the priority area by Ministry of Education, Culture, Sports, Science and Technology, Japan.


[ 1] (a) R. Blinc, J. Phys. Chem. Solids, 13, 204 (1960).
(b) M. I. McMahon, R. J. Nelmes, W. F. Kuhs, R. Derwarth, R. O. Plitz, Z. Tun, Nature, 348, 317 (1990).
(c) M. Tachikawa, T. Ishimoto, H. Tokiwa, H. Kasatani, K. Deguchi, Ferroelectrics, 268, 3 (2002).
[ 2] G. A. Jeffrey, An Introduction to Hydrogen Bonding, Oxford University Press, Oxford, U. K. (1997), 58.
[ 3] J. Mcmurry, Organic Chemistry, Cole Pub Co., Brooks (2004).
[ 4] P. M. Tolstoy, P. S-Mohammedi, S. N. Smirnov, N. S. Golubev, G. S. Denisov, H. H. Limbach, J. Am. Chem. Soc., 126, 5621 (2004).
[ 5] J. B. Togeas, J. Phys. Chem. A, 109, 5438 (2005).
[ 6] F. Madeja, M. Havenith, J. Chem. Phys., 117, 7162 (2002).
[ 7] K. Heyne, N. Huse, J. Dreyer, E. T. J. Nibbering, T. Elsaesse, S. Mukamel, J. Chem. Phys., 121, 902 (2004).
[ 8] S. Hayashi, J. Umemura, S. Kato, K. Morokuma, J. Phys. Chem., 88, 1330 (1984).
[ 9] Y. T. Chang, Y. Yamaguchi, Y. H. Miller, H. F. Schaefer III, J. Am. Chem. Soc., 109, 7245 (1987), , and references therein..
[10] Y. Kim, J. Am. Chem. Soc., 118, 1522 (1996).
[11] R. W. Gora, S. J. Grabowski, J. Leszczynski, J. Phys. Chem. A, 109, 6397 (2005).
[12] S. J. Grabowski, J. Phys. Org. Chem., 21, 694 (2008).
[13] N. Shida, P. F. Barbara, J. Almlof, J. Chem. Phys., 94, 3633 (1991).
[14] J. Dreyer, J. Chem. Phys., 122, 184406 (2005).
[15] S. Miura, M. E. Tuckerman, M. L. Klein, J. Chem. Phys., 109, 5290 (1998).
[16] M. Born, J. R. Oppenheimer, Ann. Phys., 84, 457 (1927).
[17] M. Tachikawa, K. Mori, K. Suzuki, K. Iguchi, Int. J. Quantum Chem., 70, 491 (1998).
[18] M. Tachikawa, K. Mori, H. Nakai, K. Iguchi, Chem. Phys. Lett., 290, 437 (1998).
[19] M. Tachikawa, Chem. Phys. Lett., 360, 494 (2002).
[20] T. Ishimoto, M. Tachikawa, U. Nagashima, J. Chem. Phys., 125, 144103 (2006).
[21] T. Udagawa, M. Tachikawa, J. Chem. Phys., 125, 244105 (2006).
[22] T. Udagawa, T. Ishimoto, H. Tokiwa, M. Tachikawa, U. Nagashima, J. Phys. Chem. A, 110, 7279 (2006).
[23] T. Ishimoto, M. Tachikawa, H. Tokiwa, U. Nagashima, Chem. Phys., 314, 231 (2005).
[24] Y. Itou, S. Mori, T. Udagawa, M. Tachikawa, T. Ishimoto, U. Nagashima, J. Phys. Chem. A, 111, 261 (2007).
[25] H. Nakai, K. Sodeyama, M. Hoshino, Chem. Phys. Lett., 118, 345 (2001).
[26] M. Cafiero, B. Bubin, L. Adamowicz, Phys. Chem. Chem. Phys., 5, 1491 (2003).
[27] C. Swalina, M. V. Pak, A. Chakraborty, S. Hammes-Schiffer, J. Phys. Chem. A, 110, 9983 (2006).
[28] A. Chakraborty, M. V. Pak, S. Hammes-Schiffer, J. Chem. Phys., 129, 014101 (2008).
[29] P. A. Frey, Isotope Effects in Chemistry and Biology, C.H.I.P.S., Texas (2006), 975.
[30] R. Smith, I. M. Brereton, R. Y. Chai, S. B. H. Kent, Nature Struc. Biol., 3, 946 (1996).
[31] T. J. Mcquade, A. G. Tomasselli, L. Liu, V. Karacostas, B. Moss, T. K. Sawyer, R. L. Heinrikson, W. G. Tarpley, Science, 247, 454 (1990).
[32] T. Udagawa, M. Tachikawa, Progress in Quantum Chemistry Research, NOVA Science Publishers, New York (2007), 123.
[33] A. D. Becke, J. Chem. Phys., 98, 1372 (1993).
[34] P. Pulay, Applications in Electronic Structure Theory, Plenum Press, New York (1977).
[35] M. J. Frisch, G. W. Trucks, H. B. Schlegel, et al., GAUSSIAN03, Revision B.05, Gaussian, Inc., Pittsburgh, PA (2003).
[36] C. Swalina, S. Hammes-Schiffer, J. Phys. Chem. A, 109, 10410 (2005).
[37] A. R. Ubbelohde, K. J. Gallagher, Acta Crystallographica, 8, 71 (1995).
[38] P. O. Astrand, G. Karlstrom, A. Engdahl, B. J. Nelander, Chem. Phys., 102, 3534 (1995).
[39] J. M. Perez-Jorda, A. D. Becke, Chem. Phys. Lett., 233, 134 (1995).
[40] H. Iikura, T. Tsuneda, T. Yanai, K. Hirao, J. Chem. Phys., 115, 3540 (2001).
[41] K. Kitaura, E. Ikeo, T. Asada, T. Nakano, M. Uebayashi, Chem. Phys. Lett., 313, 701 (1999).
[42] T. Ishimoto, M. Tachikawa, U. Nagashima, J. Chem. Phys., 124, 014112 (2006).