Estimating the Hydrogen Bond Energy - The ... - ACS Publications

J. Phys. Chem. A 2010, 114, 9529–9536

9529

Estimating the Hydrogen Bond Energy Katharina Wendler, Jens Thar, Stefan Zahn, and Barbara Kirchner* Lehrstuhl fu¨r Theoretische Chemie, Wilhelm-Ostwald-Institut fu¨r Physikalische und Theoretische Chemie, UniVersita¨t Leipzig, Linne´str. 2, D-04103 Leipzig, Germany ReceiVed: April 18, 2010; ReVised Manuscript ReceiVed: July 26, 2010

First, different approaches to detect hydrogen bonds and to evaluate their energies are introduced newly or are extended. Supermolecular interaction energies of 256 dimers, each containing one single hydrogen bond, were correlated to various descriptors by a fit function depending both on the donor and acceptor atoms of the hydrogen bond. On the one hand, descriptors were orbital-based parameters as the two-center or threecenter shared electron number, products of ionization potentials and shared electron numbers, and the natural bond orbital interaction energy. On the other hand, integral descriptors examined were the acceptor-proton distance, the hydrogen bond angle, and the IR frequency shift of the donor-proton stretching vibration. Whereas an interaction energy dependence on 1/r3.8 was established, no correlation was found for the angle. Second, the fit functions are applied to hydrogen bonds in polypeptides, amino acid dimers, and water cluster, thus their reliability is demonstrated. Employing the fit functions to assign intramolecular hydrogen bond energies in polypeptides, several side chain CH · · · O and CH · · · N hydrogen bonds were detected on the fly. Also, the fit functions describe rather well intermolecular hydrogen bonds in amino acid dimers and the cooperativity of hydrogen bond energies in water clusters. 1. Introduction Hydrogen bonds are a very abundant motif in nature and are observed to be crucial in, for example, many biochemical processes.1 For example, they direct the DNA structure, protein folding, and substrate-enzyme interactions. Hydrogen bonds are used also to build supermolecular capsules or reversibly associated polymers as they enable self-assembly and molecular recognition.2 They determine the properties of water, which is the most important liquid for life and displays almost 40 anomalies.3 Thus, hydrogen bonds are important in many areas and a good knowledge of hydrogen bond energies can provide insight into fundamental chemical processes and help sophisticated drug and material design. Although hydrogen bonds are widely studied, the physical origin of hydrogen bonds is not fully revealed even after 100 years of research,4,5 which leads to the following short and wide definition: “The hydrogen bond is an attractive interaction between the hydrogen from a group XH and an atom or a group of atoms Y, in the same or different molecule(s), where there is evidence of bond formation”.6 Hydrogen bonds cover a broad range from very strong, having covalent character, to weak, having energies slightly above van der Waals interactions. Further, the hydrogen bond energy cannot be measured directly and, hence, its contribution to the total energy can only be estimated. Geometric criteria are often applied to detect hydrogen bonds and assign energies. But these approaches are error-prone7 and are said to be doomed to fail in complicated situations such as bifurcated or cooperative hydrogen bonds.3 Hence, a wave function-based criterion is desirable because it is supposedly sensitive to different environments, that is, to solvent effects and cooperative factors.8–10 Especially, it is able to identify unexpected and hidden hydrogen bonds as no arbitrarily cutoff criterion is needed. In this study, the perfor* To whom correspondence should be addressed. E-mail: bkirchner@ uni-leipzig.de.

mance of geometric and wave function-based approaches is compared, as well as in situations such as water clusters in which cooperativity plays a significant role. 2. Approach Numerous different approaches were proposed to estimate the hydrogen bond strength.11–23 For the quantitative description of hydrogen bond interactions, it is common practice to define an interaction energy that is identified with the hydrogen bond energy. In the supermolecular approach, the complex including hydrogen bonds is handled as one molecule and the total energy EAB is determined. If the total energies of the unrelaxed monomers EA and EB are taken as reference, the intrinsic interaction energy E is

bA + b bA) - EB(R bB) E ) EAB(R RB) - EA(R

(1)

RB are the coordinates of the atoms of monomers where b RA and b A and B in the complex. Due to the definition in eq 1, the energy of an attractive interaction is negative. For an easier reading, the term energy denotes the absolute value of the energy |E| in the following. Necessarily, this interaction energy captures also other van der Waals interaction as dispersion or induction as well as possible repulsive interactions, electrostatic contributions and others. Nevertheless, the intrinsic interaction energy is taken as an estimate of the hydrogen bond energy in this article. We extended the shared electron number approach introduced previously.24 The idea of this method is to estimate the strength of a hydrogen bond by means of only one variable, the twocenter shared electron number. This ansatz has the main advantage that it can be applied to complex systems since the shared electron number is quantum chemically available even for large systems. Qualitatively, bond strengths are often associated with the number of bond electrons in a Lewis structure. The charge-

10.1021/jp103470e  2010 American Chemical Society Published on Web 08/13/2010

9530

J. Phys. Chem. A, Vol. 114, No. 35, 2010

Wendler et al.

density analysis by Davidson,25 Roby,26 and Heinzmann and Ahlrichs27 defines a two-center shared electron number σ (SEN) that provides a quantitative measure for the strength of a covalent bond.27,28 Ehrhardt and Ahlrichs28 proposed a “simplest possible” relationship between the two-center SEN σ and the covalent bond energy EAB:

EAB ) a

IPA + IPB σ+m 2

(2)

in which IPX is the ionization potential of atom X, and a and m are constants, the latter shall be close to zero. Reiher, Sellmann, and Hess24 first applied the SEN approach to intermolecular, noncovalent hydrogen bonds. They used a linear relationship:

Eσ ) aσ

(3)

