Structure and vibrational dynamics of hydrogen bond in hydrogenbissulfate anion in the gas phase and in the solid state: A computational study

Available online at www.sciencedirect.com

Journal of Molecular Structure 844–845 (2007) 215–221 www.elsevier.com/locate/molstruc

Structure and vibrational dynamics of hydrogen bond in hydrogenbissulfate anion in the gas phase and in the solid state: A computational study Katarina Demsˇar a

a,b

, Jernej Stare

b,c

, Janez Mavri

b,*

Department of Inorganic Chemistry, Faculty of Chemistry and Chemical Technology, University of Ljubljana, Slovenia b National Institute of Chemistry, Ljubljana, Slovenia c T-CNLS/LANSCE-LC, Los Alamos National Laboratory, NM, USA Received 1 March 2007; accepted 13 March 2007 Available online 24 March 2007 Dedicated to Professor Lucjan Sobczyk on the occasion of his 80th birthday.

Abstract Nature of a short hydrogen bond in hydrogenbissulfate ion was studied by quantum-chemical methods. Isolated triple charged anion is predicted to be unstable on the Hartree–Fock and Density Functional Theory levels using Pople basis sets with the exception of some small-size basis sets. For stable structures the predicted metric parameters are in poor agreement with the experiment and vary significantly with the basis set, indicating insufficient level of theory. The same is true for the calculated frequencies, in particular for the OAH stretching fundamental that seems to be too high. Inclusion of the crystal field effects by imposing the periodicity in conjunction with the plane wave basis sets improves agreement with the neutron diffraction data considerably. Exception is the OAH distance that was predicted to be too short. By quantization of the OAH coordinate an almost perfect agreement with the experimental OAH was found. Besides, quantization of the OAH coordinate shifts the calculated vibrational frequency toward the experimental center of the broad OAH stretching infrared absorption. Perspectives for application and extensions of presented methodology for vibrational spectroscopy and (bio)catalysis are given.  2007 Elsevier B.V. All rights reserved. Keywords: Hydrogen bond; Solid state; Vibrational spectroscopy; Anharmonicity; Quantum chemistry; Plane-wave basis set; Periodicity

1. Introduction Hydrogen bonding is an important type of molecular interaction that plays a vital role for structure and functionality of many relevant materials and biological systems. There is much interest to understand short, strong hydrogen bonds. The possibility of proton transfer in such hydrogen bonds gives rise to applications for as diverse systems as proton conductors [1], ferroelectrics [2], electrode processes, crystal engineering [3], and enzymology [4–7]. It is for the sake of high applicability in a variety of modern *

Corresponding author. Tel.: +386 1 476 0309; fax: +386 1 476 0300. E-mail address: [email protected] (J. Mavri).

0022-2860/$ - see front matter  2007 Elsevier B.V. All rights reserved. doi:10.1016/j.molstruc.2007.03.030

technologies that hydrogen bonding has been subjected to extensive experimental and computational studies. Despite enormous efforts, however, many aspects of strong and short hydrogen bond, particularly in the solid state, are not fully understood and represent a great challenge for fundamental research. Among features that are in the focus of the hydrogen bond research are the effects of polar environment or crystal field that often play an important role for the structure and energetics of hydrogen bond moieties [8–10]. Probably the most important effect that accompanies systems containing short and strong hydrogen bonds is the anharmonicity of the potential energy surface, in particular the part that corresponds to the proton motion [11]. Additionally, in such systems, the proton motion is

216

K. Demsˇar et al. / Journal of Molecular Structure 844–845 (2007) 215–221

