Predicting C-H/ interactions with nonlocal density functional theory
We examine the performance of a recently developed nonlocal density functional in predicting a model noncovalent interaction, the weak bond between an aromatic system and an aliphatic C-H group. The new functional is a significant improvement over traditional density functionals, providing results which compare favorably to high-level quantum-chemistry techniques but at considerably lower computational cost. Interaction energies in several model C-H/ systems are in generally good agreement with coupled-cluster calculations, though equilibrium distances are consistently overpredicted when using the revPBE functional for exchange. The new functional correctly predicts changes in energy upon addition of halogen substituents.
Accurate treatment of van der Waals bonding represents one of the primary challenges for modern quantum chemical calculations. Such interactions are important in a great diversity of chemical and biological systems, and have been difficult to model in a highly accurate and computationally tractable manner. State-of-the-art techniques such as coupled cluster theory (CCSD(T)) generally treat such interactions very accurately, but are computationally demanding and often depend sensitively on the basis set used. Density functional theory (DFT) methods are attractive in terms of computational efficiency, but common exchange-correlation functionals used do not include effects such as dispersion interactions. For this reason, traditional DFT calculations frequently provide a poor description of van der Waals bonding in molecular systems.
Here we consider a recently developed van der Waals density functional (vdW-DF) which explicitly includes nonlocal electron correlation without the use of any empirical parameters Dion et al. (2004); Thonhauser et al. (2006a). This method has been successfully applied to a number of aromatic systems such as dimers of benzene Puzder et al. (2006); Thonhauser et al. (2006b), napthalene, anthracene, and pyrene Chakarova and Schroder (2005), as well as a variety of bulk systems such as crystalline polyethylene Kleis et al. (2006). The vdW-DF functional was strongly influenced by the physics included in its forerunner Rydberg et al. (2000, 2003); Langreth et al. (2005).
Of key importance is the ability of this nonlocal functional to treat a range of different noncovalent interactions. Along these lines we consider the attractive bonding between an aliphatic C-H group and an aromatic system, a configuration which has received considerable attention in recent years Nishio et al. (1998); Desiraju and Steiner (1999). This so-called C-H/ interaction is frequently referred to as a weak hydrogen bond, though recent studies have demonstrated that it arises primarily from dispersion, with electrostatic attraction playing a minor role Tsuzuki et al. (1999); Tsuzuki et al. (2000a); Ringer et al. (2006). The C-H/ interaction is important in molecular packing Takahashi et al. (2000), protein interactions Nishio et al. (1995); Muraki (2002), and in a variety of biological systems Brandl et al. (2001); Watanabe et al. (2000). The wide interest in this interaction, combined with the availability of high level CCSD(T) and MP2 calculations, make it a valuable test case for the new van der Waals functional.
In this paper we consider a number of model systems which allow careful study of the C-H/ interaction. These include the methane-benzene system, trifluoro- and trichloromethane with benzene, and the methane-indole complex. All results are compared with high-level CCSD(T) and MP2 calculations. The interaction energies from vdW-DF are a significant improvement over gradient-corrected and hybrid DFT functionals, and compare favorably to CCSD(T) results extrapolated to the limit of a complete basis set. Equilibrium distances are generally overpredicted by vdW-DF . Our results suggest that this functional, which requires little time beyond a standard DFT calculation, holds great promise for modeling large organic and biological systems.
Ii Nonlocal density functional
Our approach to calculating the C-H/ interaction uses a recent density functional which explicitly includes van der Waals correlations Dion et al. (2004); Thonhauser et al. (2006a); Rydberg et al. (2003); Langreth et al. (2005); Rydberg et al. (2000). In this framework, which we will abbreviate as vdW-DF, the total energy functional is written as
where is the independent-particle kinetic energy, is the ionic pseudopotential, and is the Hartree energy. For the exchange-correlation functional, we separate the exchange energy , the local correlation , and the nonlocal component . Following previous work, we use the revPBE expression for exchange Zhang and Yang (1998) for all calculations. is approximated here by the LDA correlation. The nonlocal term, which is key for capturing long-range correlations such as dispersion, can be expressed in a general form as
The kernel depends only on the distance and the density near and . Details regarding explicit evaluation of the kernel are given in Ref. Dion et al. (2004). No empirical parameters are used in the nonlocal correlation. Though the revPBE exchange does not satisfy the Lieb-Oxford bound for any arbitrary density as does PBE, in practice it satifies the true, integrated Lieb-Oxford condition for typical systems Zhang and Yang (1998); Hammer et al. (1999). Thus our overall approach, using vdW-DF and revPBE exchange, may be considered a truly first-principles method for treating van der Waals systems.
In previous work was computed as a post-processing step using the charge density from a normal DFT calculation Puzder et al. (2006); Thonhauser et al. (2006b). Explicit evaluation of the nonlocal potential in a self-consistent way is possible through recently derived analytical expressions Thonhauser et al. (2006a). However, test calculations for the methane-benzene system indicated little difference in the interaction energies between these two approaches. Thus for all systems we continue to utilize the post-processing method, which requires very little time beyond a standard DFT calculation.
All vdW-DF calculations were done in a modified version of abinit Gonze et al. (2002). Norm-conserving Troullier-Martins pseudopotentials were used, with a kinetic energy cutoff of 50 Rydberg. Calculations were performed in a large cubic cell with at least 18 Å of space between adjacent periodic images. The nonlocal energy was calculated on a real-space grid identical to the one used for Fourier transforms in the planewave code. Calculations for hybrid DFT functionals were performed using CRYSTAL03 Saunders et al. (2003).
Molecular geometries were optimized using the B3LYP hybrid functional with an aug-cc-pVTZ basis set; this optimization gave structures that were generally within 0.01 Å of previous high-level calculations for methane-benzene and methane-indole Ringer et al. (2006). These same geometries were also converged with respect to forces derived from the self-consistent implementation of vdW-DF, indicating that the nonlocal functional vanishes appropriately in intramolecular situations. In all cases the isolated molecule geometry was kept rigid during the two-molecule calculations.
Iii Results and discussion
iii.1 Methane-benzene system
The methane-benzene complex represents one of the simplest interacting C-H/ systems. Following the results of previous experimental Schauer and Bernstein (1985); Shibasaki et al. (2006) and theoretical Tsuzuki et al. (2000b); Ringer et al. (2006) work, the molecules are oriented such that a single C-H group on the methane points directly at the center of the benzene, as shown in Fig. 1. The separation is then defined as the distance between the central C atom in methane and the center of the benzene ring. The separation axis is thus coincident with the symmetry axis of benzene, and the dimer system overall is set in a symmetry. Similar to the report by Ringer and coauthors Ringer et al. (2006), the dimer energy is essentially invariant for rotations around the dimer separation axis, varying by at most kcal/mol.
In Fig. 2 we present our results for the interaction energy of the methane-benzene dimer. The vdW-DF functional is a significant improvement over other standard DFT approaches, yielding an interaction energy comparable to the CCSD(T) results of Ringer et al Ringer et al. (2006). Their coupled-cluster interaction energy was extrapolated to the limit of a complete basis set (CBS) using Helgaker’s method with large augmented correlation-consistent triple- and quadruple- basis sets Ringer et al. (2006). Shibasaki and coauthors also report similar theoretical CCSD(T) results Shibasaki et al. (2006). Despite its relative simplicity, vdW-DF captures the C-H/ energy quite accurately as compared to these high level calculations. The equilibrium dimer separation predicted by vdW-DF is roughly 8% larger than the CCSD(T) result, which appears to be a general trend of vdW-DF when using the revPBE functional for the exchange energy Puzder et al. (2006); Thonhauser et al. (2006b, a); Kleis et al. (2006).
Three results using traditional DFT methods are also shown in Fig. 2. The Perdew-Burke-Ernzerhof (PBE) nonempirical generalized gradient (GGA) functional produces a binding that is considerably weaker than that given by vdW-DF, though its predicted equilibrium distance is only slightly higher than the nonlocal functional (see Table 1 for numerical values). Zhang and Yang’s revised PBE functional (revPBE), which is identical to PBE except for a semiempirical adjustment of one parameter in the exchange portion Zhang and Yang (1998), produces no binding. The popular hybrid functional B3LYP, which uses a semiempirical mixing of exact exchange and GGA terms Becke (1993), also produces no binding for the system. Recent calculations by Shibasaki and coauthors show similar results for a range of GGA and hybrid functionals Shibasaki et al. (2006). Thus, as we would expect for a dispersion-dominated dimer, traditional DFT methods compare poorly to high level CCSD(T) and MP2 results.
Shibasaki and coworkers have recently measured the experimental interaction energy of gas phase methane-benzene clusters by mass analyzed threshold ionization spectroscopy Shibasaki et al. (2006). For comparison with the experimental value we must subtract the vibrational zero-point energy (ZPE) from our predicted C-H/ interaction energy. The predicted change in ZPE for formation of the methane-benzene dimer is 0.296 kcal/mol using MP2 with a cc-pVTZ basis set and scaling factors to match the experimental and theoretical frequencies for the methane symmetric stretching mode Shibasaki et al. (2006). Subtracting this gives an interaction energy of 1.27 kcal/mol using vdW-DF and 1.15 kcal/mol using CCSD(T); both of these compare favorably to the experimental binding energy, which was determined to be between 1.03 and 1.13 kcal/mol Ringer et al. (2006); Shibasaki et al. (2006).
iii.2 Trifluoro- and trichloromethane
We next consider the effect of halogen substituents on the aliphatic molecule. Some experimental measurements of the role of such substituent atoms on the C-H/ interaction in solution have been reported, including evidence that benzene forms a stable dimer with trichloromethane Reeves and Schneider (1957); Zurcher (1961). Here we consider replacement of three hydrogen atoms on the methane by chlorine or flourine, both of which increase the interaction energy of the complex. The results for these dimers of trifluoro- and trichloromethane with benzene are given in Fig. 3 and Table 1.
In both systems the C-H bond points directly at the center of the benzene ring, and the remaining three H atoms on the methane are replaced. The geometries for CHF and CHCl were optimized independently, and the monomer structure was also used for dimer calculations. Trifluoro- and trichloromethane both have increased molecular polarizabilities compared to CH. Methane has a polarizability of roughly 17.9 a.u., compared to 19.0 a.u. for trifluoromethane and 57.6 a.u. for trichloromethane van Duijnen and Swart (1998). Since dispersion is the dominant component of the C-H/ interaction, we would expect an increased interaction energy for the molecules with larger polarizabilities. This is indeed observed, as seen in Fig. 3; the benzene-trichloromethane binding is much stronger than that seen in the methane system, and is comparable to the strength of the hydrogen bond in a dimer of water. We would also expect an increase in electrostatic attraction upon substitution of F or Cl, though previous results indicate this is a much smaller effect than the enhanced dispersion Tsuzuki et al. (2002).
The performance of vdW-DF in these halogenated dimers is similar to that seen in the methane-benzene system. The nonlocal functional is a considerable improvement over the prototypical GGA functional PBE, providing interaction energies comparable to CCSD(T) results extrapolated to the limit of a complete basis set by Tsuzuki and coauthors Tsuzuki et al. (2002). In these systems the vdW-DF interaction energy underestimates the previously reported CCSD(T) limit, in contrast to a slight overestimation for the methane-benzene and methane-indole complexes. The literature values for the latter two systems used slightly different basis sets as well as different methods for extracting the complete basis limit than were used for the halogenated dimers Tsuzuki et al. (2002); Ringer et al. (2006). We note again the advantage of using a simple, straightforward DFT method such as vdW-DF.
The vdW-DF functional with revPBE exchange again overstimates the equilibrium separation of the dimers by approximately 6% for trifluoromethane-benzene and 13% for trichloromethane-benzene. MP2 and CCSD(T) calculations indicate that the latter has a smaller dimer separation, as would be expected from the increased dispersion component. vdW-DF does not reproduce this trend, predicting similar separations for both halogenated dimers.
We next consider the interaction of methane with a larger aromatic system, the indole molecule. The methane is placed over the five-membered ring of indole, again with a single C-H bond pointing directly at the center of the ring (see Fig. 1). Geometries of the monomers were optimized separately and kept rigid during dimer calculations. In this configuration the interaction energy is again essentially invariant to rotations around the axis connecting the methane C and the center of the five-membered ring.
Our results for the interaction energy of the methane-indole dimer are given in Fig. 4 and Table 1. MP2 and CCSD(T) results are taken from the work of Ringer and coauthors, where an aug-cc-pVTZ basis set was used Ringer et al. (2006). The magnitude of the interaction and the dimer separation are similar to the methane-benzene dimer, though vdW-DF predicts that methane-indole is more strongly bound, in agreement with CCSD(T) and MP2 results. An analysis by symmetry-adapted perturbation theory shows a similar distribution of components in the methane-benzene and methane-indole dimers, the most significant difference being a larger dispersion component (0.589 kcal/mol more negative) in the latter. As before, vdW-DF overpredicts the CCSD(T) dimer separation by approximately 5%.
We have studied the performance of the van der Waals density functional in treating the weak bond between an aromatic system and an aliphatic C-H group. The properties of four prototypical C-H/ systems were analyzed: methane-benzene, trifluoro- and trichloromethane with benzene, and methane-indole.
The vdW-DF method is a substantial improvement over standard DFT functionals, providing in all cases an interaction energy that is comparable to high-level quantum chemistry calculations. In the case of the methane-benzene dimer, vdW-DF predicts a total interaction energy of 1.57 kcal/mol, compared to 1.454 kcal/mol from CCSD(T) and 1.79 kcal/mol from MP2 Ringer et al. (2006). Traditional DFT functionals predict little or no binding for this dimer, and generally wavefunction-based methods require large basis sets to achieve values in this range Shibasaki et al. (2006); Tsuzuki et al. (2000b).
The substitution of F and Cl for three of the methane hydrogens results in an increased binding due to larger dispersion and electrostatic interactions. The vdW-DF functional correctly reproduces the increased interaction energy in good agreement with the CCSD(T) trend, but does not match the reduction in equilibrium dimer separation. vdW-DF is also able to resolve the small binding energy difference between the methane-benzene and methane-indole dimers, but overpredicts the equilibrium separation in the latter system as well.
Overall, the performance of the new van der Waals density functional for these very weak bonds is encouraging. With further development, this computationally tractable method holds great promise for complex organic and biological systems where this manner of weak bonding plays a key role.
The authors would like to thank B. Rice, E. Byrd, and W. Mattson for useful discussions. J. H. acknowledges an NRL/ASEE Postdoctoral Fellowship at the Naval Surface Warfare Center, Indian Head. N. A. R. was supported by an NRC Research Associateship Award at the U.S. Army Research Lab. Work at the DoD was supported by the Office of Naval Research. Work at Rutgers was supported in part by NSF Grant No. DMR-0456937.
- Dion et al. (2004) M. Dion, H. Rydberg, E. Schroder, D. C. Langreth, and B. I. Lundqvist, Phys. Rev. Lett. 92, 246401 (2004).
- Thonhauser et al. (2006a) T. Thonhauser, V. R. Cooper, S. Li, A. Puzder, P. Hyldgaard, and D. C. Langreth (2006a), eprint cond-mat/0703442.
- Puzder et al. (2006) A. Puzder, M. Dion, and D. C. Langreth, J. Chem. Phys. 124, 164105 (2006).
- Thonhauser et al. (2006b) T. Thonhauser, A. Puzder, and D. C. Langreth, J. Chem. Phys. 124, 164106 (2006b).
- Chakarova and Schroder (2005) S. Chakarova and E. Schroder, J. Chem. Phys. 122, 054102 (2005).
- Kleis et al. (2006) J. Kleis, B. Lundqvist, D. Langreth, and E. Schroder (2006), eprint cond-mat/0611498.
- Rydberg et al. (2000) H. Rydberg, B. I. Lundqvist, D. C. Langreth, and M. Dion, Phys. Rev. B 67, 6997 (2000).
- Rydberg et al. (2003) H. Rydberg, M. Dion, N. Jacobson, E. Schroder, P. Hyldgaard, S. Simak, D. C. Langreth, and B. I. Lundqvist, Phys. Rev. Lett. 91, 126402 (2003).
- Langreth et al. (2005) D. Langreth, M. Dion, H. Rydberg, E. Schroder, and B. Lundqvist, Int. J. Quantum Chem. 101, 559 (2005).
- Nishio et al. (1998) M. Nishio, M. Hirota, and Y. Umezawa, The CH/ Interaction (Wiley-VCH, New York, 1998).
- Desiraju and Steiner (1999) G. R. Desiraju and T. Steiner, The Weak Hydrogen Bond (Oxford University Press, New York, 1999).
- Tsuzuki et al. (1999) S. Tsuzuki, K. Honda, T. Uchimaru, M. Mikami, and K. Tanabe, J. Phys. Chem. A 103, 8265 (1999).
- Tsuzuki et al. (2000a) S. Tsuzuki, K. Honda, T. Uchimaru, M. Mikami, and K. Tanabe, J. Am. Chem. Soc. 122, 3746 (2000a).
- Ringer et al. (2006) A. Ringer, M. Figgs, M. Sinnokrot, and C. D. Sherrill, J. Phys. Chem. A 110, 10822 (2006).
- Takahashi et al. (2000) H. Takahashi, S. Tsuboyama, Y. Umezawa, K. Honda, and M. Nishio, Tetrahedron 56, 6185 (2000).
- Nishio et al. (1995) M. Nishio, Y. Umezawa, M. Hirota, and Y. Takeuchi, Tetrahedron 51, 8665 (1995).
- Muraki (2002) M. Muraki, Protein Peptide Lett. 9, 195 (2002).
- Brandl et al. (2001) M. Brandl, M. Weiss, A. Jabs, J. Suhnel, and R. Hilgenfeld, J. Mol. Biol. 307, 357 (2001).
- Watanabe et al. (2000) T. Watanabe, T. Suzuki, Y. Umezawa, T. Takeuchi, M. Otsuka, and K. Umezawa, Tetrahedron 56, 741 (2000).
- Zhang and Yang (1998) Y. Zhang and W. Yang, Phys. Rev. Lett. 80, 890 (1998).
- Hammer et al. (1999) B. Hammer, L. Hansen, and J. K. Nørskov, Phys. Rev. B 59, 7413 (1999).
- Gonze et al. (2002) X. Gonze, J. M. Beuken, R. Caracas, F. Detraux, M. Fuchs, G. M. Rignanese, L. Sindic, M. Verstraete, G. Zerah, F. Jollet, et al., Comp. Mat. Sci 25, 478 (2002).
- Saunders et al. (2003) V. R. Saunders, R. Dovesi, C. Roetti, R. Orlando, C. M. Zicovich-Wilson, N. M. Harrison, K. Doll, B. Civalleri, I. J. Bush, P. D’Arco, et al., CRYSTAL2003 User’s Manual (2003).
- Schauer and Bernstein (1985) M. Schauer and E. R. Bernstein, J. Chem. Phys. 82, 726 (1985).
- Shibasaki et al. (2006) K. Shibasaki, A. Fujii, N. Mikami, and S. Tsuzuki, J. Phys. Chem. A 110, 4397 (2006).
- Tsuzuki et al. (2000b) S. Tsuzuki, K. Honda, T. Uchimaru, M. Mikami, and K. Tanabe, J. Am. Chem. Soc. 122, 3746 (2000b).
- Perdew et al. (1996) J. P. Perdew, K. Burke, and M. Ernzerhof, Phys. Rev. Lett. 77, 3865 (1996).
- Becke (1993) A. Becke, J. Chem. Phys. 98, 5648 (1993).
- Tsuzuki et al. (2002) S. Tsuzuki, K. Honda, T. Uchimaru, M. Mikami, and K. Tanabe, J. Phys. Chem. A 106, 4423 (2002).
- Reeves and Schneider (1957) L. W. Reeves and W. G. Schneider, Can. J. Chem. 35, 251 (1957).
- Zurcher (1961) R. F. Zurcher, Helv. Chim. Acta 44, 1380 (1961).
- van Duijnen and Swart (1998) P. van Duijnen and M. Swart, J. Phys. Chem. A 102, 2399 (1998).