in which Eσ is the hydrogen bond energy and σ is the twocenter SEN. a is a proportionality constant that was found to be between 380 kJ/(mol · e) and 637 kJ/(mol · e) depending on the method.24 The observed standard deviation was less than 3.5 kJ/mol. A new parametrization allowing for an intercept was introduced in our group recently.29 The fit improved further if the calibration set was subdivided by donor and acceptor atoms into subsets for which the parametrization was done.29 Thus, the SEN approach was applied, for example, to transitionmetal compounds,24,30,31 tetracycline-earth metal complexes,32 peptides,33 host-guest interactions in an artificial adrenaline receptor complex,29 enzymes,34 DNA and RNA,35 bipyridine derivatives,36 phenanthroline,37 ionic liquids,38,39 and rotaxanes.40 In this article, we examined the performance of various descriptors, particularly the two- and three-center SEN, in the general fit functions without and with an intercept m:

E ) wdon(adonx) + wacc(aaccx)

(4)

E ) wdon(adonx + mdon) + wacc(aaccx + macc)

(5)

where the index “don” refers to the donor subset, the index “acc” is the acceptor subset, a is the slope, and m is the intercept of a linear fit. w is a weighing parameter, and E equals the intrinsic supermolecular interaction energy. The linear fit was done in two steps: First, a and m were fitted separately for donor and acceptor subsets, afterward w were fitted using the obtained parameters for a and m. In another, the natural bond orbital (NBO) approach, a hydrogen bond is viewed as an interaction between an occupied nonbonded natural orbital of the acceptor atom nA and the unoccupied antibonding orbital of the DH bond σ*DH. The NBO interaction E2 is estimated by

E2 ) ∆Eij ) qi

F(i, j)2 εj - εi

(6)

in which qi is the donor orbital occupancy, εi and εj are diagonal elements of the Fock matrix and, thus, orbital energies, and F(i, j) is the off-diagonal NBO Fock matrix element.41 This second-order estimate was found to be 1 or 2 orders of magnitude greater for hydrogen bonds than for any other single nA f σ*DH contribution in hydrogen-bonded complexes.42 Hence,

it is used to estimate the hydrogen bond energy43 and was tested as descriptor x in eqs 4 and 5 in this article. As angle and distance are well accessible, the standard approach for the estimation of local interaction energies in complex aggregates is based on geometric criteria.3 These define the interaction of two fragments of an aggregate solely on the basis of distances and, occasionally, of angles.17,44 Also, the frequency shift of the DH stretching vibration is a widely used descriptor for the hydrogen bond strength.45–49 Indeed, linear correlations between the energy and the frequency shift were found for the red-shifted NH stretching vibration and the blueshifted CH stretching vibration in rotaxane mimetics.40,50 3. Computational Details Isolated dimer complexes were created by the combination of 65 acceptors and 28 donors. The peptides structures were taken from Deshmukh et al.51 All dimers geometries, the peptides and the water clusters were optimized using DFT theory52 with the dispersion-corrected BLYP-D53–55 functional as implemented in TURBOMOLE.56,57 The calculations were done in combination with the density fitting technique58,59 and the basis set def2-TZVPP was used for all atoms. It was established that the structures belonged to a energy minimum by frequency analysis. Furthermore, dimers fulfilling the SEN criterion (see below) were optimized applying second-order Møller-Plesset perturbation theory (MP2)60,61 in TURBOMOLE. The density fitting technique and the basis set def2TZVPP with frozen core were used in the MP2 calculations as well. For the determination of the harmonic frequencies, the SNF 4.0 program package62 was employed. The NBO analysis was performed with NBO 5.042,63 augmented in the GAUSSIAN 03d64 version. For reference examples, single point coupled cluster calculations (CCSD(T)) were carried out for several MP2 optimized dimers in MOLPRO65 with an cc-aug-pVQZ66,67 basis set for all atoms. The dimer interaction energies were calculated by applying eq 1. The final interaction energies were counterpoise corrected with the procedure by Boys and Bernadi.68 For calibration, dimers were selected in which the interaction energy was dominated by one hydrogen bond and contributions from interactions between other atoms were negligible. The arbitrarily chosen criterion stated that any second two-center SEN contact between the two monomers was less than 5% of the two-center SEN contact between the donor proton and the acceptor atom. All other dimers featuring more than one SEN contact were excluded from the analysis. Hence, the calibration set consisted of 256 BLYP-D optimized dimers and 181 MP2 optimized dimers. All had an attractive interplay. In general, acceptor denotes proton-acceptor and donor refers to proton-donor. The root-mean-square deviation s (rmsd) of the fit results was defined here as:

s)




n

∑ [Efit(i) - E(i)]2 i)1

n

(7)

in which Efit(i) denotes the fitted energy of dimer i with the found interaction energy E(i) and n is the number of all dimers in the calibration set. The correlation coefficient R between the found interaction energies and the fit energies was calculated as it is defined in general. 4. Results and Discussion 4.1. Reliability of Interaction Energies. As CCSD(T) was the most reliable method used here, the CCSD(T) interaction

Estimating the Hydrogen Bond Energy

J. Phys. Chem. A, Vol. 114, No. 35, 2010 9531

TABLE 1: Calculated Interaction Energies in kJ/mol BLYP-Da dimer FH · · · FH ClH · · · ClH ClH · · · FH OH2 · · · OH2 OH2 · · · SH2 a

MP2a

SCS-MP2a

CCSD(T)b

def2-TZVPP def2-TZVPP def2-TZVPP cc-aug-pVQZ -21.1 -8.37 -12.2 -21.5 -12.9

-16.6 -6.98 -10.3 -19.1 -10.7

-15.3 -5.03 -8.84 -17.5 -8.94

-18.9 -7.94 -11.3 -20.6 -11.9

Using TURBOMOLE. b Using MOLPRO.

TABLE 2: Linear Fit Parameters of BLYP-D Interaction Energy to Two-Center Shared Electron Numbera wdon

wacc

s

R

w int don

int wacc

sint

Rint

0.706

0.308

5.17

0.9659

0.746

0.265

4.69

0.9705

donor

adon

a int don

mdon

acceptor

aacc

int aacc

macc

C Si N P O

-832 -262 -514 -383 -429

-845 -389 -521 -90.4 -390

0.116 0.461 0.123 -2.73 -4.09

S F Cl Br

-309 -443 -305 -279

-244 -401 -244 -224

-3.16 -4.82 -8.47 -7.17

N P O Oen S F Cl Br

-335 -302 -437 -418 -298 -641 -457 -406