essentially quantum in its nature, making the computational treatment a very challenging problem that is based on the solving of the vibrational Schro¨dinger equation. A number of possible computational strategies are available for this purpose and all share the problem of high complexity and high computational resources required, most often resulting in the necessity to simplify the model to only a few vibrational degrees of freedom. For an insight into these methods the reader is referred to refs. [12–21]. It is thus evident that the prediction of vibrational spectra of strongly anharmonic systems in the condensed phase is a challenge for computational chemistry and physics. It is essential to perform calculations on sufficiently high level of theory and to properly include the effects of the environment. The shape of the proton potential is very sensitive to the applied model. Typically one has to proceed with large, flexible basis sets and inclusion of the electron correlation. Crystal environment can be considered by methods that include periodicity. Inclusion of fluctuations is possible by considering thermal averaging but it is CPU demanding [22–24]. Among crystalline systems that contain short hydrogen bonds the family of alkali hydrogenbissulfates, M3H(SO4)2, where M is an alkali metal (from this point on the hydrogenbissulfate anion will be referred to as HBS3), appears to be highly relevant. Deuterated rubidium HBS has a strong intermolecular hydrogen bond. In its deuterated form it is antiferroelectric below 25 K while the undeuterated form remains paraelectric to 4 K [25]. There were suggestions that this isotope effect is due to Ubbelohde effect [26], i.e., due to change of the O. . .O distance upon deuteration. A particularly interesting aspect of alkali HBS’s are their vibrational spectra, which feature extremely broad absorptions attributed to the OAH stretching mode. For instance, an infrared study of rubidium hydrogenbissulfate [27] places the OAH stretching absorption between 2000 cm1 and the low frequency limit of record. The series of spectra reveals no significant effect of temperature. There are however some changes in transmissions indicating the phase transitions. Ratajczak and Yaremko explained extreme broadening and structuring of asymmetric stretching of short hydrogen bonds in solids with strong interaction of the high frequency mode mOH with lattice modes [28]. In addition, Fermi resonances are often present. In other words the proton potential fluctuates with time and with space, giving rise to broadening of the mOH band. About the same conclusions are valid for the sodium HBS, Na3HBS for which an extremely short ˚ has hydrogen bond with the O. . .O distance of 2.40–2.43 A been determined by crystallographic methods [29,30]. In this article we examined the structure and vibrational properties of the HBS3 anion (Fig. 1) with several electron structure methods. Starting from basic description, we used an isolated model in conjunction with the most common quantum-chemical methods together with various basis sets. The issue of stability of the free HBS3 has been considered in terms of basis set incompleteness. As we found

Fig. 1. A model of the free hydrogenbissulfate anion with selected interatomic distances.

that the environmental effects imposed by the crystal field are crucial for the proper modeling of the HBS3 anion, we modeled the solid structure of Na3HBS imposing full periodicity and studied the influence of the crystal point group on the resulting minimized geometry and vibrational frequencies. For both isolated and periodical model we monitored geometric parameters and vibrational frequencies and we compared them with the experimental values. Moreover, for the periodical model, we constructed the one-dimensional proton potential for which we solved the vibrational Schro¨dinger equation for the asymmetric OAH motion. Perspectives for computational support of the structural parameters and vibrational spectroscopy of the hydrogen bonding in the crystalline solid state. 2. Methods and models 2.1. Gas-phase calculations All calculations on the free HBS3 anion (Fig. 1) were carried out by the Gaussian 03 program package [31] using standard ab initio and Density Functional Theory (DFT) electron structure methods (HF and B3LYP), various Pople basis sets (from small 3-21G to large 6311++G(2d,2p)) and default settings for geometry optimization. Consistency of each geometry optimization was verified by harmonic frequency analysis. Stability issue of a free HBS3 anion as a function of basis set size and level of electron correlation was considered accordingly (see below). 2.2. Crystal field calculations Crystalline Na3HBS was modeled according to the Xray and neutron crystallographic structure data of Joswig ˚, and Fuess [30] with monoclinic unit cell (a = 8.648 A ˚ ˚ b = 9.648 A, c = 9.143 A, b = 108.77, space group P21/c, Z = 4; see Fig. 2). Periodicity and symmetry operators within the unit cell have been fully considered, resulting in four equivalent HBS3 units in the cell. Optionally, lower symmetry P1 has been employed, rendering all the atoms non-equivalent, in order to improve geometry optimization and vibrational analysis in the harmonic approximation (see below). The BLYP density functional (frozen

K. Demsˇar et al. / Journal of Molecular Structure 844–845 (2007) 215–221

217

calculation was performed for each of the models considered. 2.3. Proton potentials and anharmonic OAH stretching frequencies

Fig. 2. ORTEP visualization of the unit cell of crystalline sodium hydrogenbissulfate (ref. [30]). Color legend: blue – sodium; red: oxygen; brown – sulfur; small white circles – hydrogen. (For interpretation of the references in color in this figure legend, the reader is referred to the web version of this article.)

Fig. 3. Anharmonic proton potential for the OAH stretching motion in the solid sodium hydrogenbissulfate together with the ground state and the first excited state vibrational levels and wavefunctions and the fundamental excitation frequency.

core approximation and explicit valence electrons) was employed for the electron structure calculations together with plane wave basis set of a kinetic energy cutoff of 90 or 120 Ry and Goedecker relativistic atomic pseudopotentials [32]. All the solid-state calculations were performed in the C-point approximation, which is a reasonable choice due to the fact that the unit cell is relatively large; besides, the usage of a supercell or k-point sampling strategies are prevented by the size of the unit cell. Note that the unit cell parameters were fixed to the experimentally determined values and were not changed during the optimization. The Car-Parrinello program package CPMD v. 3.9 was used [33]. Geometry optimization and harmonic frequency

