Density Functional Theory

In density functional theory (DFT) the energy of a system is given as a sum of six components:

E_DFT = E_NN + E_T + E_v + E_coul + E_exch + E_corr

The definitions for the nuclear-nuclear repulsion E_NN, the nuclear-electron attraction E_v, and the classical electron-electron Coulomb repulsion E_coul energies are the same as those used in Hartree-Fock theory. The kinetic energy of the electrons E_T as well as the non-classical electron-electron exchange energies E_exch are, however, different from those used in Hartree-Fock theory. The last term E_corr describes the correlated movement of electrons of different spin and is not accounted for in Hartree-Fock theory. Due to these differences, the exchange energies calculated exactly in Hartree-Fock theory cannot be used in density functional theory.

Various approaches exist to calculate the exchange and correlation energy terms in DFT methods. These approaches differ in using either only the electron density (local methods) or the electron density as well as its gradients (gradient corrected methods or generalized gradient approximation, GGA). Aside from these "pure" DFT methods, another group of hybrid functionals exists, in which mixtures of DFT and Hartree-Fock exchange energies are used. This class has recently been extended to include double hybrid functionals (DHDF) in which mixtures of exchange energies are combined with mixtures of correlation energies. The latter are typically obtained using a GGA correlation functional and an orbital-based perturbation theory method such as MP2.


1) Local methods
The only local exchange functional available in Gaussian is the Slater functional. Combination with the local VWN correlation functional by Vosko, Wilk and Nusair gives the Local Density Approximation (LDA) method. For open shell systems, using unrestricted wavefunctions, this is also referred to as the Local Spin Density Approximation (LSDA). This method is used in Gaussian with either the LSDA or SVWN keyword. A sample input for the frequency calculation for the methanol-1-yl radical using the SVWN functional and the 6-31G(d) basis set is as follows:
 

#SVWN/6-31G(d) freq SVWN/6-31G(d) freq methanol-1-yl radical 0 2 C1 O2 1 r2 H3 2 r3 1 a3 H4 1 r4 2 a4 3 d4 H5 1 r5 2 a5 4 d5 r2=1.35454381 r3=0.97675734 r4=1.09713107 r5=1.09208721 a3=108.95266528 a4=119.30335589 a5=113.0329959 d4=-26.58949552 d5=-149.77421707

Please observe that the acrynom VWN is used to indicate different functionals in different software packages. While in Gaussian this keyword refers to functional III described in the VWN-paper (S. J. Vosko, L. Wilk, M. Nusair, Can. J. Phys. 1980, 58, 1200), most other packages refer to functional V in the same paper. This latter functional can also be used in Gaussian with the VWN5 keyword.


2) Gradient corrected methods
The gradient-corrected exchange functionals available in Gaussian include (given with their abbreviations in parenthesis):
 

• Becke88 (B)
• Perdew-Wang (PW91)
• Modified Perdew-Wang by Barone and Adamo (MPW)
• Gill96 (G96)
• Perdew-Burke-Ernzerhof (PBE)

and the gradient-corrected correlation functionals are
 

• LYP by Lee, Yang, and Parr (LYP)
• Perdew-Wang (PW91)
• Perdew 86 (P86)
• Becke96 (B96)
• Perdew-Burke-Ernzerhof (PBE)

In all cases, the names of these functionals refer to their respective authors and the year of publication. All combinations of exchange and correlation functionals are possible, the keywords being composed of the acronyms for the two functionals. The frequently used BLYP method, for example, combines Becke's 1988 exchange functional with the correlation functional by Lee, Yang, and Parr. An input file for the BLYP frequency calculation of the methanol-1-yl radical is:
 

#P BLYP/6-31G(d) freq BLYP/6-31G(d) opt min methanol radical 0 2 C1 O2 1 r2 H3 2 r3 1 a3 H4 1 r4 2 a4 3 d4 H5 1 r5 2 a5 4 d5 r2=1.38575119 r3=0.98034945 r4=1.09658704 r5=1.09123991 a3=108.13230593 a4=118.45265645 a5=112.14316572 d4=-29.40173022 d5=-145.83450795