-302 -278 -393 -369 -276 -598 -387 -352

-5.33 -2.06 -3.21 -4.72 -2.07 -0.867 -1.48 -1.36

a σ in e (a in kJ/(mol · e), m in kJ/mol, and RMSD (s) in kJ/mol; don: donor group; acc: acceptor group; ax, wx: fit without intercept int (eq 4); aint x , wx : fit with intercept (eq 5)).

Figure 1. BLYP-D interaction energy fits to (a) the two-center shared electron number σ and (b) the product of ionization potential IPdon and σ.

energies might be taken as reference. Table 1 gives the interaction energies for the five benchmark dimers calculated by different methods. For water, a CCSD(T) interaction energy of 20.8 kJ/mol was reported,69 for the HF dimer one of about 19.1 kJ/mol70 to 19.4 kJ/mol was found.69 With another dispersion correction, a BLYP interaction energy of 21.7 kJ/ mol is obtained for water.71 As shown in Table 1, the trends were the same for all five dimers. Whereas the BLYP-D interaction energy was larger than the CCSD(T) interaction energy, the MP2 interaction energy was lower. Thus, BLYP-D was overestimating the attractive interactions. SCS-MP2 gave an even lower interaction energy than MP2 as observed also by Bachorz et al.72 and, therefore, seemed less favorable to describe these interaction energies. 4.2. Correlations and Fit Parameters. 4.2.1. WaWe Function-Based Descriptors. In Figure 1a, the fit results applying eq 4 are shown. Mostly, the deviations between fit energies (i.e., energies obtained by the fit functions) and interaction energies

calculated were less than 10 kJ/mol. The parameters of the fit functions (see Table 2) were similar to the parameters reported before.24,35 For larger SEN, the deviations became higher. This might be caused by a general deficiency of SEN to describe the short-range repulsive interaction and, thus, overestimating the attractive interaction for short hydrogen bonds. Also, in some dimers, the proton is actually transferred from the donor to the acceptor atom. Thus, a covalent bond is formed. Each fit function yielded adon and aacc for every donor or acceptor subset. Additionaly, if oxygen is the acceptor atom, it is differed between subgroup O in which oxygen is next to hydrogen or carbon and the subgroup Oen in which oxygen is bonded to an electronegative atom like nitrogen or halogenes. These subset specific parameters a obeyed a clear trend: In the periodic table groups, both adon and aacc decreased. This may imply a connection to the electronegativity or the ionization energy. Thus, following Ahlrichs’ idea for covalent bonds (see eq 2), the ionization potentials were included in the SEN fits for the first time for noncovalent bonds. Concurrently, the devision in subsets was lifted and, thus, the number of fit parameters was decreased drastically. The ionization potentials were taken from NBO analysis. On the one hand, the donor ionization potential was chosen to be the one of the binding orbital σDH between the hydrogen atom and the donor atom. On the other hand, the acceptor was characterized by the ionization potential of the lone pair orbital nA interacting most with the orbital σDH as expressed in the NBO interaction energy E2. It was found that the fit was not improved significantly by including the acceptor ionization potential.57 Hence, the following fit functions are further applied:

EIPdon ) aIPdonσ ) 598kJ/(mol · au · e)IPdonσ (s ) 6.18 kJ/mol, R ) 0.9578)

(8) EIPdon+m ) aIPdonσ + m ) 547kJ/(mol · au · e)IPdonσ - 3.51 kJ/mol (s ) 5.59kJ/mol, R ) 0.9578)

(9) The simplest fit (eq 8) with just one fit parameter gave an rmsd that was only by 1 kJ/mol higher than the one of Eσ,

9532

J. Phys. Chem. A, Vol. 114, No. 35, 2010

Wendler et al.

TABLE 3: Variation of the Intercept Inclusion in the Sum of Two-Center SEN Based BLYP-D Fit Energies (in kJ/mol) for Several Dimers with Two SEN Contactsa dimer FH · · · FH ClH · · · FH OH2 · · · OH2

Eσ -18.3 -9.34 -18.5

Eσ/σ+m -20.3 -13.2 -20.6

Eσ+m -24.1 -17.2 -24.5

EBLYP-D -21.1 -12.2 -21.5

a Eσ, see eq 4; Eσ+m, see eq 5. Eσ/σ+m equals Eσ+m for the largest SEN contact and equals Eσ+m - m for all other SEN contacts.

compare Figure 1b and Table 2. Compared to the Eσ fit, 16 degrees of freedom were left out. Including an intercept (see fit eq 9), the fit improved slightly. Hence, the inclusion of ionization potentials gives a new powerful variant of the SEN approach, needing only one or two fit parameters. Ahlrichs et al.28 proposed that the intercept should be close to zero for intramolecular bonds. Covalent bonds have very high bond energies and high SEN, dispersion has no significant influence. Intermolecular, noncovalent, bonds are weaker and, thus, van der Waals interactions become important. Even if the dimer does not form hydrogen bonds, there is an attractive van der Waals interaction between the two molecules that will be denoted as “pure van der Waals interaction” (i.e., without the hydrogen bond component). Thus, one may interpret the intercept in fit function (eq 5) as a pure van der Waals term. For the halogens, the intercept increased in the periodic table groups. Higher elements have more diffusive electron densities. Thus, dispersion and induction and, therefore, the pure van der Waals interactions are stronger. However, the intercepts did not increase in other periodic table groups. This might be caused by the collection of calibration dimers or it might hint at the limits of the proposed interpretation. The pure van der Waals interactions should not alter much if a second weak hydrogen bond is formed in a dimer. Thus, if the hypothesis is correct and the intercept is connected to pure van der Waals interactions the intercept should be included only once for multiple intermolecular hydrogen bonds. All three dimers in table 3 have two SEN contacts, a big and a very small one. Eσ/σ+m was calculated with the Eσ+m fit parameters and the intercept was omitted for all but the largest SEN contact. Indeed, fit energies in which the intercept was included only once estimate the interaction energy best. Thus, the energy without the intercept should estimate the hydrogen bond energy without the contribution of pure van der Waals interactions. Employing fit eqs 4 and 5, the three-center SEN σ3 and the NBO interaction energy E2 gave as good correlation to the BLYP-D interaction energy. Due to the correlation with E2, one might conclude that the hydrogen bond might be described well by charge-transfer interactions. The correlation to MP2 interaction energies of the descriptors discussed here was as good as to the BLYP-D energies.57 4.2.2. Geometric Descriptors. The BLYP-D and MP2 interaction energies did not depend linearly on the distance r, but on the reciprocal distance 1/r3.8 (see Figure 2a). The obtained fits were of similar quality as the SEN fits (compare Figure 1 and Table 2).57 Probably, different interactions terms gave an overall decrease of 1/r3.8, for example, a combination of electrostatic and the very fast descending charge-transfer, both of which are discussed as major constituents of the hydrogen bond energy.73,44 The fit parameters are given in Tab. 4. As can be seen, the parameter aacc increased with respect to the periodic table group. For example, the data points of dimers with sulfur or bromine accecptor atoms are found clearly right of the fit curve in Figure

