New Research In
Physical Sciences
Social Sciences
Featured Portals
Articles by Topic
Biological Sciences
Featured Portals
Articles by Topic
 Agricultural Sciences
 Anthropology
 Applied Biological Sciences
 Biochemistry
 Biophysics and Computational Biology
 Cell Biology
 Developmental Biology
 Ecology
 Environmental Sciences
 Evolution
 Genetics
 Immunology and Inflammation
 Medical Sciences
 Microbiology
 Neuroscience
 Pharmacology
 Physiology
 Plant Biology
 Population Biology
 Psychological and Cognitive Sciences
 Sustainability Science
 Systems Biology
The X3LYP extended density functional for accurate descriptions of nonbond interactions, spin states, and thermochemical properties

Contributed by William A. Goddard III, December 30, 2003
Abstract
We derive the form for an exact exchange energy density for a density decaying with Gaussianlike behavior at long range. Based on this, we develop the X3LYP (extended hybrid functional combined with Lee–Yang–Parr correlation functional) extended functional for density functional theory to significantly improve the accuracy for hydrogenbonded and van der Waals complexes while also improving the accuracy in heats of formation, ionization potentials, electron affinities, and total atomic energies [over the most popular and accurate method, B3LYP (Becke threeparameter hybrid functional combined with Lee–Yang–Parr correlation functional)]. X3LYP also leads to a good description of dipole moments, polarizabilities, and accurate excitation energies from s to d orbitals for transition metal atoms and ions. We suggest that X3LYP will be useful for predicting ligand binding in proteins and DNA.
The development of accurate functionals has made density functional theory (DFT) the method of choice for first principles predictions of fundamental processes in materials ranging from metal alloys to semiconductors, to ceramics, to new catalysts (1, 2). Despite this progress, there remain serious limitations in DFT theory. Thus, the B3LYP (Becke threeparameter hybrid functional combined with Lee–Yang–Parr correlation functional) method (3) achieves a high accuracy [mean absolute deviation (MAD) = 0.13 eV (1 eV = 1.602 × 10^{19} J)] for thermochemistry [heats of formation of the 148 molecules in the extended G2 reference set (4, 5)], but it predicts that the noble gas dimers are unstable. The PW1PW hybrid method (6, 7) leads to less accuracy (0.23 eV for G2) and far too strong bonding in noble gas dimers [≈7 times the correct answer for He_{2} (7)], and indeed leads to very strong bonding even when the functional for electron correlation (which is responsible for London dispersion forces) is omitted. Similar results are obtained for other functionals of this class [mPW (7) and PBE (8–10)], suggesting that these exchange functionals include some correlation effects, making it difficult to combine them with true correlation functionals.
The particular application motivating us to reexamine the functionals in DFT is the possibility of genomewide structurebased drug design. As the genomics revolution leads first to sequences and then to 3D structures for all proteins of all forms of life, there will be an opportunity for computation and theory to help develop new generations of drugs (agonists and antagonists) that are both very active and very specific (binding maybe to just one protein of all proteins for all forms of life, so as to minimize toxic side effects). However, for theory and computation to play this role, it is essential that the noncovalent interactions of ligands to proteins be accurately predicted. Thus, it is essential to accurately describe London dispersion forces (van der Waals attraction) along with electrostatic and hydrogen bond interactions.
The current generations of DFT methods do not provide this accuracy, but we present here an extended functional, denoted as X3LYP (extended hybrid functional combined with Lee–Yang–Parr correlation functional), that significantly improves the accuracy for van der Waals and hydrogenbonded complexes while providing an excellent description of dipole moments and polarizabilities. It also improves the accuracy in thermochemistry, ionization potentials (IPs), electron affinities (EAs), and total atomic energies over the most popular and accurate previous method, B3LYP. We expect that X3LYP will prove useful for predicting ligand binding in proteins and DNA.
Theory
Because the magnitude of the correlation energy is generally <10% of the exchange energy, it is most important that the exchange functional be accurate. In GGA (generalized gradient approximation), the exchange energy density is defined as Here, A_{x} = 3/4 (3/π)^{1/3}, F(s) is the GGA enhancement factor, and s is a dimensionless gradient (6), defined as The B88 (11) and PW91 (6) exchange functionals use the forms and where a _{1} = (48π^{2})^{1/3}, a _{2} = 6β·a _{1}, , a _{4} = 10/81  a _{3}, , and d = 4. From fitting the HF exchange energies of the noble gas atoms, Becke obtained β = 0.0042 (11).
The plot of the F ^{B88}(s) and F ^{PW91} (s) functions in Fig. 1 shows that these two functions differ significantly for large s, the region believed to affect significantly the performance of DFT in van der Waals systems (12, 13). For a system with spherically symmetric decaying density, i.e., F ^{B88}(s) assures the correct asymptotic behavior of the exchange energy density (11), i.e., (Condition 1); while Levy and Perdew showed that some scaling properties can be satisfied if the asymptotic form of the functional for large s is s ^{α}, where α ≥ 1/2 (Condition 2) (14). The Lieb–Oxford bound (12, 15), ε _{x} ≥ 1.679 · ρ(r)^{1/3}, suggests an upper bound on F(s) (9) (Condition 3).
B88 violates Conditions 2 and 3, whereas the a _{5}·s ^{d} term lets PW91 obey these conditions. However, PW91 violates Condition 1.
With modern DFT codes, most calculations on finite molecules use Gaussian basis functions, leading to a Gaussianlike long range behavior in the electron density as in Eq. 5 . The asymptotic limit for the exchange energy density of a finite system is where U _{c} is the Coulomb potential of the exchange charge, defined as Inserting Eq. 5 into Eq. 7 leads to . Since Condition 1 is fulfilled for a Gaussianlike density.
Combining Eqs. 1 and 6 with Eqs. 5 and 7 gives Using Eq. 5 , we can rewrite Eq. 2 as
Eqs. 8 and 9 determine the F ^{Gauss}(s) for a Gaussianlike asymptotic density, which Fig. 1 shows to lies between F ^{B88}(s) and F ^{PW91}(s) for s ≥ 1.5, but closer to F ^{B88}(s).
As s → 0, F ^{Gauss}(s) → (2^{5} π/3^{4})^{1/3} = 1.07466, instead of 1.0 as required to obey the limit within the local density approximation (LDA). Such a deviation may be expected for finite systems. Indeed, fitting to the unrestricted Hartree–Fock energies of the first and the secondrow atoms, Handy et al. (16) recently developed a local exchange functional, OPTX, that has an s → 0 LDA limit of 1.05151.
X3LYP Functional
Based on the F ^{Gauss}(s) behavior for s ≥ 1.5 as shown in Fig. 1, we propose an extended exchange functional: Here, we choose to obey the LDA limit as s → 0, as usual in the GGA framework. We could have allowed a more general form for F ^{X}(s) but found that combining the B88 and PW91 functionals provided sufficient flexibility. We considered that using well known functions already incorporated into computer codes could help ease the application of the X3LYP functional.
Combining F ^{X}(s) with the LYP correlation functional (17), leads to XLYP, where the mixing parameters {a_{x} _{1,} a_{x} _{2}} = {0.722, 0.347} were determined through a leastsquare fitting to the total energies of 10 atoms {H, He, Li, Be, B, C, N, O, F, Ne}; the IPs for 16 atoms {Li, Be, B, C, N, O, F, Ne, Na, Mg, Al, Si, P, S, Cl, Ar}; the EAs for 10 atoms {H, B, C, O, F, Al, Si, P, S, Cl}; and the atomization energies for 38 molecules {H_{2}, He_{2}, Li_{2}, Be_{2}, C_{2}, N_{2}, O_{2}, F_{2}, Ne_{2}, Na_{2}, Mg_{2}, Si_{2}, P_{2}, S_{2}, Cl_{2}, CN, CO, CS, NO, SO, ClO, SiO, ClF, PF, AlF, SiF, CCl, SiCl, NaCl, CH, NH, OH, HF, CO_{2}, O_{3}, SO_{3}, OCS, CS_{2}}, selected to represent the important chemistry for the first and the secondrow elements (including open and closedshell molecules; molecules with single, double, and triple bonds; ionic systems, and systems requiring multiple determinants for proper descriptions). In particular, we include He_{2} and Ne_{2} as representative van der Waals systems.
All calculations used a pruned (75,434) DFT grid from jaguar (Version 4.0, Schrödinger, Portland, OR). To facilitate comparisons with literature values, we did not use the pseudospectral capabilities in jaguar. In determining the two XLYP parameters, we used the augccpVTZ basis sets (18) for all calculations. To validate the accuracy against the G2 data set, all calculations used the 6–311+G (3df, 2p) basis set with previously reported MP2 molecular geometries and scaled Hartree–Fock vibrational frequencies to calculate zeropoint energies and finitetemperature corrections (4, 5, 19). This choice of geometries and basis sets allows a direct comparison of our results with previously published data obtained with other functionals (5, 7, 8, 10). For He_{2}, Ne_{2}, and (H_{2}O)_{2}, we used the augccpVTZ basis sets with no f functions. These bond energies were corrected for BSSE (basis set superposition error). For the transition metals we used the TZV(f) basis sets (20).
A big step forward for accurate DFT calculations was the introduction of hybrid methods (3), particularly the B3LYP, which uses the VWN functional III (21) based on correlation of the homogeneous electron gas in the random phase approximation, the Lee–Yang–Parr (17) correlation functional, plus a hybrid exchange functional of three terms: a portion of exact exchange, Slater local exchange, and the nonlocal gradient correction of Becke88. Thus, Becke obtained the hybrid parameters {a _{0}, a_{x}, a_{c} } = {0.20, 0.72, 0.19} (3) from a leastsquares fit to 56 atomization energies, 42 IPs, and 8 proton affinities (PAs) of the G21 set of atoms and molecules (4). B3LYP leads to excellent thermochemistry (0.13 eV MAD) and structures for covalently systems but does not account for London dispersion (all noble gas dimers are predicted unstable).
Following B3LYP, we introduce the extended hybrid functional, denoted as X3LYP: We determined the hybrid parameters {a _{0}, a_{x}, a_{c} } = {0.218, 0.709, 0.129} in X3LYP just as for XLYP. Thus, we normalized the mixing parameters of Eq. 10 and redetermined {a_{x} _{1}, a_{x} _{2}} = {0.765, 0.235} for X3LYP. The F ^{X}(s) function of X3LYP (Fig. 1) agrees with F ^{Gauss}(s) for larger s.
Results and Discussion
We tested the accuracy of XLYP and X3LYP for a broad range of systems and properties not used in fitting the parameters. Table 1 compares the overall performance of 17 different flavors of DFT methods, showing that X3LYP is the best or nearly best for essentially all properties, leading to an acceptable accuracy for each. A few comments follow.
Heats of Formation (ΔH_{f}). To provide a good test of the functionals for the covalent systems, we calculated the heats of formation for the 148 molecules of the extended G2 set (4, 5). X3LYP leads to MAD = 0.12 eV, which is the best DFT result. Other good results include B3LYP (MAD = 0.13 eV) (5), B3PW91 (MAD = 0.15 eV) (5), mPW1PW (MAD = 0.17) (7), and PBE1PBE (MAD = 0.21 eV) (10). All other methods lead to unacceptable errors. Table 1 makes it clear that hybrid methods dramatically improve the thermochemistry over pure GGAs. For example, PBEPBE (MAD: 0.74) compares to PBE1PBE (MAD: 0.21) (8,10); BLYP (MAD = 0.31) versus B3LYP (MAD = 0.13) (5) and XLYP (MAD = 0.33) versus X3LYP (MAD = 0.12). We should point out that the G2 energies quoted here (MAD = 0.07) include the empirical “higherorder correction” to the total energy. Without this correction, the MAD for G21 would have been 1.03 eV for MP2 and 0.50 eV for CCSD[T], partially limited by the use of a smaller basis set (7, 8).
Ionization Potential (IP). IPs (and EAs) are calculated as the total energy differences between the neutral and the corresponding ionic systems. GGAs generally dramatically improve the predictions of IPs of LDA (MAD = 0.67 eV) (19), suggesting that cations are more inhomogeneous than neutral systems. However, including exact exchange is not important for IP. Most GGA methods give a MAD over the 42 cases in G2 (19) lower than 0.2 eV, with O3LYP (16, 22) lowest at 0.139. X3LYP is the second best, giving MAD = 0.154.
Electron Affinity (EA). There is some debate in the literature concerning whether DFT methods are suitable for calculating EAs (10, 23, 24). Because the “selfinteraction error” shifts the Kohn–Sham orbital energies upwards, the anion often has a positive (unstable) highest occupied orbital energy (this cannot ionize because of the finite size of the basis functions). In any case the numerical results lead to predicted EAs with accuracy comparable to the IPs (19). Thus, over the total 25 systems in G2 (19) we obtain MADs of 0.10 eV (B3LYP), 0.13 (PBE1PBE), and 0.09 (X3LYP). Again, inclusion of exact exchange does not improve the performance over the corresponding pure DFT methods.
Proton Affinity (PA). Adding a proton to a neutral molecule leads to a significant change in the density, making it more inhomogeneous. Thus, Table 1 shows that the proton affinities (PAs) over eight systems in G2 (19) are systematically underestimated by LDA (MAD = 0.27 eV). GGAs reduce the LDA errors significantly. B3P86 and B3PW91 (MAD = 0.03 eV) show the best performance, while B3LYP, PBE1PBE, and X3LYP lead to MAD = 0.06, 0.04 (10), and 0.07 eV, respectively.
Total Energies (E_{tot}). We considered the total energies for the first 10 atoms (H to Ne). Comparing to the experimental values (25, 26), LDA (SVWN) leads to significant error (MAD = 6.67 eV), which is dramatically reduced by most GGAs. X3LYP gives good result (MAD = 0.11 eV), which can be compared to B3LYP (MAD = 0.38) and PBE1PBE (MAD = 1.09). O3LYP is the best, leading to MAD = 0.06.
Bonding Properties of RareGas Dimers. Raregas dimers provide the least ambiguous tests for the accuracy in describing van der Waals attraction (London dispersion forces). As summarized in Table 1, we calculated the equilibrium distance and bond energy of He_{2} and Ne_{2} as probes of the accuracy for van der Waals systems (27).
LDA leads to significant overbinding (by 12 times for He_{2}, with a bond too short by 0.6 Å).
Although B88 exchange is very successful in describing covalent bonds, it fails dramatically for London interactions (28, 29). Thus, every DFT method using the B88 exchange functional (pure or hybrid, with and without correlation) leads to unbound raregas dimers, far more unbound than HF.
On the other hand, every DFT method using the PW91 exchange functional severely overbinds raregas dimers, even when the correlation functional is omitted. To ameliorate this problem, Adamo and Barone modified PW91 (7) by fitting the differential exchange energies of raregas dimers to HF values, removing some of the overbinding tendency of PW91. The PBE functional gives an overall good description of raregas dimers.
The van der Waals (London) attraction between raregas atoms arises entirely from electron correlation (the instantaneous interaction between fluctuating dipoles). Thus, HF wavefunctions give a purely repulsive interaction for all interatomic distances and we expect that exchangeonly potentials should also give this behavior. However, we found that the PW91, mPW91, and PBE exchangeonly potentials lead to erratic minima, as do the corresponding hybrid models.
We conclude that X3LYP provides the best description of London forces because coupled with the LYP correlation term its gives a good description of bonding for He_{2} and Ne_{2} while without correlation, the exchangeonly (X or X3) potential is repulsive and close to the Hartree–Fock curve. Because the London forces of He and Ne are described well, we expect a good description of the London forces between electron pairs involving the first 10 atoms of the periodic table.
Note that in describing the dispersion forces, we are concerned with distances near the equilibrium configuration for van der Waals complexes, which tend to be in the range from 2.5 to 4 Å. It is at much larger distances before the neglect of overlap between the atomic distances leads to the classic 1/r ^{6} London potential (30, 31).
Bonding Properties of Water Dimer. Hydrogen bonding plays a critical role in a wide range of chemical and biological phenomena. Consequently, water dimer, the prototypical hydrogen bonded system, has been studied thoroughly experimentally and theoretically (32–34). The floppiness of (H_{2}O)_{2} has made experimental determinations of r _{e} and D _{e} difficult. Microwave measurements of the rotational moments lead to a vibrationally averaged O... O distance of R _{0} = 2.976 Å, from which it was estimated that R _{e} = 2.946 Å (32). The most widely accepted experimental bond energy of D _{e} = 0.23 ± 0.03 eV (33) was based on measurements of the thermal conductivity of water vapor followed by relatively complex interpretation including an estimate of the zeropoint energy calculated at the HF/4–21G level (33).
Fortunately, highly accurate values for the equilibrium geometry and dissociation energy of (H_{2}O)_{2} have been determined computationally (34). These ab initio CCSD(T) (Full) studies used a 275 basis function interaction optimized basis set extrapolated to infinity, leading to r _{e}(O... O) = 2.912 ± 0.005 Å and D _{e} = 0.218 ± 0.004 eV (34). We conclude that the ab initio values (34) provide the most reliable data.
All DFT methods that include GGA give H bond energies within 0.04 eV of D _{e} = 0.218 eV and bond distances within 0.04 Å of r _{e} = 2.912 Å (34). This may seem adequate, but for biological applications, greater accuracy is highly desired. The best accuracy is given by PBE1PBE and X3LYP, which lead to binding energy errors of 0.002 eV and bond distance errors of 0.016 Å (PBE1PBE) and 0.004 Å (X3LYP). We should emphasize here that (H_{2}O)_{2} was not included in the training set for these functionals, indicating that physics of hydrogen bonding is correctly described.
Electrostatics. Electrostatic interactions are quite important in the bonding of ligands to proteins and DNA in determining differential solvation also important in binding. To assess the accuracy for charge distributions from various DFT methods, we calculated the dipole moments for 70 molecules for which accurate values are available. HF generally overestimates the magnitude of dipole moments (MAD = 0.242 Debye), but all DFT methods lead to substantial improvement (e.g., MAD errors are: B3LYP, 0.083; PW1PW, 0.088; PBE1PBE, 0.084; and X3LYP, 0.088). These results were not used in optimizing the parameters.
In addition, the polarizabilities of atoms and molecules are expected to be important in electrostatic and solvation interactions. For the seven molecules for which we could find accurate experimental polarizabilities, the accuracies (MAD in au) for various DFT methods are as follows: B3LYP, 0.16; PW1PW, 0.35; PBE1PBE, 0.32; and X3LYP, 0.19. These results were not used in optimizing the parameters.
Spin States of Transition Metals. Many important biological systems involve metal atoms in various spin and oxidation states. For example, the various steps of detoxifying foreign molecules by cytochrome P450s involves highspin Fe^{III}, highspin Fe^{II}, and lowspin Fe^{II}. Consequently, it is important that the excitation energies and IPs be properly described (35, 36). We considered the nine transition metal atoms Sc–Cu and calculated the d^{n}s ^{2} to d^{n} ^{+1} s ^{1} excitation energy of the neutral atom and the d^{n}s ^{1} to d ^{n+1} s ^{0} (n = 1–9) excitation energy of the cation for various DFT functionals. Table 1 shows that X3LYP does quite well with a MAD of 0.22 eV, as compared to B3LYP (0.25) (37) and PBE1PBE (0.30). These results were not used in optimizing the parameters.
Summary
We deduce the form for the exact exchange energy density to describe a density decaying as a Gaussian at long range. We find that F ^{Gauss}(s) lies between F ^{B88}(s) and F ^{PW91}(s), serving as the basis for the extended functional, F ^{X}, which is described a linear combination of F ^{PW91} (with a sound physical basis) with F ^{B88} (which in B3LYP does best for thermochemistry but badly for London forces). Optimizing the four parameters [using the atomization energies of a small set (37) of diatomic and triatomic molecules along with atomic energies or IPs and EAs of the first 18 atoms] leads to the X3LYP functional with excellent accuracy for thermochemistry (heats of formation, IPs, EAs, proton affinities, and total atomic energies), a good description of London dispersion, an excellent description of hydrogen bond interactions, and excellent energetics for the spin states of transition metals. Thus, we recommend the use of X3LYP for application to a wide range of chemical, biological, and materials systems. X3LYP should also form a good starting point for continuing attempts to develop improved functionals for exchange and for correlation.
Acknowledgments
We thank Drs. Y. X. Cao, D. Braden, and Jason Perry at Schrödinger (Portland, OR) and Q. E. Zhang and R. P. Muller at the Materials and Process Stimulation Center (MSC) for technical support in using and modifying jaguar for these calculations. This work was supported by the Department of Energy [Accelerated Strategic Computing Initiative (ASCI)–Academic Strategic Alliances Program] with partial support from National Science Foundation Grants CHE9985574 and CTS0132002, National Natural Science Foundation of China Grant 20021002, National Natural Science Foundation of Fujian Grant 2002F010, Ministry of Science and Technology of China Grant 2001CB610506, and the Teaching and Research Award Program for Outstanding Young Teachers of the Ministry of Education of China. The computation facilities of the MSC used in these studies have been supported by grants from Army Research Office (ARO)–Defense University Research Instrumentation Program (DURIP), Office of Naval Research (ONR)–DURIP, the National Science Foundation (Major Research Instrumentation, Chemistry), and IBM–Shared University Research. In addition, the MSC is supported by grants from Department of Energy ASCI, ARO–Multidisciplinary University Research Initiative (MURI), ARO–Defense Advanced Research Planning Agency, ONR–MURI, National Institutes of Health, ONR, General Motors, ChevronTexaco, Seiko–Epson, Beckman Institute, and Asahi Kasei.
Footnotes

