Geometric Isotope Effect on Low Barrier HydrogenBonding Systems of Acetic Acid Dimer, Formic Acid Dimer, and Their Anionic Clusters by Using the MultiComponent Molecular Orbital Method
Masato KANEKO, Taro UDAGAWA and Masanori TACHIKAWA
Return
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 hydrogenbonded ferroelectric crystal systems. For example, it is believed that large isotope effect in the phase transition temperature of hydrogenbonded 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 hydrogenbonded molecular clusters such as low barrier hydrogenbonded (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 gasphase infrared spectroscopy [5  7] can identify spectral features of the hydrogenbonded 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, stateoftheart 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 lowfrequency 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 BornOppenheimer (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 fewatoms containing system. On the other hand, we have proposed the multicomponent 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 hydrogenbonded 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 HIV1 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 HIV1 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 multicomponent 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, N_{e} and N_{p} the number of the electrons and quantum nuclei, and M is the number of the classical nuclei. Z_{A} represents the charge of Ath nuclei, and M_{p} is the mass of the pth quantum nuclei. The first to third terms are the conventional electronic Hamiltonian, and the fourth term represents classical nucleusnucleus repulsion. The fifth to seventh terms are quantum nuclear terms, which are arisen by extending to multicomponent treatment.
In the HartreeFock level of multicomponent method (MC_HF), the total wavefunction Y_{0} is written as a simple product of electronic (F_{0}^{e}) and nuclear (F_{0}^{p}) wavefunctions,
The HartreeFock equations are derived by the conventional variational principle as
where f_{e}^{HF} and f_{p}^{HF} are fock operators for an electron and a nucleus, f_{i} and f_{p} are spatial MOs of an electron and a nucleus, h_{e} and h_{p} are oneelectronic and onenuclear operators, J_{i}, J_{p} and K_{i}, K_{p} 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 f_{i} and quantum nuclear MO f_{p} are obtained by iteratively solving the electronic and quantum nuclear Roothaan equations.
We have also developed the multicomponent hybrid density functional theory (MC_(HF+DFT)) [21]. The KohnSham operators used in the MC_(HF+DFT) method are given as,
where V_{XC} and V_{C} are the exchangecorrelation and the correlation functionals, (ee), (ep), and (pp) represent the case of electronelectron, electronnucleus, and nucleusnucleus, respectively. It has been already reported that the evaluation of dynamical electronnuclear correlation should be important for the quantitative discussion [19, 26  28]. In the present study, however, we only evaluate electronelectron
and neglected
and
in above equations because of the absence of the suitable functionals for (ep) and (pp) 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 631G** 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 fixednuclear point charges of heavy atoms. In this study only hydrogenbonded protons/deuterons are treated quantummechanically, since we focus on the H/D isotope effect induced by replacing hydrogenbonded proton(s) to deuteron(s). We here define the notations of isotopomers: (a) HH, nondeuteronic substituted dimer; (b) HD, singlydeuteronic substituted dimer; (c) DD, doublydeuteronic substituted dimer. The positions of the nuclear point charges are determined by ordinary optimization procedures using analytical gradients [34]. We have employed single stype 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 stype 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 R_{OH} and R_{H...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 onesided in the optimized neutral dimers. On the other hand, the difference between R_{OH} and R_{H...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 midpoint of two oxygen atoms, while central hydrogen is clearly onesided in HFoptimized geometry. This is due to the underestimation of the barrier height of the hydrogenbonding 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 MultiComponent 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 multicomponent 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 stype 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 (R_{OL}) by deuteronic substitution. We note here that the MC_MO method provide the average bond length (<R_{ave}>), but not the equilibrium bond length (R_{eq}). Thus, we have obtained the relation of R_{eq} < <R_{OD}> < <R_{OH}> 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 R_{OL} and R_{L...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 hydrogenbonded distance (R_{L...O}) and the distance between two oxygen atoms (R_{O...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 socalled secondary geometric isotope effect (GIE) can be interpreted with the idea of the Ubbelohde effect [37] in hydrogenbonded crystal. The explanation of the effect is generally based on the assumption that the hydrogenbonded proton mediates the attractive interaction of acceptor atoms more than the hydrogenbonded deuteron. Since protonic wavefunction is more delocalized than deuteronic one, the hydrogenbonded proton strongly attracts acceptor atoms more than the deuteron. This idea can be used for the geometrical change of R_{L...O} and R_{O...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 electronnuclear correlation would be required to include the lowfrequency 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 R_{OL} and R_{L...O} is less than 0.01A, that is, the central hydrogen settles in a midpoint 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 R_{OL} and R_{L...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 exchangecorrelation 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 exchangecorrelation functional. In order to theoretically analyze these anionic dimers with the MC_DFT method, new exchangecorrelation 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 E_{dimer_N} and E_{monomer_N} are the energies of neutral dimer and neutral monomer, respectively. For anionic compounds, we defined the DE as
where E_{dimer_A} and E_{monomer_A} are the energies of anionic dimer and anionic monomer, respectively.
Table 5. Calculated DE^{a} obtained with MC_HF and MC_(HF+DFT) Calculation.
^{a}Units 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 (DDE^{HD} = DE^{H}  DE^{D}) 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 DDE^{HD} 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 (a_{H}) with optimized geometry of H species for monomer and dimer (E^{H}_{monomer_N}(a_{H}), E^{H}_{dimer_A}(a_{H})); (b) using optimum deuteronic exponent values (a_{D}) with optimized geometry of D species for monomer and dimer (E^{H}_{monomer_N}(a_{D}), E^{H}_{dimer_A}(aD)). Here, we classify the contribution to DDE^{HD} into that of the distribution of wave function and of geometrical change. Additionally, we define DE(a_{H}) and DE(a_{D}) as below,
The energy differences between DE^{H}(aH) and DE^{H}(aD), defined as DDE^{Hα}, are 1.7 and 1.4 kcal/mol for MC_HF and MC_BHandHLYP calculations, respectively. This DDE^{Hα} 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^{Hα} value is almost the same as DDE^{HD}, which means that the contribution from difference of the distribution of wavefunction accounts for almost all of the DDE^{HD}. This result corresponds to our anticipation, because the interaction energy of H species is larger than that of D species even the hydrogenbonded distance of D species is shorter than that of H species. This result suggests to us that quantummechanical 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.

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.

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 hydrogenbonding system with conventional exchangecorrelation functional such as BHandHLYP.

Our study clearly demonstrates that the interaction energy of hydrogenbonded systems is sensitive to the distribution of nuclear wavefunction. DE from the change of distribution of nuclear wavefunction DE^{Hα} (1.7 kcal/mol for MC_HF and 1.4 kcal/mol for MC_BHandHLYP) covers the most part of whole DDE^{HD} (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 hydrogenbonded system. To overcome a difficulty of description of the barrier height, several methods have already been proposed such as longrange correction scheme [40]. When we analyze hydrogenbonded systems with such higheraccuracy method that can include ep 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 HIV1 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 HIV1 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 FMOMC_MO method [42] and longrange 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 GrantinAid for Scientific Research and for the priority area by Ministry of Education, Culture, Sports, Science and Technology, Japan.
References
[ 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. SMohammedi, 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. HammesSchiffer, J. Phys. Chem. A, 110, 9983 (2006).
[28] A. Chakraborty, M. V. Pak, S. HammesSchiffer, 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. HammesSchiffer, 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. PerezJorda, 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).
Return