- Email: [email protected]
Car–Parrinello simulation of the vibrational spectra of strong hydrogen bonds with isotopic substitution effects: Application to oxalic acid dihydrate
Chemical Physics Letters 558 (2013) 88–92
Contents lists available at SciVerse ScienceDirect
Chemical Physics Letters journal homepage: www.elsevier.com/locate/cplett
Car–Parrinello simulation of the vibrational spectra of strong hydrogen bonds with isotopic substitution effects: Application to oxalic acid dihydrate Mateusz Z. Brela a,⇑, Marek J. Wójcik b,⇑, Marek Boczar b, Rauzah Hashim c a
Research Group of Molecular Modeling of Catalytic Processes, Faculty of Chemistry, Jagiellonian University, Ingardena 3, 30-060 Krakow, Poland Laboratory of Molecular Spectroscopy, Faculty of Chemistry, Jagiellonian University, Ingardena 3, 30-060 Krakow, Poland c Department of Chemistry, University of Malaya, 50603 Kuala Lumpur, Malaysia b
a r t i c l e
i n f o
Article history: Received 12 October 2012 In final form 18 December 2012 Available online 3 January 2013
a b s t r a c t The nature of strong intermolecular hydrogen bonding in oxalic acid dihydrate in the crystal phase was examined by infrared spectroscopy and Car–Parrinello molecular dynamics simulation. We studied region of infrared spectra associated with the O—H modes. The spectra were calculated using harmonic approximation with crystal field and time course of the dipole moment as obtained from Car–Parrinello simulation with quantization of the O—H motion, and isotopic substitution. We obtained good agreement of the molecular dynamic simulation with experiment. To our best knowledge, this is one of the first Car– Parrinello calculations of infrared spectra including anharmonicity effects and crystal field interactions. Ó 2012 Elsevier B.V. All rights reserved.
1. Introduction Hydrogen bond is one of the most important interactions in chemical and biochemical systems. For example in Nature, hydrogen bonds are responsible for the formations and functions of many organized complex structures such as DNA [1–3] biopolymers [4–6] and cell membranes [7,8]. The properties of hydrogen bonds have been widely investigated theoretically and experimentally and have been subject of several monographs [9–13]. Many theoretical models have been proposed to describe unusual features of hydrogen bond stretching bands in the IR spectra [14–23]. Oxalic acid is produced naturally from plants and has been shown to play an important role in the prevention of chronic diseases (osteoporosis, obesity), but its biological function in human body, almost unknown, is beginning to be investigated [24]. Infrared spectra of oxalic acid have been subject of several theoretical and experimental investigations, most of them focused on hydrogen bond stretching band [25–30]. Oxalic acid dihydrate forms at room temperature white crystals with melting point at 394 K. The crystal structure of oxalic acid dihydrate was determined from neutron diffraction by Sabine et al. [31] and from X-ray by Leiserowitz [32] and Coppens [33]. The crystal of oxalic acid dihydrate belongs to the P21/n space group of symmetry. The unit cell of oxalic acid dihydrate contains two molecules of oxalic acid and four molecules of water. The crystal structure is supported by two types of hydrogen bonds (Fig⇑ Corresponding authors. E-mail address: [email protected] (M.J. Wójcik). 0009-2614/$ - see front matter Ó 2012 Elsevier B.V. All rights reserved. http://dx.doi.org/10.1016/j.cplett.2012.12.035
ure 1A). One type connects the water molecules as donors and the acid carbonyl groups as acceptors; it may be classified as normal for its length in the range of 2.869 Å. The other type connects the acid’s hydroxyl group as the donor to water oxygen as acceptor. This bond is very short; at room temperature the O O distance is 2.499 Å. The OH bond length is 1.041 Å at 298 K and it expands to 1.071 Å at 80 K. The recently investigated single crystal Raman spectra of oxalic acid dihydrate and its analog [29] have shown interesting features thus rising interest in a more detailed study of IR spectra and their theoretical modelling. The latter is of even broader interest because of the rarity of theoretical modelling of the spectra of examples of hydrogen bonds of O O distances about 2.50 Å. In this Letter we present results of ab initio Car–Parrinello molecular dynamics (CPMD) calculations [34] performed for the unit cell of the crystal of oxalic acid dihydrate. The vibrational spectrum of the oxalic acid dihydrate has been calculated, analyzed, and compared with experimental spectra. The greatest emphasis has been put on the reconstruction of the hydrogen bond vibrational bands. The mO–H band has been thoroughly examined, and the simulation of the mO–H bandshape, based on CPMD calculations, has been performed. In the first step, we optimized the geometry of the unit cell of oxalic acid dihydrate, then we calculated trajectories using the Car–Parrinello molecular dynamics (CPMD) scheme. CPMD calculations have recently become a popular tool for the interpretation of infrared spectra of hydrogen-bonded systems [35–41]. In the present Letter we took into account the periodicity of the crystal and thermal fluctuations. The second step involved quantization
89
M.Z. Brela et al. / Chemical Physics Letters 558 (2013) 88–92
Figure 1. Structure and atom labeling in crystal unit cell of oxalic acid dihydrate. The atoms are color coded: carbon = gray, oxygen = red, hydrogen = white. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)
of the proton motion, for which we used an alternative method based on an a posteriori quantization of the nuclear degrees of freedom associated with proton dynamics. The vibrational spectrum of the oxalic acid dihydrate crystal has been calculated, analyzed, and compared with experimental data and also with the results of DFT calculations performed for a crystal cell of the oxalic acid dihydrate. The largest emphasis has been placed on the reconstruction of the stretching bands of hydrogen bonded O—H groups. The tO–H bands, obtained from CPMD calculations and instantaneous one-dimensional snapshot potentials, extracted from the trajectory calculated using the Wannier functions [42] were compared with the experimental bands of the oxalic acid dihydrate crystal. We discuss also the isotopic effects on the spectra.
then calculated using the proton potentials. For each of the extracted snapshot structures, a one-dimensional (1D) proton potential function was obtained by stepwise displacement of one hydrogen-bonded proton along the hydrogen vector parallel to the pertinent O—H line. The corresponding internal coordinate (x) was defined as the distance between hydrogen and oxygen atoms. The scan typically covered the x range from 0.5 Å to 0.5 Å in relation to the equilibrium O—H distance. The hydrogen was displaced along the distance steps of 3.33% of the range. It is worth to point out that the proton potentials are ensemble averages and are hence mass independent in the limit of classical statistical mechanics. Therefore the snapshot methodology for deuterated species is applicable on the trajectory obtained by integrating equations of motion for hydrogen form.
2. Computational methods
3. Results and discussion
We considered several isotopically substituted molecules. Geometry optimizations and vibrational frequency calculations were performed to characterize the key stationary points on the potential energy surface (PES) for unit cell of oxalic acid dihydrate (Figure 1B). We have also performed a vibrational analysis in the optimized crystal structure using the (CPMD) program package version 3.13.2 [43]. The BLYP [44] density functional together with a plane-wave basis set with a kinetic energy cutoff of 300 Ry were used in our study. The Goedecker atomic pseudopotentials [45] were used for the treatment of the core electrons. The starting point for our molecular dynamics (MD) simulation was the neutron structure data obtained by Sabine et al. [31]. Three-dimensional periodicity was fully implemented in our calculation. The MD simulation was carried out in the unit cell fixed to the experimentally determined size and shape and using the Car– Parrinello scheme [34]. The fictitious orbital mass was set to 300 a.u, the temperature was 300 K and the molecular dynamics time step was set to 1 a.u. (about 0.024188 fs), while other properties were set the same as the static calculation. The total simulation time was about 10 ps. The simulation was performed on an HP Cluster Platform 3000 BL 2x220 PC/LINUX and it lasted for about 720 days using only 12 cores. The anharmonic vibrational spectrum was calculated using the Fourier transformation of the dipole autocorrelation function obtained from dipole trajectories generated by the CPMD simulation facilitated by the scripts of Kohlmeyer [46]. The dipole dynamics were calculated using the optimally localized Wannier functions [42]. Snapshot structures of the trajectory were extracted every 4000 steps (about 97 fs), yielding 125 distinct structures, representative of the fluctuating crystalline environment. Note that due to the fact that two molecules which interacted in four hydrogen bonds in the unit cell, the total number of distinct hydrogen bonded moieties extracted was 500. The contour of the O—H stretching band was
3.1. Geometry analysis The results of our study are presented in the following order: analysis of the structure, and vibrational analysis. In Figure 1A we have shown the structure and atom labeling of oxalic acid dihydrate. The geometry results of static calculations are shown in Table 1. Analysis of the Car–Parrinello trajectory for hydrogen bond formed by O–H group of oxalic acid and water shows a time-course of the OO and O—H distances in the O—H O hydrogen bond and shows that the proton is asymmetrically located nearer to the donor oxygen. The average O O distance is 2.51 Å. This suggests that the hydrogen bonds are strong. Over the entire CPMD run the O O distance spans the range between 2.44 and 2.58 Å. It is slightly longer than the experimental crystallographic value of 2.52 Å. The diffraction experiment [33] gave an O—H distance of 1.03 Å. Our DFTs calculations predicted that this bond is longer, with an average bond length 1.05 Å. Generally our CPMD static calculations reproduce the geometric parameters for hydrogen bonds in the crystal reasonably well. 3.2. Infrared spectra of crystal structure calculated in harmonic approximation and autocorrelation function obtained after CPMD Figure 2 presents crystal spectrum calculated in harmonic approximation. The experimental infrared spectrum of oxalic acid
Table 1 The calculated hydrogen bond geometry for four hydrogen bonds in the crystal unit cell. The DFT calculations were performed in the crystal field. See Figure 1A for symbol description.
O—H O O
I-HB
II-HB
III-HB
IV-HB
Average
Exp. [33]
1.01 2.69
1.09 2.50
1.09 2.45
1.01 2.56
1.05 2.55
1.05 2.52
90
M.Z. Brela et al. / Chemical Physics Letters 558 (2013) 88–92
Figure 2. Infrared spectrum of oxalic acid dihydrate crystal calculated in harmonic approximation. Insert presents fragment of crystal structure of the unit cell. The atoms are color coded: carbon = gray, oxygen = red, hydrogen = white. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)
dihydrate was presented in the paper Mohacˇek-Grošev et al. [29]. In the static calculation we obtained peaks, which represent stretching modes of OH groups. Our calculations shows that they represent hydrogen-bonded O–H stretching in oxalic acid (2300 cm 1), weakly hydrogen-bonded O–H stretching in water (2900 and 3100 cm 1) and O–H stretching of free groups in water (3300 and 3450 cm 1). The band at 2300 cm 1 represents the mixing of two kinds of modes – stretching and twisting carboxyl
Figure 3. Infrared spectra of hydrogen form of oxalic acid dihydrate. Inset presents fragment of crystal structure of the unit cell. The atoms are color coded: carbon = gray, oxygen = red, hydrogen = white. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)
group. Harmonic calculations cannot reproduce the broadening and intensity increase caused by the anharmonic intermolecular couplings present in hydrogen bonds [14–16]. In harmonic quantum chemical calculations one does not obtain bandwidths and vibrational transitions are represented by d functions. The bandshapes presented in Figure 2 were obtained by superimposing a Gaussian function with a half-width 50 cm 1 on each calculated transition; such representations have often been used, e.g. in articles [21–23].
Figure 4. Infrared spectra of the four kind isotopic substitution of oxalic acid dihydrate with crystal structure of the unit cell. The atoms are color coded: carbon = gray, oxygen = red, hydrogen = white, deuterium = green. For details see text in Section 3.4 ‘‘Isotope effects’’. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)
91
M.Z. Brela et al. / Chemical Physics Letters 558 (2013) 88–92
Table 2 Comparison between frequencies of selected vibrations in oxalic acid dihydrate crystal obtained from CPMD static calculations (harmonic frequencies obtained using BLYP in crystal field) and experimental data from paper Mohacˇek-Grošev et al. [29].
3.3. Snapshot potential energy functions We calculated one-dimensional (1D) proton potentials corresponding to the instantaneous snapshot structures extracted from the CPMD trajectory. It should be noted that typical 1D scan requires 30 points at which the energy is evaluated. We compared the bands calculated with the snapshot methodology [47,48] and band extracted from the powder spectrum obtained using Fourier transformation. We generated one-dimensional proton potentials using one linear pathway. Figure 5 shows schematically the proton potential construction and obtained potentials together with the solution of the vibrational Schrödinger equation. More than 90% of frequencies, which were calculated by using the three types of 1D potentials, span the range from 1750 to 2250 cm 1. The contour of the band is similar to experimental spectrum. It shows the small anharmonicity and structure of modes. In general, snapshot proton potentials constructed following the O—H direction are practically not shifted. This strategy also yields a narrower band contour with a half-width of about 230 cm 1. 3.4. Isotope effects We have also performed the analysis of spectra for oxalic acid dihydrate and its selected deuterated isotopomers. Figure 4 shows
Experimental frequencies (cm 1) [29]
Description
115 1175 1177 1193 1217
154R 1133
m (O O) c (O–H O)
1289 1318 1338 1345
1253
d (C–O–H)
m (C@O)
1430 1440 1504 1513
Figure 5. O—H stretching band contours of oxalic acid dihydrate: Black line contour was calculated from individual fundamental vibrational transitions as superposition of Gaussian functions with a half width of 10 cm 1. Gray d functions represent fundamental vibrational transitions. The red contour represents the spectrum calculated by Fourier transform of the autocorrelation function of the atoms position obtained from Car–Parrinello trajectory. The bands were normalized. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)
The second methodology, which we used to obtain these spectra, was the Fourier transform of the autocorrelation function. We performed calculations with two types of functions: the total dipole moment (Figures 3 and 4) and function describing atom positions to reproduce data for the analysis of a single band O–H (Figure 5). This approximation gives qualitative reproduction of the experimental spectra, though their positions are shifted to higher frequencies for all O–H bands. The region of high frequency is composed of three peaks, centered at 1750–2100 cm 1, 3100– 3350 cm 1 and small peak near 3600–3700 cm 1. Our calculations show that, they represent hydrogen bonded O–H groups in oxalic acid, weakly hydrogen-bonded O–H groups in water, and free O– H groups in water. The agreement with the experimental broad band with a maximum at 1900 cm 1 [29] which represents the O–H stretching band (strongly hydrogen-bonded), is reasonably good. The frequencies of selected vibration are shown in Table 2.
Calculated frequencies (cm 1)
m
1570 1588 1600 1616
1620
d (O H)water
1716 1722 1734 1758
1690
m (C@O)
2358 2454 2496 2676
1870 2009
m (O H)
3343 3346 3388 3390 3639 3688 3692 3721
3431 3489
m (O H)water
Stretching, d
in-plane bending, c
out-of plane bending, R – Raman.
the infrared spectra obtained using the autocorrelation functions of dipole moment (CPMD simulation) for selected deuterated isotopomers of oxalic acid dihydrate. For all these structures we performed independent CPMD simulations with the same parameters (see the computational part of this article). Analysis of calculated spectra allows us to study the interaction involving O—H stretching vibrations in the crystal. Figure 4A shows the infrared spectrum of the crystal with only one deuterium replacing hydrogen. We observed the stretching mode of O–D near 1700– 1850 cm 1. Additional effect of this substitution is increase in intensity of the band near 3600 cm 1 relative to the band near 3300 cm 1. The second structure was obtained by replacing all hydrogens except one in oxalic acid by deuterium atoms. Figure 4B. shows the result of this substitution. We obtained broad bands near 1600–1900 and 2200–2400 cm 1. The first one represents O–D vibrations of oxalic acid and the second one is composed of O—H vibrations of oxalic acid and O–D vibrations of water molecules. Figure 4C shows the spectrum of the crystal with deuterated oxalic acid and non-deuterated water molecules. In this spectrum the band 2200–2400 cm 1 of non-deuterated oxalic acid dihydrate (Figure 3) is shifted to 1600–1800 cm 1. The O—H stretching vibrations of water in this form have frequencies close to O—H frequencies of non-deuterated oxalic acid dihydrate. Figure 4D shows the spectra calculated for per-deuterated oxalic acid dihydrate. The O—D stretching bands for oxalic acid are similar as in Figure 4C but their center is shifted to slightly higher frequency. The O—D stretching modes of water molecule are observed near 2400 and 2700 cm 1. This result shows that the coupling of O—H vibrations
92
M.Z. Brela et al. / Chemical Physics Letters 558 (2013) 88–92
is observed only in non-deuterated oxalic acid dihydrate. The coupling occurs between O–H groups present in oxalic acid or water molecules and there is no significant coupling between OH groups of oxalic acid and water. This can be explained by difference of OH vibrational frequencies in water and oxalic acid. Similar effect was obtained by Latajka and Durlak in the gas phase [39,49]. This kind of phenomena gives perspectives for QM/MM studies of nuclear quantum effects for example in enzyme centers [50,51]. 4. Conclusions In our Letter we present results of the simulation of spectra of the crystal of oxalic acid dihydrate. The results can be summarized as follows: 1. Infrared spectrum calculated by CPMD method reproduces reasonably well the high-frequency bands in the experimental spectrum. Also the deuterium shifts of vibrational bands were calculated in good accordance with experimental data due to the application of dipole momentum dynamics. 2. The structure of the O–H stretching band is reproduced fairly well. In this part of our Letter, we investigated the snapshot methodology. Calculated bands describing hydrogen bonds vibrations better reproduce experimental values than results of static ab initio calculations at BLYP in the crystal field level performed for the crystal cell. 3. The identification of bands in calculated infrared spectrum is difficult in FT calculations but snapshots methodology can be useful in description of single bands. 4. Calculations for isotopically substituted molecules allow to study couplings effects and intermolecular interactions in crystal of oxalic acid dihydrate. We focused our research only on stretching bands. We observe couplings of similar modes. 5. Improvement in calculation results proves that CPMD method is adequate for investigations of systems with hydrogen bonds as it takes into account most of mechanisms determining the hydrogen bonds dynamics (anharmonicity, couplings between vibrational modes and intermolecular interactions in crystal).
Acknowledgments This Letter was supported by the International Ph.D.-studies program at the Faculty of Chemistry Jagiellonian University within the Foundation for Polish Science MPD Program co-financed by the EU European Regional Development Fund. Results presented in this letter were partially obtained using PL-Grid Infrastructure and resources provided by ACC Cyfronet AGH. The authors thank Prof. Dušan Hadzˇi for valuable comments. M.J.W. also thank the University of Malaya and the UM.C/625/1/ HIR/MOHE/05 Grant, which supported his visit to Malaysia. References [1] D. Zanuy, C. Aleman, J. Phys. Chem. B 112 (2008) 3222.
[2] P.O. Löwdin, Ann. N. Y. Acad. Sci. 158 (1969) 86. [3] V.V. Prabhu, L. Young, K.W. Awati, W. Zhuang, E.W. Prohofsky, Phys. Rev. B 41 (1990) 7839. [4] M.M. Deshmukh, L.J. Bartolotti, S.R. Gadre, J. Comp. Chem. 32 (2011) 2996. [5] T. Kondo, A. Koschella, B. Heublein, D. Klemm, T. Heinze, Carbohydr. Res. 343 (2008) 2600. [6] V. Vicente, J. Martin, J. Jimenez-Barbero, J.L. Chiara, C. Vicent, Chemistry 10 (2004) 4240. [7] M. Hermansson, G. von Heijne, J. Mol. Biol. 334 (2003) 803. [8] L.S. Arneson, J.F. Katz, M. Liu, A.J. Sant, J. Immunol. 167 (2001) 6939. [9] G.C. Pimentel, A.L. McClellan, The Hydrogen Bond, Freeman Co, San Francisco, 1960. [10] P. Schuster, G. Zundel, C. Sandorfy (Eds.), The Hydrogen Bond. Recent Developments in Theory and Experiments, vols. I–III. North Holland Publishing Company, Amsterdam, 1976. [11] D. Hadzˇi (Ed.), Theoretical Treatments of Hydrogen Bonding, J. Wiley, Chichester, 1997. [12] S.J. Grabowski (Ed.), Hydrogen Bonding: New Insights, Springer, Dordrecht, 2006. [13] P. Gilli, G. Gilli, The Nature of the Hydrogen Bond, Oxford University Press, Oxford, 2009. [14] Y. Maréchal, A. Witkowski, J. Chem. Phys. 48 (1968) 3697. [15] A. Witkowski, M. Wójcik, Chem. Phys. 1 (1973) 9. [16] M.J. Wójcik, Int. J. Quant. Chem. 10 (1976) 747. [17] O. Henri-Rousseau, P. Blaise, Adv. Chem. Phys. 103 (1998) 1. [18] A.M. Yaremko, H. Ratajczak, J. Baran, A.J. Barnes, E.V. Mozdor, B. Silvi, Chem. Phys. 306 (2004) 57. [19] P. Blaise, M.J. Wójcik, O. Henri-Rousseau, J. Chem. Phys. 122 (2005) 064306. [20] C. Sandorfy, J. Mol. Struct. 790 (2006) 50. [21] M.Z. Brela, J. Stare, G. Pirc, M. Sollner Dolenc, M. Boczar, M.J. Wójcik, J. Mavri, J. Phys. Chem. B 116 (2012) 4510. [22] M. Boczar, J. Kwiendacz, M.J. Wójcik, J. Chem. Phys. 128 (2008) 164506. [23] M.J. Wójcik, J. Kwiendacz, M. Boczar, Ł. Boda, Y. Ozaki, Chem. Phys. 372 (2010) 72. [24] R. Errakhi, P. Meimoun, A. Lehner, G. Vidal, J. Briand, F. Corbineau, J.P. Rona, F. Bouteau, J. Exp. Bot. 59 (2008) 3121. [25] A. Witkowski, M. Wójcik, Chem. Phys. 21 (1977) 385. [26] J. de Villepin, A. Novak, F. Romain, Spectrochim. Acta 34A (1978) 1009. [27] J. de Villepin, A. Novak, F. Romain, Spectrochim. Acta 34A (1978) 1019. [28] M. Banno, K. Ohta, K. Tominaga, J. Phys. Chem. A 112 (2008) 4170. [29] V. Mohacˇek-Grošev, J. Grdadolnik, J. Stare, D. Hadzˇi, J. Raman Spectrosc. 40 (2009) 1605. [30] M. Boczar, R. Kurczab, M.J. Wójcik, Vib. Spectrosc. 52 (2010) 39. [31] T.M. Sabine, G.W. Cox, B.M. Craven, Acta Cryst. B25 (1969) 2437. [32] L. Leiserowitz, Acta Cryst. B32 (1976) 775. [33] P. Coppens, Acta Cryst. A37 (1981) C123. [34] R. Car, M. Parrinello, Phys. Rev. Lett. 55 (1985) 2471. [35] S.S. Iyengar, M.K. Petersen, T.J.F. Day, C.J. Burnham, V.E. Teige, G.A. Voth, J. Chem. Phys. 123 (2005) 084309. [36] J. Mavri, A. Jezierska, J.J. Panek, A. Koll, J. Chem. Phys. 126 (2007) 205101. [37] J. Stare, J.J. Panek, J. Eckert, J. Grdadolnik, J. Mavri, D. Hadzˇi, J. Phys. Chem. A 112 (2008) 1576. [38] P. Durlak, Z. Latajka, S. Berski, J. Chem. Phys. 131 (2009) 024308. [39] P. Durlak, Z. Latajka, Chem. Phys. Lett. 477 (2009) 249. [40] J. Mavri, G. Pirc, J. Stare, J. Chem. Phys. 132 (2010) 224506. [41] J. Kwiendacz, M. Boczar, M.J. Wójcik, Chem. Phys. Lett. 501 (2011) 623. [42] G.H. Wannier, J. Math. Phys. 19 (1978) 131. [43] CPMD, copyright IBM 1990–2008. Copyright MPI für Festkörperfoshung Stuttgart 1997–2001. [44] C.T. Lee, W.T. Yang, R.G. Parr, Phys. Rev. B 37 (1988) 785. [45] S. Goedecker, K. Maschke, Phys. Rev. A 45 (1992) 88. [46] H. Forbert, A. Kohlmeyer, Fourier, v 1.1, 2002–2008. [47] R. Vianello, J. Stare, J. Mavri, J. Grdadolnik, J. Zidar, Z.B. Maksic, J. Phys. Chem. B 115 (2011) 5999. [48] G. Pirc, J. Stare, J. Mavri, J. Chem. Phys. 132 (2010) 224506. [49] P. Durlak, Z. Latajka, J. Mol. Model. 17 (2011) 2159. [50] R. Vianello, M. Repicˇ, J. Mavri, Eur. J. Org. Chem. 2012 (2012) 7057. [51] N.V. Plotnikov, S.C.L. Kamerlin, A. Warshel, J. Phys. Chem. B 115 (2011) 7950.
Contents lists available at SciVerse ScienceDirect
Chemical Physics Letters journal homepage: www.elsevier.com/locate/cplett
Car–Parrinello simulation of the vibrational spectra of strong hydrogen bonds with isotopic substitution effects: Application to oxalic acid dihydrate Mateusz Z. Brela a,⇑, Marek J. Wójcik b,⇑, Marek Boczar b, Rauzah Hashim c a
Research Group of Molecular Modeling of Catalytic Processes, Faculty of Chemistry, Jagiellonian University, Ingardena 3, 30-060 Krakow, Poland Laboratory of Molecular Spectroscopy, Faculty of Chemistry, Jagiellonian University, Ingardena 3, 30-060 Krakow, Poland c Department of Chemistry, University of Malaya, 50603 Kuala Lumpur, Malaysia b
a r t i c l e
i n f o
Article history: Received 12 October 2012 In final form 18 December 2012 Available online 3 January 2013
a b s t r a c t The nature of strong intermolecular hydrogen bonding in oxalic acid dihydrate in the crystal phase was examined by infrared spectroscopy and Car–Parrinello molecular dynamics simulation. We studied region of infrared spectra associated with the O—H modes. The spectra were calculated using harmonic approximation with crystal field and time course of the dipole moment as obtained from Car–Parrinello simulation with quantization of the O—H motion, and isotopic substitution. We obtained good agreement of the molecular dynamic simulation with experiment. To our best knowledge, this is one of the first Car– Parrinello calculations of infrared spectra including anharmonicity effects and crystal field interactions. Ó 2012 Elsevier B.V. All rights reserved.
1. Introduction Hydrogen bond is one of the most important interactions in chemical and biochemical systems. For example in Nature, hydrogen bonds are responsible for the formations and functions of many organized complex structures such as DNA [1–3] biopolymers [4–6] and cell membranes [7,8]. The properties of hydrogen bonds have been widely investigated theoretically and experimentally and have been subject of several monographs [9–13]. Many theoretical models have been proposed to describe unusual features of hydrogen bond stretching bands in the IR spectra [14–23]. Oxalic acid is produced naturally from plants and has been shown to play an important role in the prevention of chronic diseases (osteoporosis, obesity), but its biological function in human body, almost unknown, is beginning to be investigated [24]. Infrared spectra of oxalic acid have been subject of several theoretical and experimental investigations, most of them focused on hydrogen bond stretching band [25–30]. Oxalic acid dihydrate forms at room temperature white crystals with melting point at 394 K. The crystal structure of oxalic acid dihydrate was determined from neutron diffraction by Sabine et al. [31] and from X-ray by Leiserowitz [32] and Coppens [33]. The crystal of oxalic acid dihydrate belongs to the P21/n space group of symmetry. The unit cell of oxalic acid dihydrate contains two molecules of oxalic acid and four molecules of water. The crystal structure is supported by two types of hydrogen bonds (Fig⇑ Corresponding authors. E-mail address: [email protected] (M.J. Wójcik). 0009-2614/$ - see front matter Ó 2012 Elsevier B.V. All rights reserved. http://dx.doi.org/10.1016/j.cplett.2012.12.035
ure 1A). One type connects the water molecules as donors and the acid carbonyl groups as acceptors; it may be classified as normal for its length in the range of 2.869 Å. The other type connects the acid’s hydroxyl group as the donor to water oxygen as acceptor. This bond is very short; at room temperature the O O distance is 2.499 Å. The OH bond length is 1.041 Å at 298 K and it expands to 1.071 Å at 80 K. The recently investigated single crystal Raman spectra of oxalic acid dihydrate and its analog [29] have shown interesting features thus rising interest in a more detailed study of IR spectra and their theoretical modelling. The latter is of even broader interest because of the rarity of theoretical modelling of the spectra of examples of hydrogen bonds of O O distances about 2.50 Å. In this Letter we present results of ab initio Car–Parrinello molecular dynamics (CPMD) calculations [34] performed for the unit cell of the crystal of oxalic acid dihydrate. The vibrational spectrum of the oxalic acid dihydrate has been calculated, analyzed, and compared with experimental spectra. The greatest emphasis has been put on the reconstruction of the hydrogen bond vibrational bands. The mO–H band has been thoroughly examined, and the simulation of the mO–H bandshape, based on CPMD calculations, has been performed. In the first step, we optimized the geometry of the unit cell of oxalic acid dihydrate, then we calculated trajectories using the Car–Parrinello molecular dynamics (CPMD) scheme. CPMD calculations have recently become a popular tool for the interpretation of infrared spectra of hydrogen-bonded systems [35–41]. In the present Letter we took into account the periodicity of the crystal and thermal fluctuations. The second step involved quantization
89
M.Z. Brela et al. / Chemical Physics Letters 558 (2013) 88–92
Figure 1. Structure and atom labeling in crystal unit cell of oxalic acid dihydrate. The atoms are color coded: carbon = gray, oxygen = red, hydrogen = white. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)
of the proton motion, for which we used an alternative method based on an a posteriori quantization of the nuclear degrees of freedom associated with proton dynamics. The vibrational spectrum of the oxalic acid dihydrate crystal has been calculated, analyzed, and compared with experimental data and also with the results of DFT calculations performed for a crystal cell of the oxalic acid dihydrate. The largest emphasis has been placed on the reconstruction of the stretching bands of hydrogen bonded O—H groups. The tO–H bands, obtained from CPMD calculations and instantaneous one-dimensional snapshot potentials, extracted from the trajectory calculated using the Wannier functions [42] were compared with the experimental bands of the oxalic acid dihydrate crystal. We discuss also the isotopic effects on the spectra.
then calculated using the proton potentials. For each of the extracted snapshot structures, a one-dimensional (1D) proton potential function was obtained by stepwise displacement of one hydrogen-bonded proton along the hydrogen vector parallel to the pertinent O—H line. The corresponding internal coordinate (x) was defined as the distance between hydrogen and oxygen atoms. The scan typically covered the x range from 0.5 Å to 0.5 Å in relation to the equilibrium O—H distance. The hydrogen was displaced along the distance steps of 3.33% of the range. It is worth to point out that the proton potentials are ensemble averages and are hence mass independent in the limit of classical statistical mechanics. Therefore the snapshot methodology for deuterated species is applicable on the trajectory obtained by integrating equations of motion for hydrogen form.
2. Computational methods
3. Results and discussion
We considered several isotopically substituted molecules. Geometry optimizations and vibrational frequency calculations were performed to characterize the key stationary points on the potential energy surface (PES) for unit cell of oxalic acid dihydrate (Figure 1B). We have also performed a vibrational analysis in the optimized crystal structure using the (CPMD) program package version 3.13.2 [43]. The BLYP [44] density functional together with a plane-wave basis set with a kinetic energy cutoff of 300 Ry were used in our study. The Goedecker atomic pseudopotentials [45] were used for the treatment of the core electrons. The starting point for our molecular dynamics (MD) simulation was the neutron structure data obtained by Sabine et al. [31]. Three-dimensional periodicity was fully implemented in our calculation. The MD simulation was carried out in the unit cell fixed to the experimentally determined size and shape and using the Car– Parrinello scheme [34]. The fictitious orbital mass was set to 300 a.u, the temperature was 300 K and the molecular dynamics time step was set to 1 a.u. (about 0.024188 fs), while other properties were set the same as the static calculation. The total simulation time was about 10 ps. The simulation was performed on an HP Cluster Platform 3000 BL 2x220 PC/LINUX and it lasted for about 720 days using only 12 cores. The anharmonic vibrational spectrum was calculated using the Fourier transformation of the dipole autocorrelation function obtained from dipole trajectories generated by the CPMD simulation facilitated by the scripts of Kohlmeyer [46]. The dipole dynamics were calculated using the optimally localized Wannier functions [42]. Snapshot structures of the trajectory were extracted every 4000 steps (about 97 fs), yielding 125 distinct structures, representative of the fluctuating crystalline environment. Note that due to the fact that two molecules which interacted in four hydrogen bonds in the unit cell, the total number of distinct hydrogen bonded moieties extracted was 500. The contour of the O—H stretching band was
3.1. Geometry analysis The results of our study are presented in the following order: analysis of the structure, and vibrational analysis. In Figure 1A we have shown the structure and atom labeling of oxalic acid dihydrate. The geometry results of static calculations are shown in Table 1. Analysis of the Car–Parrinello trajectory for hydrogen bond formed by O–H group of oxalic acid and water shows a time-course of the OO and O—H distances in the O—H O hydrogen bond and shows that the proton is asymmetrically located nearer to the donor oxygen. The average O O distance is 2.51 Å. This suggests that the hydrogen bonds are strong. Over the entire CPMD run the O O distance spans the range between 2.44 and 2.58 Å. It is slightly longer than the experimental crystallographic value of 2.52 Å. The diffraction experiment [33] gave an O—H distance of 1.03 Å. Our DFTs calculations predicted that this bond is longer, with an average bond length 1.05 Å. Generally our CPMD static calculations reproduce the geometric parameters for hydrogen bonds in the crystal reasonably well. 3.2. Infrared spectra of crystal structure calculated in harmonic approximation and autocorrelation function obtained after CPMD Figure 2 presents crystal spectrum calculated in harmonic approximation. The experimental infrared spectrum of oxalic acid
Table 1 The calculated hydrogen bond geometry for four hydrogen bonds in the crystal unit cell. The DFT calculations were performed in the crystal field. See Figure 1A for symbol description.
O—H O O
I-HB
II-HB
III-HB
IV-HB
Average
Exp. [33]
1.01 2.69
1.09 2.50
1.09 2.45
1.01 2.56
1.05 2.55
1.05 2.52
90
M.Z. Brela et al. / Chemical Physics Letters 558 (2013) 88–92
Figure 2. Infrared spectrum of oxalic acid dihydrate crystal calculated in harmonic approximation. Insert presents fragment of crystal structure of the unit cell. The atoms are color coded: carbon = gray, oxygen = red, hydrogen = white. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)
dihydrate was presented in the paper Mohacˇek-Grošev et al. [29]. In the static calculation we obtained peaks, which represent stretching modes of OH groups. Our calculations shows that they represent hydrogen-bonded O–H stretching in oxalic acid (2300 cm 1), weakly hydrogen-bonded O–H stretching in water (2900 and 3100 cm 1) and O–H stretching of free groups in water (3300 and 3450 cm 1). The band at 2300 cm 1 represents the mixing of two kinds of modes – stretching and twisting carboxyl
Figure 3. Infrared spectra of hydrogen form of oxalic acid dihydrate. Inset presents fragment of crystal structure of the unit cell. The atoms are color coded: carbon = gray, oxygen = red, hydrogen = white. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)
group. Harmonic calculations cannot reproduce the broadening and intensity increase caused by the anharmonic intermolecular couplings present in hydrogen bonds [14–16]. In harmonic quantum chemical calculations one does not obtain bandwidths and vibrational transitions are represented by d functions. The bandshapes presented in Figure 2 were obtained by superimposing a Gaussian function with a half-width 50 cm 1 on each calculated transition; such representations have often been used, e.g. in articles [21–23].
Figure 4. Infrared spectra of the four kind isotopic substitution of oxalic acid dihydrate with crystal structure of the unit cell. The atoms are color coded: carbon = gray, oxygen = red, hydrogen = white, deuterium = green. For details see text in Section 3.4 ‘‘Isotope effects’’. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)
91
M.Z. Brela et al. / Chemical Physics Letters 558 (2013) 88–92
Table 2 Comparison between frequencies of selected vibrations in oxalic acid dihydrate crystal obtained from CPMD static calculations (harmonic frequencies obtained using BLYP in crystal field) and experimental data from paper Mohacˇek-Grošev et al. [29].
3.3. Snapshot potential energy functions We calculated one-dimensional (1D) proton potentials corresponding to the instantaneous snapshot structures extracted from the CPMD trajectory. It should be noted that typical 1D scan requires 30 points at which the energy is evaluated. We compared the bands calculated with the snapshot methodology [47,48] and band extracted from the powder spectrum obtained using Fourier transformation. We generated one-dimensional proton potentials using one linear pathway. Figure 5 shows schematically the proton potential construction and obtained potentials together with the solution of the vibrational Schrödinger equation. More than 90% of frequencies, which were calculated by using the three types of 1D potentials, span the range from 1750 to 2250 cm 1. The contour of the band is similar to experimental spectrum. It shows the small anharmonicity and structure of modes. In general, snapshot proton potentials constructed following the O—H direction are practically not shifted. This strategy also yields a narrower band contour with a half-width of about 230 cm 1. 3.4. Isotope effects We have also performed the analysis of spectra for oxalic acid dihydrate and its selected deuterated isotopomers. Figure 4 shows
Experimental frequencies (cm 1) [29]
Description
115 1175 1177 1193 1217
154R 1133
m (O O) c (O–H O)
1289 1318 1338 1345
1253
d (C–O–H)
m (C@O)
1430 1440 1504 1513
Figure 5. O—H stretching band contours of oxalic acid dihydrate: Black line contour was calculated from individual fundamental vibrational transitions as superposition of Gaussian functions with a half width of 10 cm 1. Gray d functions represent fundamental vibrational transitions. The red contour represents the spectrum calculated by Fourier transform of the autocorrelation function of the atoms position obtained from Car–Parrinello trajectory. The bands were normalized. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)
The second methodology, which we used to obtain these spectra, was the Fourier transform of the autocorrelation function. We performed calculations with two types of functions: the total dipole moment (Figures 3 and 4) and function describing atom positions to reproduce data for the analysis of a single band O–H (Figure 5). This approximation gives qualitative reproduction of the experimental spectra, though their positions are shifted to higher frequencies for all O–H bands. The region of high frequency is composed of three peaks, centered at 1750–2100 cm 1, 3100– 3350 cm 1 and small peak near 3600–3700 cm 1. Our calculations show that, they represent hydrogen bonded O–H groups in oxalic acid, weakly hydrogen-bonded O–H groups in water, and free O– H groups in water. The agreement with the experimental broad band with a maximum at 1900 cm 1 [29] which represents the O–H stretching band (strongly hydrogen-bonded), is reasonably good. The frequencies of selected vibration are shown in Table 2.
Calculated frequencies (cm 1)
m
1570 1588 1600 1616
1620
d (O H)water
1716 1722 1734 1758
1690
m (C@O)
2358 2454 2496 2676
1870 2009
m (O H)
3343 3346 3388 3390 3639 3688 3692 3721
3431 3489
m (O H)water
Stretching, d
in-plane bending, c
out-of plane bending, R – Raman.
the infrared spectra obtained using the autocorrelation functions of dipole moment (CPMD simulation) for selected deuterated isotopomers of oxalic acid dihydrate. For all these structures we performed independent CPMD simulations with the same parameters (see the computational part of this article). Analysis of calculated spectra allows us to study the interaction involving O—H stretching vibrations in the crystal. Figure 4A shows the infrared spectrum of the crystal with only one deuterium replacing hydrogen. We observed the stretching mode of O–D near 1700– 1850 cm 1. Additional effect of this substitution is increase in intensity of the band near 3600 cm 1 relative to the band near 3300 cm 1. The second structure was obtained by replacing all hydrogens except one in oxalic acid by deuterium atoms. Figure 4B. shows the result of this substitution. We obtained broad bands near 1600–1900 and 2200–2400 cm 1. The first one represents O–D vibrations of oxalic acid and the second one is composed of O—H vibrations of oxalic acid and O–D vibrations of water molecules. Figure 4C shows the spectrum of the crystal with deuterated oxalic acid and non-deuterated water molecules. In this spectrum the band 2200–2400 cm 1 of non-deuterated oxalic acid dihydrate (Figure 3) is shifted to 1600–1800 cm 1. The O—H stretching vibrations of water in this form have frequencies close to O—H frequencies of non-deuterated oxalic acid dihydrate. Figure 4D shows the spectra calculated for per-deuterated oxalic acid dihydrate. The O—D stretching bands for oxalic acid are similar as in Figure 4C but their center is shifted to slightly higher frequency. The O—D stretching modes of water molecule are observed near 2400 and 2700 cm 1. This result shows that the coupling of O—H vibrations
92
M.Z. Brela et al. / Chemical Physics Letters 558 (2013) 88–92
is observed only in non-deuterated oxalic acid dihydrate. The coupling occurs between O–H groups present in oxalic acid or water molecules and there is no significant coupling between OH groups of oxalic acid and water. This can be explained by difference of OH vibrational frequencies in water and oxalic acid. Similar effect was obtained by Latajka and Durlak in the gas phase [39,49]. This kind of phenomena gives perspectives for QM/MM studies of nuclear quantum effects for example in enzyme centers [50,51]. 4. Conclusions In our Letter we present results of the simulation of spectra of the crystal of oxalic acid dihydrate. The results can be summarized as follows: 1. Infrared spectrum calculated by CPMD method reproduces reasonably well the high-frequency bands in the experimental spectrum. Also the deuterium shifts of vibrational bands were calculated in good accordance with experimental data due to the application of dipole momentum dynamics. 2. The structure of the O–H stretching band is reproduced fairly well. In this part of our Letter, we investigated the snapshot methodology. Calculated bands describing hydrogen bonds vibrations better reproduce experimental values than results of static ab initio calculations at BLYP in the crystal field level performed for the crystal cell. 3. The identification of bands in calculated infrared spectrum is difficult in FT calculations but snapshots methodology can be useful in description of single bands. 4. Calculations for isotopically substituted molecules allow to study couplings effects and intermolecular interactions in crystal of oxalic acid dihydrate. We focused our research only on stretching bands. We observe couplings of similar modes. 5. Improvement in calculation results proves that CPMD method is adequate for investigations of systems with hydrogen bonds as it takes into account most of mechanisms determining the hydrogen bonds dynamics (anharmonicity, couplings between vibrational modes and intermolecular interactions in crystal).
Acknowledgments This Letter was supported by the International Ph.D.-studies program at the Faculty of Chemistry Jagiellonian University within the Foundation for Polish Science MPD Program co-financed by the EU European Regional Development Fund. Results presented in this letter were partially obtained using PL-Grid Infrastructure and resources provided by ACC Cyfronet AGH. The authors thank Prof. Dušan Hadzˇi for valuable comments. M.J.W. also thank the University of Malaya and the UM.C/625/1/ HIR/MOHE/05 Grant, which supported his visit to Malaysia. References [1] D. Zanuy, C. Aleman, J. Phys. Chem. B 112 (2008) 3222.
[2] P.O. Löwdin, Ann. N. Y. Acad. Sci. 158 (1969) 86. [3] V.V. Prabhu, L. Young, K.W. Awati, W. Zhuang, E.W. Prohofsky, Phys. Rev. B 41 (1990) 7839. [4] M.M. Deshmukh, L.J. Bartolotti, S.R. Gadre, J. Comp. Chem. 32 (2011) 2996. [5] T. Kondo, A. Koschella, B. Heublein, D. Klemm, T. Heinze, Carbohydr. Res. 343 (2008) 2600. [6] V. Vicente, J. Martin, J. Jimenez-Barbero, J.L. Chiara, C. Vicent, Chemistry 10 (2004) 4240. [7] M. Hermansson, G. von Heijne, J. Mol. Biol. 334 (2003) 803. [8] L.S. Arneson, J.F. Katz, M. Liu, A.J. Sant, J. Immunol. 167 (2001) 6939. [9] G.C. Pimentel, A.L. McClellan, The Hydrogen Bond, Freeman Co, San Francisco, 1960. [10] P. Schuster, G. Zundel, C. Sandorfy (Eds.), The Hydrogen Bond. Recent Developments in Theory and Experiments, vols. I–III. North Holland Publishing Company, Amsterdam, 1976. [11] D. Hadzˇi (Ed.), Theoretical Treatments of Hydrogen Bonding, J. Wiley, Chichester, 1997. [12] S.J. Grabowski (Ed.), Hydrogen Bonding: New Insights, Springer, Dordrecht, 2006. [13] P. Gilli, G. Gilli, The Nature of the Hydrogen Bond, Oxford University Press, Oxford, 2009. [14] Y. Maréchal, A. Witkowski, J. Chem. Phys. 48 (1968) 3697. [15] A. Witkowski, M. Wójcik, Chem. Phys. 1 (1973) 9. [16] M.J. Wójcik, Int. J. Quant. Chem. 10 (1976) 747. [17] O. Henri-Rousseau, P. Blaise, Adv. Chem. Phys. 103 (1998) 1. [18] A.M. Yaremko, H. Ratajczak, J. Baran, A.J. Barnes, E.V. Mozdor, B. Silvi, Chem. Phys. 306 (2004) 57. [19] P. Blaise, M.J. Wójcik, O. Henri-Rousseau, J. Chem. Phys. 122 (2005) 064306. [20] C. Sandorfy, J. Mol. Struct. 790 (2006) 50. [21] M.Z. Brela, J. Stare, G. Pirc, M. Sollner Dolenc, M. Boczar, M.J. Wójcik, J. Mavri, J. Phys. Chem. B 116 (2012) 4510. [22] M. Boczar, J. Kwiendacz, M.J. Wójcik, J. Chem. Phys. 128 (2008) 164506. [23] M.J. Wójcik, J. Kwiendacz, M. Boczar, Ł. Boda, Y. Ozaki, Chem. Phys. 372 (2010) 72. [24] R. Errakhi, P. Meimoun, A. Lehner, G. Vidal, J. Briand, F. Corbineau, J.P. Rona, F. Bouteau, J. Exp. Bot. 59 (2008) 3121. [25] A. Witkowski, M. Wójcik, Chem. Phys. 21 (1977) 385. [26] J. de Villepin, A. Novak, F. Romain, Spectrochim. Acta 34A (1978) 1009. [27] J. de Villepin, A. Novak, F. Romain, Spectrochim. Acta 34A (1978) 1019. [28] M. Banno, K. Ohta, K. Tominaga, J. Phys. Chem. A 112 (2008) 4170. [29] V. Mohacˇek-Grošev, J. Grdadolnik, J. Stare, D. Hadzˇi, J. Raman Spectrosc. 40 (2009) 1605. [30] M. Boczar, R. Kurczab, M.J. Wójcik, Vib. Spectrosc. 52 (2010) 39. [31] T.M. Sabine, G.W. Cox, B.M. Craven, Acta Cryst. B25 (1969) 2437. [32] L. Leiserowitz, Acta Cryst. B32 (1976) 775. [33] P. Coppens, Acta Cryst. A37 (1981) C123. [34] R. Car, M. Parrinello, Phys. Rev. Lett. 55 (1985) 2471. [35] S.S. Iyengar, M.K. Petersen, T.J.F. Day, C.J. Burnham, V.E. Teige, G.A. Voth, J. Chem. Phys. 123 (2005) 084309. [36] J. Mavri, A. Jezierska, J.J. Panek, A. Koll, J. Chem. Phys. 126 (2007) 205101. [37] J. Stare, J.J. Panek, J. Eckert, J. Grdadolnik, J. Mavri, D. Hadzˇi, J. Phys. Chem. A 112 (2008) 1576. [38] P. Durlak, Z. Latajka, S. Berski, J. Chem. Phys. 131 (2009) 024308. [39] P. Durlak, Z. Latajka, Chem. Phys. Lett. 477 (2009) 249. [40] J. Mavri, G. Pirc, J. Stare, J. Chem. Phys. 132 (2010) 224506. [41] J. Kwiendacz, M. Boczar, M.J. Wójcik, Chem. Phys. Lett. 501 (2011) 623. [42] G.H. Wannier, J. Math. Phys. 19 (1978) 131. [43] CPMD, copyright IBM 1990–2008. Copyright MPI für Festkörperfoshung Stuttgart 1997–2001. [44] C.T. Lee, W.T. Yang, R.G. Parr, Phys. Rev. B 37 (1988) 785. [45] S. Goedecker, K. Maschke, Phys. Rev. A 45 (1992) 88. [46] H. Forbert, A. Kohlmeyer, Fourier, v 1.1, 2002–2008. [47] R. Vianello, J. Stare, J. Mavri, J. Grdadolnik, J. Zidar, Z.B. Maksic, J. Phys. Chem. B 115 (2011) 5999. [48] G. Pirc, J. Stare, J. Mavri, J. Chem. Phys. 132 (2010) 224506. [49] P. Durlak, Z. Latajka, J. Mol. Model. 17 (2011) 2159. [50] R. Vianello, M. Repicˇ, J. Mavri, Eur. J. Org. Chem. 2012 (2012) 7057. [51] N.V. Plotnikov, S.C.L. Kamerlin, A. Warshel, J. Phys. Chem. B 115 (2011) 7950.