Hydrogen-bond strengths in large complexes: Efficient calculations using locally dense basis sets

doi:10.1016/j.cplett.2006.12.080

Chemical Physics Letters

Volume 435, Issues 4–6, 19 February 2007, Pages 201-207

https://doi.org/10.1016/j.cplett.2006.12.080 Get rights and content

Abstract

The calculation of hydrogen bond strengths for large systems remains a computationally costly task. Here we show that the appropriate assignment of locally dense basis sets (LDBS) can greatly reduce the cost of such calculations with little or no reduction in accuracy. The use and performance of the LDBS approach with the B971 density functional and the MP2 method are demonstrated on eight small hydrogen bonded systems. Application of the LDBS approach to the adenine–thymine and cytosine–guanine nucleobase complexes shows that results of comparable accuracy to those obtained with balanced basis sets can be achieved with one to two orders of magnitude lower compute times. Similar results are obtained for complexes of tert-butylmethyl ether with five substituted phenols.

Graphical abstract

Using larger basis sets on the 1° region (atoms in boxes) than on the rest of the molecule improves predicted H-bond strengths.

Introduction

Hydrogen bonding (H-bonding) is an extremely important interaction in chemistry and plays a central role in determining the structures of biomolecules (DNA, proteins) and supramolecular constructs [1]. Unlike dissociation enthalpies of covalent bonds, which are well-reproduced by conventional density functional theory techniques and moderate basis sets [2], predicting accurate H-bond strengths for extended systems remains a computational challenge [3]. This is the case because small absolute errors resulting from a correlation- or basis set-limited computational treatment are large relative to the strength of H-bonds themselves.

Halkier et al. performed an excellent comprehensive study of the convergence of hydrogen bond strengths determined with the Hartree–Fock (HF), second-order Møller–Plesset (MP2) and coupled-cluster singles and doubles with perturbative triples (CCSD(T)) approaches as a function of size of aug-cc-pVXZ basis sets (X = D, T, Q, 5) [4]. They showed that H-bond energies converge to the basis set limit from above when counterpoise (CP) corrections [5] for basis set superposition effects are included but converge from below when CP corrections are omitted. With MP2 and CCSD(T), the convergence is slow and unsystematic without CP corrections. However, with these corrections the errors from an incomplete description of the electronic cusp conditions dominate and the convergence becomes monotonic. Tsuzuki et al. found similar results on a different set of hydrogen-bonded complexes [6].

The slow convergence of correlated wavefunction energies with basis size necessitates the use of large basis sets and CP corrections in order to obtain accurate H-bond strengths. Fortunately, the systematic convergence of CP-corrected H-bond strengths can be exploited to reduce the computational costs associated with calculating H-bond strengths close to the basis set limit. One such approach involves the use of basis set extrapolations. For example, Martin’s [7] simple two-point basis set extrapolation of the form $Δ E_{X + 1, X} = \frac{(X + 3 / 2)^{4}}{(X + 3 / 2)^{4} - (X + 1 / 2)^{4}} Δ E_{X + 1} - \frac{(X + 1 / 2)^{4}}{(X + 3 / 2)^{4} - (X + 1 / 2)^{4}} Δ E_{X}$ greatly improves agreement with basis set limit values. This technique extrapolates relative energies (ΔE), such as H-bond energies, to the basis set limit based on values obtained using correlation consistent basis sets and the cardinal numbers (X) of those basis sets. However, for large chemical systems it is often impractical to perform even the aug-cc-pVTZ calculation needed for the lowest level (TZ–DZ or 3-2) extrapolation. One approach to deal with this problem is to use locally dense basis sets (LDBS).