Figure 2. BLYP-D interaction energies in dependence of a (a) distance rH · · · A and (b) angle RDH · · · A (Colors correspond to different donor subsets).

TABLE 4: Linear Fit Parameters of BLYP-D Interaction 3.8 a 3.8 Energy to Distance 1/r H · · · A in 1/pm wdon

wacc

s

R

wint don

int wacc

sint

0.303

0.709

5.46

0.9664

0.138

0.869

4.44

adon

aint don

donor

mdon

C Si N P O

-9.70 -12.7 2.14 -1.26 -1.34 0.0413 -11.3 -10.2 -1.10 -5.33 -5.43 0.0727 -16.4 -16.9 1.23

S F Cl Br

-11.1 -14.9 -13.1 -12.0

-10.5 -13.7 -12.4 -10.2

-0.842 -4.22 -2.06 -5.05

Rint 0.9736

acceptor

aacc

int aacc

N P O Oen S F Cl Br

-15.6 -24.8 -11.2 -12.2 -22.0 -6.58 -12.3 -15.6

-18.4 -35.5 -13.3 -15.3 -30.3 -10.1 -18.1 -23.8

macc

8.62 9.39 4.96 9.11 8.87 4.78 3.89 4.70

a in 109 kJ/(mol · pm3.8), m in kJ/mol, and RMSD s in kJ/mol; don: donor group; acc: acceptor group; ax, wx: fit without intercept int (eq 4); aint x , wx : fit with intercept (eq 5). a

2a. These atoms have larger atomic radii than their congeners. Hence, the hydrogen bond distance, taken from nucleus to nucleus, might be increased without lowering the interaction energy. The acceptor (and its radius) was more important for the fit than the donor as was expressed in a wacc of more than 0.7. The donor parameters adon decreased in the group and, thus, reflected trends found also in the SEN based fits. The BLYP-D interaction energies are illustrated in dependence of the angle in Figure 2b. The colors were assigned to the donor subsets. Obviously, there was no correlation between

Estimating the Hydrogen Bond Energy

J. Phys. Chem. A, Vol. 114, No. 35, 2010 9533

∆r TABLE 5: Linear Fit Parameters of Frequency Shift ∆νfit to the Distance Change ∆rDH in pm (eq 10) (adon in 1019/m2) and Linear Fit Parameters of BLYP-D Interaction Energy to Frequency Shift ∆ν of the Hydrogen-Donor Stretching Vibration νDH ∆ν in cm-1 (eq 11) (adon in kJ · cm/mol)

donor

C

Si

N

P

O

S

F

Cl

Br

a∆r don a∆ν don

-10.7 0.130

-5.17 -0.00799

-10.6 0.0529

-7.05 -0.0453

-9.21 0.0705

-13.0 0.0727

-20.4 0.0719

-12.2 0.0631

-10.1 0.0568

energy and angle in any subset. A coarse rule of thumb might be apparent: the smaller the angle is, the smaller the interaction energy is. Still, it was not tested if a correlation existed for special system like, for exmaple, water or proteins, as was discussed in literature.18,74 However, one should note that almost all hydrogen bonds had an angle RDH · · · A above 150° in accordance with the angle cutoff criteria in many molecular dynamics simulations.75 4.2.3. Frequency Shift. The band shift was detected for the hydrogen-donor stretching vibration νDH. It was calculated by: ∆νDH ) νDH,dimer - νDH,donor. The frequencies of the donor molecule were calculated for the donor in an optimized geometry. All frequencies were obtained using SNF 4.0. In Table 5, the fit parameters are given according to:

∆νfit ) a∆r with don∆r ∆r ) rDH,dimer - rDH,donor

(10)

E∆ν ) a∆ν don∆ν

(11)

and

The results of the first fit (eq 10) are illustrated in Figure 3a, the rmsd was 66.1 cm-1 and R was 0.9834. As can be seen, the linearity was kept rather well. The data points with large deviations in Figure 3a correspond to dimers with oxygen as donor atom. Probably, the acceptor atom had a strong influence for oxygen donor dimers. Red-shifted and blue-shifted hydrogen bonds behaved similar. Thus, they can be described with the same fit parameters. Most of the deviations between BLYP-D interaction energy and fit energy (see eq 11) were about 7 kJ/mol (see Figure 3b), the rmsd was 6.57 kJ/mol and R was 0.9614. A good correlation was found for some donor subsets such as, for example, oxygen or nitrogen donor dimers, which showed high frequency shifts. Subsets with donor atoms participating only in weak hydrogen bonds such as, for example, carbon did not show correlation. Indeed, the stretching vibrations νDH associated with weak hydrogen bonds were often coupled to other vibrations, especially for CH vibrations. Negative slopes (see Table 5) and, thus, overall blue-shift trends were found for silicon and phosphorus donor subsets that were already reported by Schlegel et al.76 If the dimers showing blue-shifts or red-shifts were fitted separately, the overall rmsd decreased to 5.98 kJ/mol.57 Thus, this rmsd was slightly larger than the rmsd of the two-center SEN based fit without intercept, which was 5.17 kJ/mol (see table 2). It is very likely not sufficient to consider only harmonic frequencies and neglect anharmonicity because the corresponding errors of the frequencies may be too large for a reasonable fit. The correlation of the change in absolute intensity to the frequency shift gave a very large rmsd of 7.94 kJ/mol.57 However, the intensities themselves were connected to a significant uncertainty.