One-dimensional single point proton potentials were calculated for periodical models in the following way. The proton was pointwise displaced along a straight line from one side of the hydrogen bond to another, keeping the OAH. . .O moiety linear all the time. Since the deviation from linearity within the OAH. . .O moiety is of the order of a few tenths of a degree in all cases, the impact of the assumed linearity to the potential function can be considered to be negligible. For each point the energy was calculated using the corresponding model chemistry. During the scan all other nuclei were fixed to their optimized positions. This corresponds to the physical fact that the protonic motion is fast relative to other degrees of freedom (or at least to a considerable part of them) and feels the environment effectively to be frozen. Such a choice of the scan strategy is justified by the previously published experience with protonic vibrational modes of hydrogen-bonded systems [18,34]. The proton was dis˚ steps with the restriction that its distance placed in 0.05 A ˚. to the nearer oxygen atom was not shorter than 0.80 A The periodical model includes four Na3HBS formula units in the simulation cell (Fig. 2), hence it is possible to move up to four hydrogen atoms at a time in the corresponding hydrogen bonds. We found that the motions of the four protons in the unit cell practically did not interfere, i.e., the energy profile was effectively the same regardless of the number of displaced protons and their relative phases. This is supported by the fact that the experimental distance between the nearest protons is ˚ , effectively reducing the many-quantum-body 4.588 A problem to be only one-dimensional in its nature. Having acquired the proton potential function, we solved the corresponding vibrational Schro¨dinger equation by the Fourier Grid Hamiltonian (FGH) method [15]. The FGH approach is based on discretization of the coordinate space into uniform intervals that comprise the orthogonal vibrational basis set. The Hamiltonian matrix elements are calculated by double Fourier transform which allows for a simple formulation for both the kinetic and potential parts of the Hamiltonian. We utilized a special version of the FGH method that has been optimized for application in the internal coordinate formalism [35]. The gridded vibrational basis set spanned the OAH dis˚ and was uniformly discrettance range from 0.6 to 1.8 A ized to 201 points. The reduced mass of the proton motion was set to 0.995 amu. The solution of the onedimensional vibrational Schro¨dinger equation yields discrete vibrational levels (together with their corresponding wavefunctions and energies), giving rise to the OAH stretching frequency beyond the harmonic approximation (Fig. 3).

218

K. Demsˇar et al. / Journal of Molecular Structure 844–845 (2007) 215–221

3. Results and discussion

3.2. Crystal field calculations

3.1. Gas-phase calculations