Another frequently used GGA functional is BP86 composed of the Becke 1988 exchange functional and the Perdew 86 correlation functional. The PW91 functional combines exchange and correlation functionals developed by the same authors in 1991. The keywords used in Gaussian for a particular GGA functional are combinations of the acronyms for exchange and correlation functionals. Use of the PW91 exchange and correlation functionals thus requires the keyword PW91PW91.


3) Hybrid functionals
The basic idea behind the hybrid functionals is to mix exchange energies calculated in an exact (Hartree-Fock-like) manner with those obtained from DFT methods in order to improve performance. Frequently used methods are:

• Becke-Half-and-Half-LYP (BHandHLYP) uses a 1:1 mixture of DFT and exact exchange energies:
  E_XC = 0.5*E_X(HF) + 0.5*E_X(B88) + E_C(LYP)
   
• Becke-3-LYP (B3LYP) uses a different mixing scheme involving three mixing parameters:
  E_XC = 0.2*E_X(HF) + 0.8*E_X(LSDA) + 0.72*DE_X(B88) + 0.81*E_C(LYP) + 0.19*E_C(VWN)
  In this latter case, the B88 gradient correction to the local LSDA exchange energies carries its own scaling factor of 0.72. The three scaling factors have been derived through fitting the parameters to a set of thermochemical data (G1 set).
   
• PBE0 uses a 1:3 mixture of DFT and exact exchange energies:
  E_XC = 0.25*E_X(HF) + 0.75*E_X(PBE) + E_C(mPW91)
   
• B98 (Becke98) is based on a 10-parameter equation scaling components of exact exchange, GGA exchange and GGA correlation. All parameters have been optimized simultaneously to fit thermochemica data collected in the extended G2 data set:
  E_XC = 0.2198*E_X(HF) + E_X + E_C
   

The Becke-half-and-half-LYP method is used with the BHandHLYP keyword, the very popular Becke-3-LYP method is used with the B3LYP keyword, and the PBE0 method with the PBE1PBE keyword. A sample input file for the B3LYP frequency calculation on the methanol-1-yl radical is:
 

#P B3LYP/6-31G(d) freq Becke-3-LYP/6-31G(d) opt min methanol radical 0 2 C1 O2 1 r2 H3 2 r3 1 a3 H4 1 r4 2 a4 3 d4 H5 1 r5 2 a5 4 d5 r2=1.37007444 r3=0.96903862 r4=1.08886325 r5=1.08381987 a3=108.87657926 a4=118.49962766 a5=112.61039767 d4=-29.19863769 d5=-146.73450089

It is important to note that implementations of the B3LYP and BHLYP functionals in different programs (Gaussian, Jaguar, MOLPRO, ORCA, Q-Chem, Turbomole)) vary somewhat. Other hybrid functionals available in Gaussian are: B3P86 uses the P86 correlation functional instead of LYP, but retains the three parameters derived for B3LYP; B3PW91 uses the PW91 correlation functional instead of LYP, but retains the three parameters derived for B3LYP; B1B96, a one parameter functional developed by Becke using the Becke 96 correlation functional; MPW1PW91 developed by Barone and Adamo using a modified version of the PW91 exchange functional in combination with the original PW91 correlation functional and a mixing ratio of exact and DFT exchange of 0.25 : 0.75.


4) User-defined models
Calculations of reaction and activation energies in radical reactions or transition metal mediated reactions often show that the mixing ratio of HF- and DFT-exchange is not optimal for a given purpose. The definition of user-defined models tailored to a certain purpose is therefore sometimes necessary. This can be achieved by first choosing gradient corrected exchange and correlation functionals from the list of available choices and then specifying the mixing ratio of the various exchange and correlation energy terms. For this latter step Gaussian uses the general formula:

E_XC = P₂E_X^HF + P₁(P₄E_X^Slater + P₃dE_X^non-local) + P₆E_C^local + P₅dE_C^non-local