Figure 4. BLYP-D geometrically optimized structure of alanine-alanine dimer and obtained Eσ fit energies in kJ/mol (blue: N, orange: C, red: O, white: H).57

4.3. Intermolecular Hydrogen Bonds in Amino Acid Dimers. The reliability and the applicability of all fit functions were tested on the supermolecular interaction energies of four amino acid dimers. The results were similar in each case and will be discussed for alanine-alanine in the following.57 The sum of all hydrogen bond energies Eσ of -24.5 kJ/mol, as denoted in Figure 4, was 36% smaller than the BLYP-D interaction energy of -38.2 kJ/mol. But, even Eσ+m, which included intercepts or any other wave function based fit, underestimated the dimer interaction energy. In contrast, E1/r3.8

Figure 3. (a) Frequency shift ∆νDH dependence on distance change ∆rDH, and (b) BLYP-D interaction energy dependence on the frequency shift ∆νDH.

9534

J. Phys. Chem. A, Vol. 114, No. 35, 2010

overestimated the interaction energy by 4.2 kJ/mol even if the CH · · · O hydrogen bonds were not regarded. In contrast to the distance-based fit, the SEN based fits enabled the detection of hydrogen bonds without any cutoff criteria. So, two CH · · · O contacts were found in alanine-alanine. However, there were other SEN contacts that are not shown in Figure 4, for example, between the two NH groups or between carbon and oxygen atoms. These interactions were not parametrized and, thus, it was not possible to assign an interaction energy. In particular, the NH · · · NH interaction is classified as a secondary repulsive interaction.44 The SEN method failed to describe these repulsive interactions as the NH · · · NH interaction was characterized by a positive two-center SEN just like attractive interactions. The other orbital-based descriptors as the NBO interaction energy E2 or the three-center SEN σ3 failed to distinguish between the attractive and repulsive interactions, too. E1/r3.8 cannot display repulsive interactions due to the design of the distance fit function. Additionally, the pure van der Waals interaction may be different from the situation in the small dimers used for calibration. Thus, the deviation of the predicted interaction energy might be understandable. The frequency shift fits were not applied because the vibrations were strongly coupled. Furthermore, the assignment of any frequency shift was not possible on this level of theory as both NH · · · N and NH · · · O hydrogen bonds occurred at the same NH group displaying naturally only one overall shift.77,78 Nevertheless, the hydrogen bond fit energies of all fit functions showed the general expected order4 of: ENH · · · O ≈ -6 kJ/mol > ENH · · · N ≈ ECH · · · O ≈ -3 kJ/mol. Hence, the fits introduced here seemed suitable for the estimation of hydrogen bond energies. 4.4. Intramolecular Hydrogen Bonds in Peptides. Intramolecular hydrogen bonds are very interesting, but they are not accessible directly from experiment. Hence, theoretical results are the only reference data. The nine examined peptide structures were taken from literature,51 optimized with the standard procedure used here and analyzed by each previous introduced fit function. The NH · · · N interaction were generally lower in energy than the NH · · · O hydrogen bonds.57 The NH · · · N Eσ fit energies were between -5.13 and -7.88 kJ/mol and the NH · · · O Eσ fit energies were between -6.62 and -14.1 kJ/mol.57 Thus, the NH · · · O fit energies were spread over a broader range than the NH · · · N fit energies. The energies are in the range proposed in literature.3,13,79–83 These fitted hydrogen bond energies were higher than for the amino acid dimers. This might show that the fit functions described cooperativity to a certain degree. Appealingly, the increase in the fit energies was the highest for EIPdon yielding 80% higher fit energies for the polypeptide than for the amino acid dimer. Also, Eσ (66%) grew more than E1/r3.8 (40%). Thus, one may conclude that the distance-based fit covered less of the energy increase as would be expected for phenomenon based on charge-flow like π cooperativity.84 In contrast to methods known in literature, the wave function based fits presented here can be used to detect hydrogen bonds in the whole peptide structure without any cutoff criteria. Indeed, there were SEN contacts between side chain protons and backbone nitrogen or oxygen atoms in five of the nine studied peptides57 that are illustrated in Figure 5. In general, the Eσ fit energies were in the range from -0.97 to -5.27 kJ/mol.57 The CH · · · O (-2.05 and -5.27 kJ/mol) interactions were stronger than the CH · · · N (-0.97 to -2.97 kJ/mol).57 This also corresponded to the distances of rCH · · · O, which was between 237

Wendler et al.

Figure 5. BLYP-D geometrically optimized structures of polypeptides with the detected side chain interactions (blue: N, orange: C, red: O, white: H): (a) AcGIGGGNH2, (b) AcGVGGGNH2, (c) AcAAAAANH2, (d) AcLLLLLNH2, (e) AcVVVVVNH2 (standard perspective), and (f) AcVVVVVNH2 (other perspective).

and 279 pm and, thus, smaller than rCH · · · N, which was 271-317 pm. Whereas CH · · · O interactions were discussed in literature,4,13,80,85–89 CH · · · N interactions were not described yet. The CH · · · N formed five- or six-membered rings, nitrogen interacted always with side chain β protons. The CH · · · O formed bigger six- or seven-membered rings, oxygen interacted with side chain β or γ protons. Anyway, these energies were lower than NH · · · O and NH · · · N hydrogen bond energies but they might influence the position of the side chains.88 The possible great advantage of the here introduced wave function based fit functions is that they enable the detection of hydrogen bonds without any cutoff criterion on the fly and in the whole system studied. Additionally, these fits do not depend on the relative energies of selected structures as other methods,13,16,90 but only on the system in an optimized geometry and its ground state wave function. So, the calculation of SEN or the NBO analyses are done easily. 4.5. Hydrogen Bond Energy Cooperativity in Water Clusters. The nonadditivity of hydrogen bond properties is caused by mutual polarization of the involved groups.44 Using QCE theory, it was found that the neglect of cooperativity leads to a much worse description of liquid water than disregarding dispersion.3,91 Thus, cooperativity is a very important aspect of hydrogen bonds. Five different water clusters, from trimers to a hexamer, were studied using the different fit functions, a selection is illustrated in Figure 6.57 Within the clusters, the hydrogen bond energies were fitted as described. Applying the supermolecular approach, the interaction energy of one molecule i to all other molecules was obtained. This interaction energy should equal the sum of the