Calculated geometric parameters listed in Table 2 feature only minor differences between the applied basis set size, convergence criterion, or crystal symmetry. The calculated O. . .O distance is in the range between 2.439 and ˚ , which is in very good agreement with the experi2.443 A ˚ . Even when the P21/c symmetry mental value of 2.432 A constraints are turned off and P1 is used instead, the difference between the four OAH. . .O moieties is negligible (note that for the P21/c symmetry only one value is given for each geometric parameter, since all four hydrogen bonds in the unit cell are equivalent; in contrast, they become formally inequivalent with the P1 symmetry). The calculated ˚ appears to be too short OAH distance of about 1.094 A as compared to the neutron-diffraction-determined value ˚ . This disaccord can be explained through the of 1.156 A nuclear quantum effect of the proton (see below). As far as geometric parameters are concerned, the choice of the applied symmetry, basis set size and convergence criterion does not seem to be crucial – in any case the experimental structure can be reasonably reproduced no matter what combination of the calculation settings is used. This is not true however for the calculated harmonic frequencies, particularly in the sense that harmonic frequencies reflect the consistency of a geometry optimization. Application of vibrational analysis to the structure obtained by the CPMD default convergence criterion and a plane-wave cutoff of 90 Ry yields as many as 13 modes with imaginary frequencies, 9 of which are located at around 1000i cm1, and 4 at around 200i cm1 (see Table 2). Similar results are obtained with the larger plane-wave cutoff of 120 Ry, except that the imaginary frequencies are smaller in magnitude (9 at around 400i cm1, and 4 at around 100i cm1). Notable is the fact that the former 9 modes include only contributions of the four sodium atoms that are resident on the inversion centers; in contrast, the latter 4 modes predominantly include motion of sodium ions located at general positions and also small contributions of oxygen, sulfur and hydrogen nuclei. This suggests that the so obtained optimized structure does not correspond to a true minimum on the potential energy hypersurface; in particular the positions of sodium atoms

In Table 1 are shown the calculated geometric parameters and characteristic harmonic OAH stretching frequencies for the isolated HBS3 ion (Fig. 1) calculated at various levels of theory. It is interesting that the HBS3 ion is predicted to be stable only for relatively small basis sets. Namely, when using the 6-31+G(d,p) basis set and larger (together with the B3LYP method) the ion is not stable and dissociates to sulfate and hydrogensulfate anions during the geometry optimization. The fact that the HBS3 ion is stable when using sufficiently small basis sets can be explained by the basis set incompleteness which results in the basis set superposition error (BSSE). Namely, because of the basis set incompleteness, electrons from either entity of the molecular complex (sulfate and hydrogensulfate ion in the present case) tend to be more delocalized into orbitals of the other part. This results in an unphysical attractive force between the entities, called the Pulay force [36]; the smaller the basis set, the larger the BSSE and Pulay force. As a result, the attraction due to BSSE may compensate the Coulomb repulsion  between the charged components ðSO2 4 and HSO4 Þ, making the complex virtually stable. Further support for the assumed role of BSSE in the sought stability of the HBS3 ion can be deduced from the fact that when using the noncorrelated Hartree–Fock method instead of the correlated B3LYP the critical basis set size, beyond which the ion dissociates, is decreased from 6-31+G(d,p) to 6-31G(d), which clearly reflects the well known feature that BSSE is smaller with non-correlated methods [37]. We can conclude that at sufficiently high level of theory the free HBS3 is not stable, suggesting the necessity to avoid gas phase models. In addition, as seen from Table 1, the calculated geometric parameters and harmonic OAH stretching frequencies vary significantly with the applied level of theory even when the calculation predicts the HBS3 ion to be stable, indicating that the description in terms of the basis set size and flexibility is poor. For that sake we did not proceed with calculations of the OAH stretching frequency beyond the harmonic approximation.