The values for the parameters P₁ - P₆ are provided through the IOP switches IOP(3/76=xxxxxyyyyy), IOP(3/77=xxxxxyyyyy), and IOP(3/78=xxxxxyyyyy) with P₁ = xxxxx/10000 from IOP(3/76) etc. The use of this facility will be illustrated with the MPW1K functional by Truhlar and coworkers (B. J. Lynch, P. L. Fast, M. Harris, D. G. Truhlar, J. Phys. Chem. A 2000, 104, 4811-4815) developed to optimize reaction and activation energies of free radical reactions. This method is very similar to the mPW1PW91 hybrid functional developed by Adamo and Barone (C. Adamo, V. Barone, J. Chem. Phys. 1998, 108, 664-675) , but uses a mixing ratio of exact and DFT exchange energy terms of 0.428 : 0.572. An input file for the methanol-1-yl radical frequency calculation at the MPW1K/6-31G(d) level is:
 

#P MPWPW91/6-31G(d) freq IOp(3/76=1000004280) IOp(3/77=0572005720) IOp(3/78=1000010000) MPW1K/6-31G(d) opt min methanol radical 0 2 C1 O2 1 r2 H3 2 r3 1 a3 H4 1 r4 2 a4 3 d4 H5 1 r5 2 a5 4 d5 r2=1.35355619 r3=0.95690342 r4=1.08170148 r5=1.07685321 a3=109.38001654 a4=118.56662303 a5=112.99211737 d4=-29.64928765 d5=-147.42512309

The gradient corrected MPWPW91 functional is chosen here as the basis of the hybrid method and the amount of exact exchange is specified through parameter P₂. The composition of the functional is also described in the output file in a compact format:
 

IExCor= 908 DFT=T Ex=mPW+HF Corr=PW91 ExCW=0 ScaHFX=  0.428000
ScaDFX=  0.572000  0.572000  1.000000  1.000000

The value given after ScaHFX corresponds to P₂ while the four values listed after ScaDFX are P₃ - P₆. Please observe that theses lines are only printed when using the extended output option of Gaussian. From the general formula for user-defined methods it is apparent that the parameter P₁ can be used to specify DFT exchange contributions in methods with identical parameters P₃ and P₄. The above specification of the MPW1K functional is thus exactly identical to the following:
 

#P MPWPW91/6-31G(d) freq IOp(3/76=0572004280)

MPW1K/6-31G(d) opt min methanol radical

0 2

Before using this facility with a particular version of Gaussian, please make sure that your IOP directives are read correctly. In some older versions of Gaussian, the options P₁ - P₆ are set with IOP(5/45) - IOP(5/47) instead of IOP(3/76) - IOP(3/78). In this latter case, the values of P1 - P6 are given in four-digit format. The keyword line for the PBE0 frequency calculation thus becomes:
 

#P PBEPBE/6-31G(d) freq
 IOP(5/45=10000250) IOP(5/46=07500750) IOP(5/47=10001000)

Absolute and relative CPU-times for DFT- and Hartree-Fock-calculations depend dramatically on the size of the system under investigation, the computer hardware and the operating system used. The following table contains CPU times (in seconds) for all frequency calculations on the methanol-1-yl radical described above with Gaussian 03, Rev. B.03. The platforms used are the CIP-room Pentium 4/2 GHz computers under LINUX. In all cases the frequency calculation has been performed on the structure optimized at the given theoretical level. The 6-31G(d) basis set has been used in all cases.
 

method	CPU (P4/2GHz)  [sec]	vib(1)  [cm-1]	vib(2)  [cm-1]	vib(3)  [cm-1]
USVWN/6-31G(d)	31	3690	1455	1036
UBLYP/6-31G(d)	69	3599	1463	1042
UB3LYP/6-31G(d)	71	3761	1508	1072
UMPW1K/6-31G(d)	73	3950	1552	1102
UBHandHLYP/6-31G(d)	74	3946	1565	1110
UHF/6-31G(d)	11	4123	1627	1156
obsvd.	-	3650	1459	1056

Despite all platform-dependent variations, it is clear that hybrid DFT calculations are similarly laborious as GGA calculations, while local DFT methods are somewhat and UHF calculations much cheaper. In matching calculated (harmonic) gas phase vibrational frequencies with experimentally observed anharmonic values measured in a matrix (R. D. Johnson III, J. W. Hudjens, J. Phys. Chem., 1996, 100, 19874), please remember that uniform scaling is often performed to improve correspondence. The values reported here are unscaled. A general trend visible in this example is the increase of vibrational frequency wavenumbers with an increase in the amount of Hartree-Fock exchange contributions (varying between 0 and 100% from the second entry to the bottom in Table 1).