Estimating the Hydrogen Bond Energy

J. Phys. Chem. A, Vol. 114, No. 35, 2010 9535 Acknowledgment. We would like to thank the RZ Leipzig for computational resources. J. T. thanks the Fonds der Chemischen Industrie for financial support. This work was supported by DFG, in particular by the projects KI-768/4-1, KI-768/5-1, and KI-786/5-2 of the priority program SPP 1191 “Ionic Liquids”.

Figure 6. BLYP-D geometrically optimized structures of water clusters with estimated hydrogen bond energies Eσ/σ+m in kJ/mol (red: O, white: H): (a) dimer, (b) trimer, (c) pentamer, and (d) hexamer.

fitted hydrogen bond energies from molecule i to all other molecules in the cluster. The Eσ/σ+m fit function predicted the absolute values of the BLYP-D interaction energies most exactly with an rmsd of 2 kJ/mol.57 For MP2 interaction energies, the rmsd was about 4 kJ/mol.57 The hydrogen bond energy was found to be twice as large in the equilateral hexamer as in the dimer whereas it is only sligthly increased in the trimer compared to the dimer (see Figure 6). The highest difference in hydrogen bond energies within one cluster was found in a pentamer (see Figure 6). As can be seen, the hydrogen bond energies vary by more than 17 kJ/mol. The hydrogen bond distance from the isolated water molecule to the molecule in the cluster center is shorter than the distance between the center molecule and the two twice-hydrogen bonded water molecules. Still, the isolated water molecule was predicted to form a weaker hydrogen bond (-20.4 kJ/mol) to the center molecule than the two others (-23.3 kJ/mol). Thus, this energetic order of Eσ/σ+m fit energies might be attributed to the electron distribution.57 As the molecule was isolated on the other side, the electron density was probably lower there than on the other molecules. In the trimer, one hydrogen bond had a slightly lower Eσ/σ+m fit energy (Figure 6) which was consistent with the calculated BLYP-D interaction energies. Applying the E1/r3.8 fit, it is found that another one of the three hydrogen bonds was stronger than the other two which were equally strong according to E1/r3.8. Thus, the distance-based fit gave the wrong energetic order of the hydrogen bonds. In all other cases, the trends predicted by the distance-based fit were the same as for the SEN-based fits. 5. Summary and Conclusion All presented fits were capable to estimate the interaction energy of the calibration dimers with an rmsd of -4 to -8 kJ/ mol and thus offer a simple approach to hydrogen bond energies. For several fit functions, fit parameters are given here. It is suggested that the shared-electron number fit intercepts were linked to the pure van der Waals interaction between the dimers. Extraordinarily, the product of the ionization potential of the hydrogen-donor binding orbital and the shared-electron number correlated well with the interaction energy. With the distance based E1/r3.8 fit, a very straightforward estimate was proposed for hydrogen bond energies which did perform similar to the wave function descriptor fits. Cooperativity was described rather well by all fits. However, the orbitalbased shared-electron number fits enabled the detection of hydrogen bonds of e.g. peptide side chains without any arbitrary cutoff criterion and in whole structure studied. Geometric criteria might never permit such assignments. For intramolecular hydrogen bonds in peptides, the fit energies were in agreement with literature values. Very alike results were found using the MP2 method instead of dispersion-corrected BLYP-DFT.