Table 1 Optimized geometry parameters of hydrogen bond and harmonic OAH stretching frequencies in hydrogenbissulfate anion calculated for gas phase model with the Hartree–Fock and the Density Functional Theory B3LYP method together with experimental values determined by neutron diffraction and infrared spectroscopy (see Refs. [30,27]) ˚) ˚) Method/basis set ROO (A ROH (A mOH (cm1) HF/3-21G HF/6-31G HF/6-31G(d) and beyond – up to the HF/6-311++G(2d,2p) B3LYP/3-21G B3LYP/6-31G B3LYP/6-31G(d) B3LYP/6-31G(d,p) B3LYP/6-31 + G(d,p) and beyond – up to the B3LYP/6-311++G(2d,2p) Experiment

2.445 2.534 Not stable 2.502 2.548 2.816 2.741 Not stable 2.432

1.141 1.056

1292 2971

1.251 1.097 1.014 1.738

1152 1893 3065 2960

1.156

(Very broad) 500–2000

K. Demsˇar et al. / Journal of Molecular Structure 844–845 (2007) 215–221

219

Table 2 Optimized values of selected geometric hydrogen bond parameters and harmonic frequencies of sodium hydrogenbissulfate in the solid state, calculated by using the BLYP density functional and different plane-wave basis sets, convergence criteria and crystal symmetry constraints ˚) ˚) Symmetry Cutoff Convergence ROO (A ROH (A mOH (cm1) Number and range (i cm1) (Ry)

criteriona

P21/c 90 Default 90 Tight P21/c P21/c 120 Default P1 90 Default Experimental (neutron diffraction) 1D anharmonic calculation along ROH

of imaginary modes 2.442 1.094 1847, 1852, 1857, 1897 2.439 1.094 1881, 1882, 1901, 1927 2.442 1.095 1846, 1876, 1905, 1987 2.443, 2.441, 2.441, 2.441 1.095, 1.093, 1.093, 1.094 1852, 1874, 1888, 1912 2.432 1.156 Very broad 500–2000 – 1.151 1571

9 at 1000; 4 at 200 9 at 1000 9 at 400; 4 at 100 0 – –

Experimental data and results of 1D anharmonic calculation based on the proton potential in the solid state (see text and Fig. 3) are also added. a Default convergence criterion used with the CPMD program: 5 · 104 a.u. for the maximum force component; tight convergence criterion: 5 · 105 a.u.

are questionable. Tightening of the geometry convergence criterion by a factor of 10 improves the calculation in that the lower 4 imaginary modes vanish; however the imaginary modes associated with the sodium ions located at special positions persist. Since these ions were not allowed to move during geometry optimization due to symmetry constraints, one may find the cause for inconsistent geometry optimization in the assumed symmetry. For that sake we performed geometry optimization in the absence of any symmetry within the unit cell, i.e., the P1 space group was applied, allowing the formerly fixed sodium atoms to move from their initial positions. The result was a structure free of imaginary modes, suggesting that the experimentally observed crystal symmetry needs to be broken in order to obtain a globally optimized geometry. Since the P21/c symmetry constraints prevent some of the sodium ions from leaving their positions, geometry optimization under P21/c can be regarded as constrained. In contrast, when all the crystallographic symmetry constraints are omitted and P1 is used, this corresponds to a full geometry optimization. Despite the fact that the discussed sodium ions were allowed to move freely in the P1 optimization course, their displacements from original positions were of the order of ˚ . Additionally, no notable inequiva few hundreths of an A alence in the geometry of the four hydrogen bond moieties was introduced because of the symmetry lowering. We can conclude that the applied pseudopotentials in conjunction with plane-wave basis sets of different size do not reproduce perfectly the experimental structure, giving rise to imaginary vibrational modes. Despite that the mismatch of the geometries was found to be marginal, the vibrational analysis is sensitive enough to distinguish between the P21/ c and P1 optimized geometries. The fact that the imaginary frequencies obtained with the larger plane-wave basis set of 120 Ry were notably smaller in magnitude than those obtained with the smaller basis set of 90 Ry (see Table 2) suggests that the assumed P21/c structure is closer to the global P1 minimum when the larger basis set is used. Indeed, comparison of displacements of the sodium ions from their initial special positions during the course of a P1 optimization with different plane-wave cutoffs shows that with the larger cutoff of 120 Ry the average displace-

˚ while with the smaller cutoff the average ment is 0.025 A ˚ . This suggests that with the displacement is 0.052 A increasing basis set the fully optimized structure approaches the experimentally determined P21/c symmetry, which confirms the reliability of the applied model. We found that the values of the harmonic OAH stretching frequencies show minor difference between the applied calculation settings (see Fig. 3 and Table 2). In any case, a set of four different values was obtained for the mOH modes. Note that the four frequencies are not equal, even when the P21/c symmetry is imposed. This can be interpreted by the crystal field splitting effect. While the splitting effect is usually very small for most of the modes, it is of the order of a few tens of wavenumbers for the OAH stretching mode. The calculated harmonic stretching frequencies vary slightly with the applied level of calculation, but the difference is not significant, particularly in the sense that the harmonic approximation is known to be insufficient for the description of the proton motion in short hydrogen bonds. Indeed, the calculated values of about 1900 cm1 are at the higher end of the extremely broad absorption observed for the hydrogen bond of a HBS3 ion in the solid state [27]. The absorption attributed to the OAH stretching mode extends from about 2000 cm1 to the low frequency limit of recording at 500 cm1. Thus the calculated harmonic frequencies do not match well with the experimental band center. An advanced model is therefore needed to reasonably reproduce the observed red-shifted OAH stretching frequency. For that sake we constructed the one-dimensional proton stretching potential energy function and numerically solved the vibrational time-independent Schro¨dinger equation corresponding to that potential. The energy gap between the first excited and the ground vibrational level is an estimate for the fundamental excitation frequency of the OAH stretching vibration. The ground level was found to be 2.34 kcal/mol above the potential energy minimum while the first excited state was at 6.83 kcal/mol above the minimum. This corresponds to the fundamental OAH stretching frequency of 1571 cm1 (see Fig. 3 and Table 2). This value, although probably still higher than the assumed center of the observed broad absorption in the infrared spectrum (tentatively located at around 1000 cm1 by visual

220

K. Demsˇar et al. / Journal of Molecular Structure 844–845 (2007) 215–221

estimation), is about 300 cm1 lower than the frequency calculated in the harmonic approximation and thus qualitatively in a better agreement with the experiment. From Table 2 it is evident that the predicted OAH dis˚ is considerably shorter than the one tance of 1.094 A ˚ ). However, expecobtained by neutron diffraction (1.156 A tation value of the OAH distance, calculated from the ground state vibrational wavefunction, was equal to ˚ , providing excellent agreement with the experimen1.151 A tal value. This is not surprising since the proton potential, from which the wavefunction was obtained, is essentially asymmetric and has a relatively wide low energy range at ˚ , allowing the proton wavefunction to be ROH > 1.094 A delocalized to longer OAH distances, thus enlarging the expectation value of ROH. We conclude that quantization of the nuclear motion for the highly anharmonic degree of freedom improves the results. In this study we considered periodicity and we solved the vibrational Schro¨dinger equation for the OAH motion beyond the harmonic approximation. Relative to the harmonic approximation a good agreement with the OAH bond length and an improved agreement with the OAH stretching frequency was found. Weak point of this study is that we did not reproduce the experimentally observed broad stretching band but a single line. Obviously it is necessary to introduce to the computational scheme thermal fluctuations and quantization of additional degrees of freedom. Thermal fluctuations can be easily included by conventional Car-Parrinello molecular dynamics simulation followed by a posteriori quantization of the proton motion [24], but this is not practical for this size of the unit cell. Quantization of additional degrees of freedom is more demanding since for example the O. . .O motion involves degrees of freedom that include several unit cells; moreover, the number of energy evaluations required to construct the potential energy surface increases exponentially with the number of considered degrees of freedom. It is worthy to stress that the methodology applied in this study is not limited to hydrogen-bonded systems in solid phase but can be applied to hydrogen bonds in solution and in enzyme active centers. 4. Conclusions We have calculated structural and vibrational properties of a short hydrogen bond of the hydrogenbissulfate anion using different environment models (gas phase vs. crystal field), various levels of theory and advanced methodology for treating the proton nuclear quantum effects. The high-level gas-phase models predict the hydrogenbissulfate ion to be unstable (a stable ion is obtained only with small and incomplete basis sets due to the basis set superposition error). In contrast, periodic plane-wave models proved to reliably reproduce the structure of the solid sodium hydrogenbissulfate, particularly the O. . .O distance. Due to the fact that the potential energy function of the hydrogen-bonded proton considerably deviates from harmonici-

ty, computational treatment beyond the harmonic approximation is essential for a reliable description of vibrational dynamics of the hydrogenbissulfate ion. We have shown that an approach based on the solving of a one-dimensional vibrational Schro¨dinger equation for the single point proton potential gives a notable improvement over the plain harmonic treatment in that the calculated OAH stretching frequency is shifted to a lower value and the OAH distance is enlarged, bringing the calculated results closer to experimental findings. Excellent agreement between the experimental and quantum-averaged OAH distance was found. We are aware that extensions of the treatment presented in this work are possible. The quantum vibrational problem can be extended to several, fully coupled degrees of freedom, giving rise to a better description of the vibrational dynamics in the hydrogen bond. Moreover, dynamical effects that give rise to band broadening can be modeled by inclusion of the thermal averaging. The presented methodology is also applicable to systems in the solid state as well as to solutions or to systems in fluctuating macromolecular environment, thus being interesting for the application to catalytic materials and biocatalysis. Acknowledgements Financial support from the Slovenian Ministry of Higher Education, Science and Technology (Grant P1-0012) is gratefully acknowledged. This work has also benefited from the use of facilities at the Center of Non-Linear Studies at Theory Division and the Lujan Neutron Scattering Center at LANSCE, which are funded by the Department of Energy’s Office of Basic Energy Sciences. Los Alamos National Laboratory is operated by Los Alamos National Security LLC under DOE Contract DE AC52-06NA25396. References [1] P. Colomban, Proton Conductors. Solids, Membranes and Gels – Materials and Devices, Cambridge University Press, 1992. [2] R. Blinc, R. Pirc, Collective behavior of hydrogen bonds in ferroelectrics and proton glasses, in: D. Hadzˇi (Ed.), Theoretical Treatments of Hydrogen Bonding, Wiley, Weinheim, 1997, pp. 229– 264. [3] G.R. Desiraju, Crystal Design: Structure and Function, Wiley, 2003. [4] A. Warshel, Computer Modelling of Chemical Reactions in Enzymes and Solutions, Wiley, New York, 1991. ˚ qvist, A. Warshel, J. Mol. Biol. 224 (1992) 7. [5] J. A [6] W.W. Cleland, M.M. Kreevoy, Science 264 (1994) 1887. [7] M.A. Massiah, C. Viragh, P.M. Reddy, I.M. Kovach, J. Johnson, T.L. Rosenberry, A.S. Mildvan, Biochemistry 40 (2001) 5682. [8] O. Lehtonen, J. Hartikainen, K. Rissanen, O. Ikkala, L.-O. Pietila¨, J. Chem. Phys. 116 (2002) 2417. [9] C.C. Wilson, C.A. Morrison, Chem. Phys. Lett. 362 (2002) 85. [10] J. Panek, J. Stare, D. Hadzˇi, J. Phys. Chem. A 108 (2004) 7417. [11] G.S. Denisov, J. Mavri, L. Sobczyk, Potential energy shape for the proton motion in hydrogen bonds reflected in infrared and NMR spectra, in: S.J. Grabowski (Ed.), Hydrogen Bonding – New Insights, Springer, 2006, pp. 377–416. [12] J.M. Bowman, J. Chem. Phys. 68 (1978) 608.

K. Demsˇar et al. / Journal of Molecular Structure 844–845 (2007) 215–221 [13] G. Chaban, J.O. Jung, B. Gerber, J. Chem. Phys. 111 (1999) 1823. [14] S. Webb, T. Iordanov, S. Hammes-Schiffer, J. Chem. Phys. 117 (2002) 4106. [15] G.G. Balint-Kurti, R.N. Dixon, C.C. Marston, Int. Rev. Phys. Chem. 11 (1992) 317. [16] J. Stare, J. Mavri, Comp. Phys. Comm. 143 (2002) 222. [17] D.A. Clabo, W.D. Allen, R.B. Remington, Y. Yamaguchi, H.F. Schaefer III, Chem. Phys. 123 (1988) 187. [18] N. Dosˇlic´, O. Ku¨hn, Z. Phys. Chem. 217 (2003) 1507. [19] P. Blaise, M.J. Wojcik, O. Henri-Rousseau, J. Chem. Phys. 122 (2005) 064306. [20] M. Boczar, L. Boda, M.J. Wojcik, J. Chem. Phys. 124 (2006) 084306. [21] M. Boczar, L. Boda, M.J. Wojcik, J. Chem. Phys. 125 (2006) 084709. [22] J. Mavri, J. Grdadolnik, J. Phys. Chem. A. 105 (2001) 2039. [23] J. Mavri, J. Grdadolnik, J. Phys. Chem. A. 105 (2001) 2045. [24] J. Stare, J. Panek, J. Eckert, J. Grdadolnik, J. Mavri, D. Hadzˇi (submitted for publication). [25] K. Gesi, J. Phys. Soc. Jpn. 61 (1992) 162. [26] J.M. Robertson, A.R. Ubbelohde, Proc. Roy. Soc. A 170 (1939) 222. [27] U. Mikac, D. Hadzˇi, R. Blinc, Ferroelectrics 239 (2000) 375. [28] H. Ratajczak, A.M. Yaremko, Chem. Phys. Lett. 243 (1995) 348. [29] M. Catti, G. Ferrari, G. Ivaldi, Acta Crystallogr. B 35 (1979) 525. [30] W. Joswig, H. Fuess, Acta Crystrallogr. B 38 (1982) 2798. [31] M.J. Frisch, G.W. Trucks, H.B. Schlegel, G.E. Scuseria, M.A. Robb, J.R. Cheeseman, J.A. Montgomery Jr., T. Vreven, K.N. Kudin, J.C.

[32] [33] [34] [35] [36] [37]

221

Burant, J.M. Millam, S.S. Iyengar, J. Tomasi, V. Barone, B. Mennucci, M. Cossi, G. Scalmani, N. Rega, G.A. Petersson, H. Nakatsuji, M. Hada, M. Ehara, K. Toyota, R. Fukuda, J. Hasegawa, M. Ishida, T. Nakajima, Y. Honda, O. Kitao, H. Nakai, M. Klene, X. Li, J.E. Knox, H.P. Hratchian, J.B. Cross, C. Adamo, J. Jaramillo, R. Gomperts, R.E. Stratmann, O. Yazyev, A.J. Austin, R. Cammi, C. Pomelli, J.W. Ochterski, P.Y. Ayala, K. Morokuma, G.A. Voth, P. Salvador, J.J. Dannenberg, V.G. Zakrzewski, S. Dapprich, A.D. Daniels, M.C. Strain, O. Farkas, D.K. Malick, A.D. Rabuck, K. Raghavachari, J.B. Foresman, J.V. Ortiz, Q. Cui, A.G. Baboul, S. Clifford, J. Cioslowski, B.B. Stefanov, G. Liu, A. Liashenko, P. Piskorz, I. Komaromi, R.L. Martin, D.J. Fox, T. Keith, M.A. AlLaham, C.Y. Peng, A. Nanayakkara, M. Challacombe, P.M.W. Gill, B. Johnson, W. Chen, M.W. Wong, C. Gonzalez, J.A. Pople, Gaussian 03, Revision C.02, Gaussian, Inc., Wallingford, CT (2004). C. Hartwigsen, S. Goedecker, J. Hutter, Phys. Rev. B58 (1998) 3641. Car-Parrinello Molecular Dynamics, version 3.9.2; Copyright IBM Corporation 1990–2004. D. Babic´, S.D. Bosanac, N. Dosˇlic´, Chem. Phys. Lett. 358 (2002) 337. J. Stare, G.G. Balint-Kurti, J. Phys. Chem. A 107 (2003) 7204. P. Pulay, Mol. Phys. 17 (1969) 197. J.G.C.M. van Duijneveldt-van de Rijdt, F.B. van Duijneveldt, Ab initio methods applied to hydrogen-bonded systems, in: D. Hadzˇi (Ed.), Theoretical Treatments of Hydrogen Bonding, Wiley, 1997, pp. 13–47.