↵ ‡ To whom correspondence should be addressed. Email: wag{at}wag.caltech.edu.

Abbreviations: B3LYP, Becke threeparameter hybrid functional combined with Lee–Yang–Parr correlation functional; DFT, density functional theory; EA, electron affinity; GGA, generalized gradient approximation; IP, ionization potential; LDA, local density approximation; MAD, mean absolute deviation; X3LYP, extended hybrid functional combined with Lee–Yang–Parr correlation functional.
 Copyright © 2004, The National Academy of Sciences
References

↵
Chong, D. P., ed. (1997) Recent Advances in Density Functional Methods (World Scientific, Singapore), Pts. I and II.

↵
Barone, V. & Bencini, A., eds. (1999) Recent Advances in Density Functional Methods (World Scientific, Singapore), Pt. III.
 ↵
 ↵
 ↵

↵
Perdew, J. P. (1991) in Electronic Structure of Solids'91, eds. Ziesche, P. & Eschrig, H. (Akademie, Berlin), p. 11.
 ↵
 ↵
 ↵
 ↵
 ↵
 ↵
 ↵
 ↵

↵
Lieb, E. H. & Oxford, S. (1981) Int. J. Quantum Chem. 19 , 427439.
 ↵
 ↵
 ↵
 ↵
 ↵

↵
Vosko, S. H., Wilk, L. & Nusair, M. (1980) Can. J. Phys. 58 , 12001211.
 ↵
 ↵
 ↵
 ↵
 ↵
 ↵
 ↵
 ↵
 ↵
 ↵
 ↵
 ↵
 ↵
 ↵
 ↵

↵
Koch, W. & Hertwig, R. H. (1998) in Encyclopedia of Computational Chemistry, ed. von Ragué Schleyer, P. (Wiley, Chichester, U.K.).
 ↵
Citation Manager Formats
More Articles of This Classification
Physical Sciences
Related Content
 No related articles found.