References and Notes (1) Jeffrey, G.; Saenger, W. Hydrogen Bonding in Biochemical Structures; Springer-Verlag: 1991. (2) Desiraju, G. R. Chem. Commun. 1997, 2, 1475–1482. (3) Kirchner, B. Phys. Rep. 2007, 440, 1–111. (4) Desiraju, G. R.; Steiner, T. The Weak Hydrogen Bond; Oxford University Press: 1999. (5) Goswani, M.; Arunan, E. Phys. Chem. Chem. Phys. 2009, 11, 8974– 8983. (6) Arunan, E. Curr. Sci. 2007, 92, 17–18. (7) Lakshmi, B.; Samuelson, A. G.; Jose, K. V. J.; Gadre, S. R.; Arunan, E. New. J. Chem. 2005, 29, 371–377. (8) Raghavendra, B.; Mandal, P. K.; Arunan, E. Phys. Chem. Chem. Phys. 2006, 8, 5276–5286. (9) Schwo¨bel, J.; Ebert, R. U.; Ku¨hne, R.; Schu¨u¨rmann, G. J. Phys. Chem. A 2009, 113, 10104–10112. (10) Mohajeri, A.; Nobandegani, F. F. J. Phys. Chem. A 2008, 112, 281– 295. (11) Iogansen, A. V. Spectrochim. Acta A 1999, 55, 1585–1612. (12) Kabsch, W.; Sander, C. Biopolymers 1983, 22, 2577–2637. (13) Zhang, Y.; Wang, C.-S. J. Comput. Chem. 2009, 30, 1251–1260. (14) Estacio, S. G.; Cabral do Couto, P.; Costa Cabral, B. J.; Minas da Piedade, M. E.; Martinho Simoes, J. A. J. Phys. Chem. A 2004, 108, 10834– 10843. (15) Schlund, S.; Schmuck, C.; Engels, B. J. Am. Chem. Soc. 2005, 127, 11115–11124. (16) Deshmukh, M. M.; Gadre, S. R.; Bartolotti, L. J. J. Phys. Chem. A 2006, 110, 12519–12523. (17) Luzar, A.; Chandler, D. Phys. ReV. Lett. 1996, 76, 928–931. (18) Lipsitz, R. S.; Sharma, Y.; Brooks, B. R.; Tjandra, N. J. Am. Chem. Soc. 2002, 124, 10621–10626. (19) Gu, Y.; Kar, T.; Scheiner, S. J. Am. Chem. Soc. 1999, 121, 9411– 9422. (20) Platts, J. A. Phys. Chem. Chem. Phys. 2000, 2, 973–980. (21) Platts, J.; Butina, D.; Abraham, M.; Hersey, A. J. Chem. Inf. Comp. Sci. 1999, 39, 835–845. (22) Grabowski, S. J. Chem. Phys. Lett. 2001, 338, 361–366. (23) Abraham, M. H.; Ibrahim, A.; Zissimos, A. M. J. Chrom. A 2004, 1037, 29–47. (24) Reiher, M.; Sellmann, D.; Hess, B. A. Theor. Chem. Acc. 2001, 106, 379–392. (25) Davidson, E. R. J. Chem. Phys. 1967, 46, 3320–3324. (26) Roby, K. R. Mol. Phys. 1974, 27, 81–104. (27) Heinzmann, R.; Ahlrichs, R. Theor. Chim. Acta 1976, 42, 33–45. (28) Ehrhardt, C.; Ahlrichs, R. Theor. Chim. Acta 1985, 68, 231–245. (29) Thar, J. Investigations of Hydrogen Bonds in Supermolecular HostGuest Complexes, Master’s thesis, Friedrich-Wilhelms-University Bonn: 2005. (30) Reiher, M.; Salomon, O.; Sellmann, D.; Hess, B. A. Chem.sEur. J. 2001, 7, 5195–5202. (31) Schenk, S.; Le Guennic, B.; Kirchner, B.; Reiher, M. Inorg. Chem. 2008, 47, 3634–3650. (32) Schneider, S.; Schmitt, M. O.; Brehm, G.; Reiher, M.; Matousek, P.; Towrie, M. Photochem. Photobiol. Sci. 2003, 2, 1107–1117. (33) Reiher, M.; Kirchner, B. J. Phys. Chem. A 2003, 107, 4141–4146. (34) Schmidt, M.; Zahn, S.; Carella, M.; Ohlenschla¨ger, O.; Go¨rlach, M.; Kothe, E.; Weston, J. ChemBioChem 2008, 9, 2135–2146. (35) Thar, J.; Kirchner, B. J. Phys. Chem. A 2006, 110, 4229–4237. (36) Zahn, S.; Reckien, W.; Kirchner, B.; Staats, H.; Matthey, J.; Ltzen, A. Chem.sEur. J. 2009, 15, 2572–2580. (37) Reiher, M.; Brehm, G.; Schneider, S. J. Phys. Chem. A 2004, 108, 734–742. (38) Kossmann, S.; Thar, J.; Kirchner, B.; Hunt, P. A.; Welton, T. J. Chem. Phys. 2006, 124, 174506. (39) Nockemann, P.; Thijs, B.; Driesen, K.; Janssen, C. R.; Van Hecke, K.; Van Meervelt, L.; Kossmann, S.; Kirchner, B.; Binnemans, K. J. Phys. Chem. B 2007, 111, 5254–5263. (40) Kirchner, B.; Spickermann, C.; Reckien, W.; Schaley, C. A. J. Am. Chem. Soc. 2010, 132, 484–494. (41) Weinhold, F., Ed.; NBO 5.0 Program Manual; University of Wisconsin: 1996. (42) Reed, A. E.; Curtiss, L. A.; Weinhold, F. Chem. ReV. 1988, 88, 899–926.

9536

J. Phys. Chem. A, Vol. 114, No. 35, 2010

(43) Vallejos, M. M.; Angelina, E. L.; Peruchena, N. M. J. Phys. Chem. A 2010, 114, 2855–2863. (44) Steiner, T. Angew. Chem., Int. Ed. 2002, 41, 48–76. (45) Belkova, N. V.; Shubina, E. S.; Epstein, L. M. Acc. Chem. Res. 2005, 38, 624–631. (46) Castro, R. A. E.; Cantilho, J.; Nunes, S. C. C.; Eusebio, M. E. S.; Redinha, J. S. Spectrochim. Acta A 2009, 819–826. (47) Corcelli, S. A.; Skinner, J. L. J. Phys. Chem. A 2005, 109, 6154– 6165. (48) Chen, Y.-F.; Viswanathan, R.; Dannenberg, J. J. J. Phys. Chem. B 2007, 111, 8329–8334. (49) Jesus, A. J. L.; Rosado, M. T. S.; Reva, I.; Fausto, R.; Eusebio, M. E. S.; Redinha, J. S. J. Phys. Chem. A 2008, 112, 4669–4678. (50) Reckien, W.; Kirchner, B.; Peyerimhoff, S. D. J. Phys. Chem. A 2006, 110, 12963–12970. (51) Deshmukh, M. M.; Gadre, S. R. J. Phys. Chem. A 2009, 113, 7927– 7932. (52) Treutler, O.; Ahlrichs, R. J. Chem. Phys. 1995, 102, 346–353. (53) Becke, A. D. Phys. ReV. A 1988, 38, 3098–3100. (54) Lee, C.; Yang, W.; Parr, R. G. Phys. ReV. B 1988, 37, 785–789. (55) Grimme, S. J. Comput. Chem. 2006, 27, 1787–1799. (56) Ahlrichs, R.; Ba¨r, M.; Ha¨ser, M.; Horn, H.; Ko¨lmel, C. Chem. Phys. Lett. 1989, 162, 165–169. (57) Scholze, K. Hydrogen Bond Detection with Quantum Chemical Methods and Frequency Analysis, Master’s thesis, University of Leipzig: 2009. All structures and other data are available upon request. Please contact: [email protected]. ¨ hm, H.; Ha¨ser, M.; Ahlrichs, R. Chem. (58) Eichkorn, K.; Treutler, O.; O Phys. Lett. 1995, 240, 283–289. (59) Eichkorn, K.; Weigend, F.; Treutler, O.; Ahlrichs, R. Theor. Chem. Acc. 1997, 97, 119–124. (60) Weigend, F.; Ha¨ser, M. Theor. Chem. Acc. 1997, 97, 331–340. (61) Weigend, F.; Ha¨ser, M.; Patzelt, H.; Ahlrichs, R. Chem. Phys. Lett. 1998, 294, 143–152. (62) Neugebauer, J.; Reiher, M.; Kind, C.; Hess, B. A. J. Comput. Chem. 2002, 23, 895–910. (63) Weinhold, F.; Landis, C. R. Chem. Educ. Res. Pract. Eur. 2001, 2, 91–104. (64) Frisch, M. J.; Trucks, G. W.; Schlegel, H. B.; Scuseria, G. E.; Robb, M. A.; Cheeseman, J. R.; Montgomery, J. A., Jr.; Vreven, T.; Kudin, K. N.; Burant, J. C.; Millam, J. M.; Iyengar, S. S.; Tomasi, J.; Barone, V.; Mennucci, B.; Cossi, M.; Scalmani, G.; Rega, N.; Petersson, G. A.; Nakatsuji, H.; Hada, M.; Ehara, M.; Toyota, K.; Fukuda, R.; Hasegawa, J.; Ishida, M.; Nakajima, T.; Honda, Y.; Kitao, O.; Nakai, H.; Klene, M.; Li, X.; Knox, J. E.; Hratchian, H. P.; Cross, J. B.; Bakken, V.; Adamo, C.; Jaramillo, J.; Gomperts, R.; Stratmann, R. E.; Yazyev, O.; Austin, A. J.; Cammi, R.; Pomelli, C.; Ochterski, J. W.; Ayala, P. Y.; Morokuma, K.; Voth, G. A.; Salvador, P.; Dannenberg, J. J.; Zakrzewski, V. G.; Dapprich, S.; Daniels, A. D.; Strain, M. C.; Farkas, O.; Malick, D. K.; Rabuck, A. D.; Raghavachari, K.; Foresman, J. B.; Ortiz, J. V.; Cui, Q.; Baboul, A. G.; Clifford, S.; Cioslowski, J.; Stefanov, B. B.; Liu, G.; Liashenko, A.; Piskorz, P.; Komaromi, I.; Martin, R. L.; Fox, D. J.; Keith, T.; Al-Laham, M. A.; Peng, C. Y.; Nanayakkara, A.; Challacombe, M.; Gill, P. M. W.; Johnson,