The LDBS approach became popular in NMR calculations after Chestnut’s group demonstrated its utility in computing chemical shifts in large molecules. They showed that such calculations could be performed more quickly and with little loss in accuracy by applying large basis sets only on the NMR chromophores and small basis sets on the other atoms in the molecule [8], [9]. Later, we showed that this approach could be used in a wide range of thermochemical applications including bond dissociation enthalpies, activation energies, proton/electron affinities, solvation energies, and for a limited set of H-bond strengths with little or no loss in accuracy [10]. In our work, we used a large basis set to describe the atoms involved in the region or regions of chemical interest and treated the rest of the system with smaller basis sets or a series of smaller basis sets. We also provided some guidance for applying LDBS to chemical systems of varying complexity. LDBS calculations have much shorter run times than fully balanced large basis set calculations and can alleviate convergence difficulties and hardware limitations.

In this work, we present a more detailed study of the suitability of the LDBS approach for calculating H-bond strengths with the B971 [11] density functional and MP2 methods. In this context, we also examine the effectiveness of the LDBS approach used in conjunction with Martin’s simple, two-point basis set extrapolation scheme [7] with CP corrections. Test calculations are performed on eight small H-bonded complexes. The LDBS approach is then demonstrated by calculating the H-bond strengths of the adenine–thymine (A–T) and cytosine–guanine (C–G) DNA nucleobase pairs and on complexes involving tert-butylmethyl ether with five substituted phenols.

Section snippets

Method of calculation

To test the effectiveness of the LDBS approach in calculating H-bond strengths, we began by considering a set of eight hydrogen-bonded complexes spanning a wide range of interaction energies: ammonia dimer ((NH₃)₂), water dimer ((H₂O)₂), formamide–water (OHCNH₂–OH₂), methanol dimer ((H₃COH)₂), water–ammonia (OH₂–NH₃), hydrogen cyanide–hydrogen fluoride (HCN–HF), formamide dimer ((OHCNH₂)₂), and formic acid dimer ((HCOOH)₂). In order to make comparisons of H-bond strengths as meaningful as

Density functional theory H-bond strengths

We begin by considering the CP-corrected balanced basis set B971 BEs for the set of eight small H-bonded complexes. These data are presented in Table 1. The B971 binding energies converge rapidly to the basis set limit, in a manner similar to Hartree–Fock [4]. The aug-cc-pVTZ BEs are in excellent agreement with the BEs at the basis set limit, with a mean absolute error (MAE) of only 0.03 kcal/mol. Increasing the basis set to aug-cc-pVQZ lowers the MAE to 0.00 kcal/mol.

Also listed in Table 1 are

Conclusions

We showed that the use of locally dense basis sets with CP-corrections and basis set extrapolation is an efficient method to obtain BEs of hydrogen-bonded complexes in good agreement with basis set limit values. To be very effective, LDBS calculations should be used for geometry optimizations of the H-bonded complexes as well. Our tests indicate that there is excellent agreement between the LDBS and balanced basis set geometries provided that appropriate complex partitioning is performed.

Acknowledgements

We are grateful to T. Kubar and P. Hobza for kindly providing us with calculated structures for the A–T and C–G nucleobase complexes. G.A.D. thanks the Centre for Excellence in Integrated Nanotools (University of Alberta) for access to computational resources.

References (23)

J.M.L. Martin
Chem. Phys. Lett.
(1996)
E.R. Johnson et al.
Chem. Phys. Lett.
(2006)
M. Feyereisen et al.
Chem. Phys. Lett.
(1993)
H. Fenniri et al.
J. Am. Chem. Soc.
(2002)
E.R. Johnson et al.
J. Phys. Chem. A
(2003)
P. Jurečka et al.
Phys. Chem. Chem. Phys.
(2006)
A. Halkier et al.
J. Chem. Phys.
(1999)
S.B. Boys et al.
Mol. Phys.
(1970)
S. Tsuzuki et al.
J. Chem. Phys.
(1999)
D.B. Chestnut et al.
J. Comput. Chem.
(1989)

D.B. Chestnut et al.

J. Comput. Chem.

(1993)

