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 nuclearnuclear repulsion E_{NN}, the nuclearelectron attraction E_{v}, and the classical electronelectron Coulomb repulsion E_{coul} energies are the same as those used in HartreeFock theory. The kinetic energy of the electrons E_{T} as well as the nonclassical electronelectron exchange energies E_{exch} are, however, different from those used in HartreeFock theory. The last term E_{corr} describes the correlated movement of electrons of different spin and is not accounted for in HartreeFock theory. Due to these differences, the exchange energies calculated exactly in HartreeFock 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 HartreeFock exchange energies are used.
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 methanol1yl radical using the SVWN functional and the 631G(d) basis set is as follows:
#SVWN/631G(d) freq
SVWN/631G(d) freq methanol1yl 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 VWNpaper (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 gradientcorrected exchange functionals available in Gaussian include (given with their abbreviations in parenthesis):
and the gradientcorrected correlation functionals are
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 methanol1yl radical is:
#P BLYP/631G(d) freq
BLYP/631G(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 (HartreeFocklike) manner with those obtained from DFT methods in order to improve performance. Frequently used methods are:
The BeckehalfandhalfLYP method is used with the BHandHLYP keyword, the very popular Becke3LYP 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 methanol1yl radical is:
#P B3LYP/631G(d) freq
Becke3LYP/631G(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, Turbomole, Jaguar, QChem) 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) Userdefined models
Calculation of reaction and activation energies in either radical reactions or transition metal mediated reactions often shows that the mixing ratio of HF and DFTexchange is not optimal for a given purpose. The definition of userdefined 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 ratios of the various exchange and correlation energy terms. For this latter step Gaussian uses the general formula:
E_{XC} = P_{2}E_{X}^{HF} + P_{1}(P_{4}E_{X}^{Slater} + P_{3}dE_{X}^{nonlocal}) + P_{6}E_{C}^{local} + P_{5}dE_{C}^{nonlocal}
The values for the parameters P_{1}  P_{6} are provided through the IOP switches IOP(5/45=xxxxyyyy), IOP(5/46=xxxxyyyy), and IOP(5/47=xxxxyyyy) with P_{1} = xxxx/1000 from IOP(5/45) 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) developed to optimize reaction and activation energies of free radical reactions. This method is very similar to the MPW1PW91 hybrid functional, but uses a mixing ratio of exact and DFT exchange energy terms of 0.428 : 0.572. An input file for the methanol1yl radical frequency calculation at the MPW1K/631G(d) level is:
#P MPWPW91/631G(d) freq IOP(5/45=10000428) IOP(5/46=05720572) IOP(5/47=10001000) MPW1K/631G(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_{2}. 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 ScaHFX= 0.4280
ScaDFX= 0.5720 0.5720 1.0000 1.0000
The value given after ScaHFX corresponds to P_{2} while the four values listed after ScaDFX are P_{3}  P_{6}. Please observe that theses lines are only printed when using the extended output option of Gaussian. Before using this facility with a particular version of Gaussian, please make sure that your IOPdirectives are read correctly. In some versions of Gaussian, the options P_{1}  P_{6} are set with IOP(3/76)  IOP(3/78) instead of IOP(5/45)  IOP(5/47). In this latter case, the values of P1  P6 are given in fivedigit format. The keyword line for the PBE0 frequency calculation thus becomes:
#P PBEPBE/631G(d) freq
IOP(3/76=1000002500) IOP(3/77=0750007500) IOP(3/78=1000010000)
Absolute and relative CPUtimes for DFT and HartreeFockcalculations 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 methanol1yl radical described above with Gaussian 03, Rev. B.03. The platforms used are the CIProom 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 631G(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/631G(d)  31  3690  1455  1036 
UBLYP/631G(d)  69  3599  1463  1042 
UB3LYP/631G(d)  71  3761  1508  1072 
UMPW1K/631G(d)  73  3950  1552  1102 
UHF/631G(d)  11  4123  1627  1156 
obsvd.    3650  1459  1056 
Despite all platformdependend variations, it is clear that hybrid DFT calculations are similarly laborious as GGA calculations, while local DFT methods are somewhat and HFcalculations 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 (see our earlier discussion of this point). 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 HartreeFock exchange contributions (varying between 0 and 100% from the second entry to the bottom in Table 1).
5) Further technical considerations
Evaluation of the exchangecorrelation energies of all density functional methods implemented in Gaussian involves a gridbased numerical integration step. The computational effort required for this step strongly depends on the selected grid size. The larger the number of integration points per atom, the higher is the computational cost and the better the numerical accuracy of the calculation. Several predefined grids are available in Gaussian, which can be used through the keyword
int=(gridsize) or int=(Grid=gridsize)
Alternatively, the grid size can also be specified with IOP settings, using either IOP(5/44=gridsize) or IOP(3/75=gridsize), depending on the version of Gaussian. The following are acceptable values for gridsize:
gridsize  corresponding IOP setting 
(radial shells * angular points) 
CPU [sec]  E_{tot}  frequencies (cm^{1}) 