Wendler et al. B.; Chen, W.; Wong, M. W.; Gonzalez, C.; Pople, J. A. Gaussian 03, reVision D; Gaussian, Inc.: Wallingford, CT, 2004. (65) Werner, H.-J. MOLPRO, 2006. (66) Dunning, J.; Thom, H. J. Chem. Phys. 1989, 90, 1007–1023. (67) Dunning, J.; Thom, H.; Peterson, K. A.; Wilson, A. K. J. Chem. Phys. 2001, 114, 9244–9253. (68) Boys, S. F.; Bernardi, F. Mol. Phys. 1970, 19, 553–566. (69) Boese, A. D.; Martin, J. M. L.; Klopper, W. J. Phys. Chem. A 2007, 111, 11122–11133. (70) Friedrich, J.; Perlt, E.; Roatsch, M.; Spickermann, C.; Kirchner, B. submitted. (71) Arey, J. S.; Aeberhard, P. C.; Lin, I.-C.; Ro¨thlisberger, U. J. Phys. Chem. B 2009, 113, 4726–4732. (72) Bachorz, R. A.; Bischoff, F. A.; Ho¨fener, S.; Klopper, W.; Ottiger, P.; Leist, R.; Frey, J. A.; Leutwyler, S. Phys. Chem. Chem. Phys. 2008, 10, 2758–2766. (73) Morokuma, K. Acc. Chem. Res. 1977, 10, 294–300. (74) Kumar, R.; Schmidt, J. R.; Skinner, J. L. J. Chem. Phys. 2007, 126, 204107. (75) Chandler, D. J. Chem. Phys. 1978, 68, 2959–2970. (76) Li, X.; Liu, L.; Schlegel, H. B. J. Am. Chem. Soc. 2002, 124, 9639– 9647. (77) Kozich, V.; Dreyer, J.; Werncke, W. J. Chem. Phys. 2009, 130, 034505. (78) Dreyer, J. J. Chem. Phys. 2007, 127, 054309. (79) Reiher, M. Found. Chem. 2003, 2, 17–163. (80) Wieczorek, R.; Dannenberg, J. J. J. Am. Chem. Soc. 2003, 125, 8124–8129. (81) Sheu, S.-Y.; Yang, D.-Y.; Selzle, H. L.; Schlag, E. W. Proc. Natl. Acad. Sci. U. S. A. 2003, 100, 12683–12687. (82) Ishimoto, T.; Tokiwa, H.; Teramae, H.; Nagashima, U. J. Chem. Phys. 2005, 122, 094905. (83) Margulis, C. J.; Stern, H. A.; Berne, B. J. J. Phys. Chem. B 2002, 106, 10748–10752. (84) Gilli, G.; Bellucci, F.; Ferretti, V.; Bertolasi, V. J. Am. Chem. Soc. 1989, 111, 1023–1028. (85) Jeffrey, G. A.; Maluszynska, H. Int. J. Biol. Macromol. 1982, 4, 173–185. (86) Derewenda, Z. S.; Lee, L.; Derewenda, U. J. Mol. Bio. 1995, 241, 83–93. (87) Gould, R. O.; Gray, A. M.; Taylor, P.; Walkinshaw, M. D. J. Am. Chem. Soc. 1985, 107, 5921–5927. (88) Vargas, R.; Garza, J.; Dixon, D. A.; Hay, B. P. J. Am. Chem. Soc. 2000, 122, 4750–4755. (89) Vargas, R.; Garza, J.; Friesner, R. A.; Stern, H.; Hay, B. P.; Dixon, D. A. J. Phys. Chem. A 2001, 105, 4963–4968. (90) Deshmukh, M. M.; Bartolotti, L. J.; Gadre, S. R. J. Phys. Chem. A 2008, 112, 312–321. (91) Spickermann, C.; Lehmann, S. B. C.; Kirchner, B. J. Chem. Phys. 2008, 128, 244506.

JP103470E