Cited by (12)

Structural and energetic quantum chemical investigations into how the bioactive thiazolidinedione and rhodanine scaffolds interact with cytosine to form part of DNA
2018, Computational and Theoretical Chemistry
Citation Excerpt :
The strength of the hydrogen bond lies between 1 and 2 kJ/mol (very weak) and 161.5 kJ/mol (very strong) [5,6]. Hydrogen bonding interactions are also mainly electrostatic with some contribution of charge transfer and polarization in nature [7–13]. Many bioactive molecules could act as both hydrogen bond acceptor and donor, and thus investigations including the effects of hydrogen bond interactions on the properties, activity and nature of the interactions of biologically active molecules are of great interest and much computational and experimental research works has been focused on them [14–35].
The interaction within biologically active complex configurations composed of thiazolidinedione-cytosine (TC) and rhodanine-cytosine (RC) have been fully investigated using B3LYP and M062X methods in conjunction with various basis sets, including 6-311++G(d,p), 6-311++G(2d,2p), 6-311++G(df,pd) and AUG-cc-pVDZ. Dispersion corrections to interaction energies using M06-2X-GD3 and double hybrid density functionals (B2PLYP-GD2, B2PLYP-GD3 and mPW2PLYP-GD2) have also been taken into account. The results showed that the interaction energies, including long range dispersion corrections obtained using double hybrid density functionals, are about 1–2.5 kcal/mol higher than uncorrected ones. Various combinations of all possible sites around the monomers which could lead to formation of stable complex configurations were considered, and then six double hydrogen bonded interactions which are formed in a cyclic pattern were further investigated for both complex configurations between cytosine-thiazolidinedione and cytosine-rhodanine. The intermolecular hydrogen bonds which are involved in forming of cyclic interaction patterns are of NH⋯O(N) and CH⋯O(N) types.
The TCS′2S1 and RCS′2S1 complex configurations were found to have the largest interaction energies within all the complex configurations studied for cytosine-thiazolidinedione and cytosine-rhodanine respectively at all computational levels used for the calculation of interaction energies. These two most stable complex configurations formed between the S1 site of the cytosine monomer and the S′2 site of the thiazolidinedione or rhodanine.
The BSSE-corrected interaction energy for six cyclic double hydrogen bonded complex configurations, calculated at M062X/6-311++G(df,pd) level range from −11.26 to −16.22 kcal/mol for TC complexes and −11.04 to −15.59 kcal/mol for RC complexes.
The stability order within thiazolidinedione-cytosine (TC) complexes is TCS′2S1 > TCS′1S1 > TCS′2S2 > TCS′1S2 > TCS′3S2 > TCS′3S1 in all levels studied. The same trend was observed for the stability order of rhodanine-cytosine (RC) complexes.
The NBO results reveal that the charge transfer occurred from cytosine to thiazolidinedione or rhodanine within complex configurations which are formed at S′1S1, S′1S2, S′2S1 and S′2S2 sites, whereas in the case of complex configurations of cytosine with thiazolidinedione or rhodanine at S′3S1 and S′3S2 sites this order was reversed.
The atom in molecule theory (AIM) and natural bond orbital analysis (NBO) were used to understand the nature of the interactions within the complex configurations studied. This showed good consistency with the interaction energies. In addition, good correlations between geometrical parameters, interaction energies, topological AIM data and NBO results were obtained. Based on AIM and NBO results, it is concluded that the H-bond interactions are electrostatic and partially covalent. Finally, the relationship between interaction energies and vibrational frequency changes at involved H-bonds was analyzed.
A quantum chemical insight to intermolecular hydrogen bonding interaction between cytosine and nitrosamine: Structural and energetic investigations
2016, Journal of Molecular Structure
Citation Excerpt :
The Hydrogen bonding type in many of DNA molecules is as N–H⋯O and N–H⋯N. Because of theirs importance, they have been widely investigated. The results of investigations indicated that the nature of hydrogen bonding interactions is mainly electrostatic with some contribution of charge transfer and polarization [4–10]. Hydrogen bonding also stabilizes the structures of biological molecules such as proteins and nucleic acids.
Hydrogen bond interactions which are formed during complex formation between cytosine and nitrosamine have been fully investigated using B3LYP, B3PW91 and MP2 methods in conjunction with various basis sets including 6–311++G (d,p), 6–311++G (2d,2p), 6–311++G (df,pd) and AUG-cc-pVDZ. Three regions around the most stable conformer of cytosine in the gas phase with six possible double H-bonded interactions were considered.
Two intermolecular hydrogen bonds of type N_C–N-H_NA and O–H(N–H)_C–O_NA were found on the potential energy surface in a cyclic system with 8-member in CN1, CN3, CN5 and 7-member in CN2, CN4, CN6 systems.
Results of binding energy calculation at all applied methods reveal that the CN1 structure is the most stable one which is formed by interaction of nitrosamine with cytosine in S1 region. The BSSE-corrected binding energy for six complex system is ranging from −23.8 to −43.6 kJ/mol at MP2/6–311++G (df,pd) level and the stability order is as CN1 > CN2 > CN3 > CN4 > CN5 > CN6 in all studied levels of theories. The NBO results reveal that the charge transfer occurred from cytosine to nitrosamine in CN1, CN3, CN5 and CN6 whereas this matter in the case of CN2 and CN4 was reversed.
The relationship between BEs with red shift of H-bond involved bonds vibrational frequencies, charge transfer energies during complex formation and electron densities at H-bond BCPs were discussed. In addition activation energetic properties related to the proton transfer process between cytosine and nitrosamine have been calculated at MP2/6–311++G (df,pd) level.
AIM results imply that H-bond interactions are electrostatic with partially covalent characteristic in nature.
A systematic study on hydrogen bonding interaction between formamide and cytosine
2010, Journal of Molecular Structure: THEOCHEM
A systematic computational study was carried out to characterize the hydrogen bonding of complexes formed between formamide and cytosine by DFT calculations. The computations were performed mainly with the B3LYP/6-311++G(d,p) level. Seven stable cyclic structures are found on the potential energy surface, in which four structures have two normal hydrogen bonds and the others have only one normal hydrogen bond with a very weak hydrogen bond that can be neglected. In the four structures with two normal hydrogen bonds, two have seven-membered rings, and the others have an eight-membered ring. The eight-membered ring is preferred to the seven-membered one by analyzing the hydrogen bond lengths and the interaction energies. The infrared spectrum frequencies, vibrational frequency shifts and charge number are also reported.
DFT Investigation of a Charge-Transfer Complex Formation Between p-Phenylenediamine and 3,5-Dinitrosalicylic Acid
2019, Journal of Structural Chemistry
Understanding the Many-Body Basis Set Superposition Error: Beyond Boys and Bernardi
2018, Journal of Chemical Theory and Computation
Performance of conventional and dispersion-corrected density-functional theory methods for hydrogen bonding interaction energies
2013, Physical Chemistry Chemical Physics

View all citing articles on Scopus

¹: Present address: Department of Chemistry, Dalhousie University, Halifax, Canada NS B3H 4J3.

View full text

Hydrogen-bond strengths in large complexes: Efficient calculations using locally dense basis sets

Abstract

Graphical abstract

Introduction

Section snippets

Method of calculation

Density functional theory H-bond strengths

Conclusions

Acknowledgements

Chem. Phys. Lett.

Chem. Phys. Lett.

Chem. Phys. Lett.

J. Am. Chem. Soc.

J. Phys. Chem. A

Phys. Chem. Chem. Phys.

J. Chem. Phys.

Mol. Phys.

J. Chem. Phys.

J. Comput. Chem.

J. Comput. Chem.