CoarseGrid  2  (35*110), pruned  33.8  115.0520127  110.9325 0.0015 0.0009 0.0008 24.6825 59.6546 
SG1Grid  1  (50*194), pruned  45.1  115.0520371  30.5802 0.0015 0.0015 0.0016 10.5877 29.2740 
finegrid  4  (75*302), pruned  77.1  115.0520326  20.7294 8.3913 0.0009 0.0009 0.0014 11.1426 
ultrafine  5  (99*590), pruned  192.7  115.0520324  0.0013 0.0013 0.0011 2.6708 4.0642 14.0308 
mmm,nnn  mmmnnn  (mmm*nnn), not pruned       
The CPU times and total energies listed here (in sec.) together with the grid specifications are those for the UB3LYP/631G(d) frequency calculation on the methanol radical mentioned before on the CIP room Pentium 4 LINUX PCs. It is clear from these results that neither the total electronic energy nor the lowfrequency vibrations converge quickly with increasing grid size. Thus, for each application it must be decided, whether the substantial computational effort for a very high grid size is justified. It is also important to note that the energy variation with grid size is significant enough that all calculations for stationary points on the same potential energy surface must be performed with the same grid in order to be comparable!
In addition to the grid selection for integral evaluation using the int=(gridsize) keyword, a particular grid can also be selected for the calculation of second derivatives using the CPHF=gridsize keyword. The available options for gridsize are identical in both cases. If no explicit grid is specified through the CPHF keyword, the following relations exist between the grids used for integral and second derivative calculations:
gridsize used for integral evaluation 
grid used for second derivatives 

CoarseGrid  CoarseGrid 
SG1Grid  CoarseGrid 
finegrid  CoarseGrid 
ultrafine  SG1Grid 
(mmm*nnn)  (mmm*nnn) 
These combinations are usually appropriate for most systems. The default grid size for integral evaluation in the most recent versions of Gaussian is finegrid. Better grids are required for kinetic isotope effect calculations, the optimization of stationary points on very flat potential energy surfaces, and IRCfollowing on flat potential energy surfaces. Increasing the CPHF grid size over what is specified per default through the integral grid selection is rarely helpful. Smaller grids can be quite useful to speed up geometry optimizations of large systems, especially when still far away from the next stationary point.
6) Treating open shell systems
For practical reasons the DFT methods also use a DFT "wavefunction" constructed as an antisymmetrized product of one electron functions (now termed KohnSham orbitals). As in HartreeFock theory, electrons of unlike spin can either use the same spatial orbitals (restricted DFT methods) or different spatial orbitals (unrestricted DFT methods). Specification of a restricted open shell wavefunction for the Becke3LYP hybrid functional requires the ROB3LYP keyword, while the unrestricted version is specified using UB3LYP. If no explicit choice is made in the input file, open shell systems are treated with the unrestricted approach by default. It should be noted that unrestricted DFT methods (GGA or hybrid functionals alike) are much less affected by the phenomenon of spin contamination than HartreeFock theory and also give much better predictions of ESR spectra of organic radicals. The UB3LYP/631G(d) calculation of the allyl radical is achieved with the following input file:
#UBecke3LYP/631G(d) scf=(tight)
UB3LYP/631G(d) ally radical (C2v)
0 2
H
C,1,r2
C,2,r3,1,a3
C,2,r3,1,a3,3,180.
H,3,r5,2,a5,1,0.
H,4,r5,2,a5,1,0.
H,3,r7,2,a7,1,180.
H,4,r7,2,a7,1,180.
r2=1.09054601
r3=1.38611149
r5=1.08492921
r7=1.08690606
a3=117.45043254
a5=121.64506508
a7=121.12827199

Comparison of the expectation value of the <S^{2}> operator calculated at the UHF/631G(d) level (0.9729) with that obtained at UB3LYP/631G(d) (0.7818) shows that the spin contamination problem is much reduced with DFT methods. From the construction principle of the hybrid DFT methods it is clear that a larger HFexchange component will also give rise to larger amounts of spin contamination. Using the allyl radical at the UB3LYP/631G(d) structure as an example, the <S^{2}> value amounts to 0.7650 at the BLYP/631G(d) level and to 0.8228 at the BHandHLYP/631G(d) level. The precise values of the expectation value of the <S^{2}> operator will also depend to some extend on the basis set used.
Another important point to consider in DFT calculations concerns the self interaction error resulting from an artificial interaction of the electron density with itself. While the definition of the exchange and Coulomb energies in HartreeFock theory leads to a perfect cancellation of these self interaction terms, the cancellation is incomplete in DFT methods. As nicely explained by Koch and Holthausen in their monograph on DFT methods, this can be easily demonstrated using the hydrogen atom. The exact energy for this oneelectron system is 0.50000 Hartree. As there is no second electron in the system, no correlation energies have to be calculated and HartreeFock theory is able to provide the exact answer in the complete basis set limit. Using a basis set considered "large" for molecular systems such as the ccpVQZ basis set, a value of 0.499946 Hartree is indeed obtained. Due to the parameterized exchange and correlation functionals, however, much more negative values can be obtained in particular with hybrid DFT functionals. Using the same basis set, the UB3LYP functional predicts a value of 0.502346 Hartree!!