In recent years, Density Functional Theory [21, 22, 23, 24] has emerged as an accurate alternative firstprinciples approach to quantum mechanical molecular investigations. DFT currently accounts for approximately 90% of all quantum chemical calculations being performed, not only because of its proven chemical accuracy, but also because of its relatively cheap computational expense. These two features suggest that DFT is likely to remain a leading method in the quantum chemist’s toolkit well into the future. QChem contains fast, efficient and accurate algorithms for all popular density functional theories, which make calculations on quite large molecules possible and practical.
DFT is primarily a theory of electronic ground state structures based on the electron density, , as opposed to the manyelectron wavefunction There are a number of distinct similarities and differences to traditional wavefunction approaches and modern DFT methodologies. Firstly, the essential building blocks of the many electron wavefunction are singleelectron orbitals are directly analogous to the KohnSham (see below) orbitals in the current DFT framework. Secondly, both the electron density and the manyelectron wavefunction tend to be constructed via a SCF approach that requires the construction of matrix elements which are remarkably and conveniently very similar.
However, traditional approaches using the many electron wavefunction as a foundation must resort to a postSCF calculation (Chapter 5) to incorporate correlation effects, whereas DFT approaches do not. PostSCF methods, such as perturbation theory or coupled cluster theory are extremely expensive relative to the SCF procedure. On the other hand, the DFT approach is, in principle, exact, but in practice relies on modeling the unknown exact exchange correlation energy functional. While more accurate forms of such functionals are constantly being developed, there is no systematic way to improve the functional to achieve an arbitrary level of accuracy. Thus, the traditional approaches offer the possibility of achieving an arbitrary level of accuracy, but can be computationally demanding, whereas DFT approaches offer a practical route but the theory is currently incomplete.
The Density Functional Theory by Hohenberg, Kohn and Sham [25, 26] stems from the original work of Dirac [27], who found that the exchange energy of a uniform electron gas may be calculated exactly, knowing only the charge density. However, while the more traditional DFT constitutes a direct approach and the necessary equations contain only the electron density, difficulties associated with the kinetic energy functional obstructed the extension of DFT to anything more than a crude level of approximation. Kohn and Sham developed an indirect approach to the kinetic energy functional which transformed DFT into a practical tool for quantum chemical calculations.
Within the KohnSham formalism [26], the ground state electronic energy, , can be written as
(4.30) 
where is the kinetic energy, is the electron–nuclear interaction energy, is the Coulomb selfinteraction of the electron density and is the exchangecorrelation energy. Adopting an unrestricted format, the alpha and beta total electron densities can be written as
(4.31) 
where and are the number of alpha and beta electron respectively and, are the KohnSham orbitals. Thus, the total electron density is
(4.32) 
Within a finite basis set [28], the density is represented by
(4.33) 
The components of Eq. (4.28) can now be written as
(4.34)  
(4.35)  
(4.36)  
(4.37) 
Minimizing with respect to the unknown KohnSham orbital coefficients yields a set of matrix equations exactly analogous to the UHF case
(4.38)  
(4.39) 
where the Fock matrix elements are generalized to
(4.40)  
(4.41) 
where and are the exchangecorrelation parts of the Fock matrices dependent on the exchangecorrelation functional used. The PopleNesbet equations are obtained simply by allowing
(4.42) 
and similarly for the beta equation. Thus, the density and energy are obtained in a manner analogous to that for the HartreeFock method. Initial guesses are made for the MO coefficients and an iterative process applied until self consistency is obtained.
There are an increasing number of exchange and correlation functionals and hybrid DFT methods available to the quantum chemist, many of which are very effective. In short, there are nowadays five basic working types of functionals (five rungs on the Perdew’s “Jacob‘s Ladder"): those based on the local spin density approximation (LSDA) are on the first rung, those based on generalized gradient approximations (GGA) are on the second rung. Functionals that include not only density gradient corrections (as in the GGA functionals), but also a dependence on the electron kinetic energy density and/or the Laplacian of the electron density, occupy the third rung of the Jacob‘s Ladder and are known as “metaGGA" functionals. The latter lead to a systematic, and often substantial improvement over GGA for thermochemistry and reaction kinetics. Among the metaGGA functionals, a particular attention deserve the VSXC functional [29], the functional of Becke and Roussel for exchange [30], and for correlation [31] (the BR89B94 metaGGA combination [31]). The latter functional did not receive enough popularity until recently, mainly because it was not representable in an analytic form. In QChem, BR89B94 is implemented now selfconsistently in a fully analytic form, based on the recent work [32]. The one and only nonempirical metaGGA functional, the TPSS functional [33], was also implemented recently in QChem [34]. Each of the above mentioned “pure” functionals can be combined with a fraction of exact (HartreeFock) nonlocal exchange energy replacing a similar fraction from the DFT local exchange energy. When a nonzero amount of HartreeFock exchange is used (less than a 100%), the corresponding functional is a hybrid extension (a global hybrid) of the parent “pure” functional. In most cases a hybrid functional would have one or more (up to 21 so far) linear mixing parameters that are fitted to experimental data. An exception is the hybrid extension of the TPSS metaGGA functional, the nonempirical TPSSh scheme, which is also implemented now in QChem [34].
The forth rung of functionals (“hyperGGA” functionals) involve occupied KohnSham orbitals as additional nonlocal variables [35, 36, 37, 38]. This helps tremendously in describing cases of strong inhomogeneity and strong nondynamic correlation, that are evasive for global hybrids at GGA and metaGGA levels of the theory. The success is mainly due to one novel feature of these functionals: they incorporate a 100% of exact (or HF) exchange combined with a hyperGGA model correlation. Employing a 100% of exact exchange has been a long standing dream in DFT, but most previous attempts were unsuccessful. The correlation models used in the hyperGGA schemes B05 [35] and PSTS [38], properly compensate the spuriously high nonlocality of the exact exchange hole, so that cases of strong nondynamic correlation become treatable.
In addition to some GGA and metaGGA variables, the B05 scheme employs a new functional variable, namely, the exactexchange energy density:
(4.43) 
where
(4.44) 
This new variable enters the correlation energy component in a rather sophisticated nonlinear manner [35]: This presents a huge challenge for the practical implementation of such functionals, since they require a HartreeFock type of calculation at each grid point, which renders the task impractical. Significant progress in implementing efficiently the B05 functional was reported only recently [39, 40]. This new implementation achieves a speedup of the B05 calculations by a factor of 100 based on resolutionofidentity (RI) technique (the RIB05 scheme) and analytical interpolations. Using this methodology, the PSTS hyperGGA was also implemented in QChem more recently [34]. For the time being only singlepoint SCF calculations are available for RIB05 and RIPSTS (the energy gradient will be available soon).
In contrast to B05 and PSTS, the forthrung functional MCY employs a 100% global exact exchange, not only as a separate energy component of the functional, but also as a nonlinear variable used the MCY correlation energy expression [36, 37]. Since this variable is the same at each grid point, it has to be calculated only once per SCF iteration. The form of the MCY correlation functional is deduced from known adiabatic connection and coordinate scaling relationships which, together with a few fitting parameters, provides a good correlation match to the exact exchange. The MCY functional [36] in its MCY2 version [37] is now implemented in QChem, as described in Ref. [34].
The fifthrung functionals include not only occupied KohnSham orbitals, but also unoccupied orbitals, which improves further the quality of the exchangecorrelation energy. The practical application so far of these consists of adding empirically a small fraction of correlation energy obtained from MP2like postSCF calculation [41, 42]. Such functionals are known as “doublehybrids”. A more detailed description of some these as implemented in QChem is given in Sections 4.3.9 and 4.3.4.3. Finally, the socalled rangeseparated (or longrange corrected, LRC) functionals that employ exact exchange for the longrange part of the functional are showing excellent performance and considerable promise (see Section 4.3.4). In addition, many of the functionals can be augmented by an empirical dispersion correction, “D” (see Section 4.3.6).
In summary, QChem includes the following exchange and correlation functionals:
LSDA functionals:
SlaterDirac (Exchange) [27]
VoksoWilkNusair (Correlation) [43]
PerdewZunger (Correlation) [44]
Wigner (Correlation) [45]
PerdewWang 92 (Correlation) [46]
ProynovKong 2009 (Correlation) [47]
GGA functionals:
Becke86 (Exchange) [48]
Becke88 (Exchange) [49]
PW86 (Exchange) [50]
refit PW86 (Exchange) [51]
Gill96 (Exchange) [52]
GilbertGill99 (Exchange) [53]
LeeYangParr (Correlation) [54]
Perdew86 (Correlation) [55]
GGA91 (Exchange and correlation) [56]
mPW1PW91 (Exchange and Correlation) [57]
mPW1PBE (Exchange and Correlation)
mPW1LYP (Exchange and Correlation)
revPBE (Exchange) [60]
LB94 scheme (Exchange) [61]
LFA schemes (Exchange) [62]
AK13 (Exchange) [63]
PBE0 (25% HartreeFock exchange + 75% PBE exchange + 100% PBE correlation) [64]
PBE50 (50% HartreeFock exchange + 50% PBE exchange + 100% PBE correlation)
B3LYP (Exchange and correlation within a hybrid scheme) [65]
B3PW91 (B3 Exchange + PW91 correlation)
B3P86 (B3 Exchange + PW86 correlation)
B5050LYP (50% HartreeFock exchange + 5% Slater exchange + 42% Becke exchange + 100% LYP correlation) [66]
BHHLYP (50% HartreeFock exchange + 50% Becke exchange + 100% LYP correlation) [65]
O3LYP (Exchange and correlation) [67]
X3LYP (Exchange and correlation) [68]
CAMB3LYP (Range separated exchange and LYP correlation) [69]
Becke97 (Exchange and correlation within a hybrid scheme) [70, 59]
Becke971 (Exchange and correlation within a hybrid scheme) [71, 59]
Becke972 (Exchange and correlation within a hybrid scheme) [72, 59]
B97D (Exchange and correlation and empirical dispersion correction) [73]
HCTH (Exchange correlation within a hybrid scheme) [71, 59]
HCTH120 (Exchange correlation within a hybrid scheme) [74, 59]
HCTH147 (Exchange correlation within a hybrid scheme) [74, 59]
HCTH407 (Exchange correlation within a hybrid scheme) [75, 59]
The B97X functionals developed by Chai and HeadGordon [76] (Exchange and correlation within a hybrid scheme, with longrange correction, see further in this manual for details)
The B97XD3 functional (ExchangeCorrelation within a hybrid scheme, with longrange correction and dispersion correction, see futher in this manual for details) [77]
BOP (Becke88 exchange plus the “oneparameter progressive” correlation functional, OP) [80]
PBEOP (PBE Exchange plus the OP correlation functional) [80]
SOGGA (Exchange plus the PBE correlation functional) [81]
SOGGA11 (Exchange and Correlation) [82]
SOGGA11X (Exchange and Correlation within a hybrid scheme, with reoptimized SOGGA11 parameters) [83]
LRCPBEPBE (Longrange corrected PBE exchange and PBE correlation) [84]
LRCPBEhPBE (Longrange corrected hybrid PBE exchange and PBE correlation) [85]
Note: The OP correlation functional used in BOP has been parameterized for use with Becke88 exchange, whereas in the PBEOP functional, the same correlation ansatz is reparameterized for use with PBE exchange. These two versions of OP correlation are available as the correlation functionals (B88)OP and (PBE)OP. The BOP functional, for example, consists of (B88)OP correlation combined with Becke88 exchange.
MetaGGA functionals involving the kinetic energy density (), and or the Laplacian of the electron density:
VSXC (Exchange and Correlation) [29]
TPSS (Exchange and Correlation in a single nonempirical scheme) [33, 34]
TPSSh (Exchange and Correlation within a nonempirical hybrid scheme) [86]
BMK (Exchange and Correlation within a hybrid scheme) [87]
M05 (Exchange and Correlation within a hybrid scheme) [88, 89]
M052X (Exchange and Correlation within a hybrid scheme) [90, 89]
M06HF (Exchange and Correlation within a hybrid scheme) [92, 89]
M06 (Exchange and Correlation within a hybrid scheme) [93, 89]
M062X (Exchange and Correlation within a hybrid scheme) [93, 89]
M08HX (Exchange and Correlation within a hybrid scheme) [94]
M08SO (Exchange and Correlation within a hybrid scheme) [94]
M11L (Exchange and Correlation) [95]
M11 (Exchange and Correlation within a hybrid scheme, with longrange correction) [96]
B95 (Correlation) [97]
B1B95 (Exchange and Correlation) [97]
PK06 (Correlation) [98]
The M05D and M06D3 functionals (ExchangeCorrelation within a hybrid scheme, with longrange correction and dispersion correction, see futher in this manual for details) [99, 77]
HyperGGA functionals:
B05 (A full exactexchange KohnSham scheme of Becke that accounts for static correlation via realspace corrections) [35, 39, 40]
mB05 (Modified B05 method that has simpler functional form and SCF potential) [100]
PSTS (HyperGGA functional of PerdewStaroverovTaoScuseria) [38]
MCY2 (The adiabatic connectionbased MCY2 functional) [36, 37, 34]
Fifthrung, doublehybrid (DH) functionals:
B97X2 (Exchange and Correlation within a DH generalization of the LC corrected B97X scheme) [42]
B2PLYP (another DH scheme proposed by Grimme, based on GGA exchange and correlation functionals) [73]
XYG3 and XYGJOS (an efficient DH scheme based on generalization of B3LYP) [101]
PBE0DH [102] and PBE02 [103] functionals (nonempirical DH scheme based on the PBE functional).
In addition to the above functional types, QChem contains the Empirical Density Functional 1 (EDF1), developed by Adamson, Gill and Pople [104]. EDF1 is a combined exchange and correlation functional that is specifically adapted to yield good results with the relatively modestsized 631+G* basis set, by direct fitting to thermochemical data. It has the interesting feature that exact exchange mixing was not found to be helpful with a basis set of this size. Furthermore, for this basis set, the performance substantially exceeded the popular B3LYP functional, while the cost of the calculations is considerably lower because there is no need to evaluate exact (nonlocal) exchange. We recommend consideration of EDF1 instead of either B3LYP or BLYP for density functional calculations on large molecules, for which basis sets larger than 631+G* may be too computationally demanding.
EDF2, another Empirical Density Functional, was developed by Ching Yeh Lin and Peter Gill [105] in a similar vein to EDF1, but is specially designed for harmonic frequency calculations. It was optimized using the ccpVTZ basis set by fitting into experimental harmonic frequencies and is designed to describe the potential energy curvature well. Fortuitously, it also performs better than B3LYP for thermochemical properties.
A few more words deserve the hybrid functionals [65], where several different exchange and correlation functionals can be combined linearly to form a hybrid functional. These have proven successful in a number of reported applications. However, since the hybrid functionals contain HF exchange they are more expensive that pure DFT functionals. QChem has incorporated two of the most popular hybrid functionals, B3LYP [106] and B3PW91 [30], with the additional option for users to define their own hybrid functionals via the $xc_functional keyword (see userdefined functionals in Section 4.3.19, below). Among the latter, a recent new hybrid combination available in QChem is the ’B3tLap’ functional, based on Becke’s B88 GGA exchange and the ’tLap’ (or ’PK06’) metaGGA correlation [98, 107]. This hybrid combination is on average more accurate than B3LYP, BMK, and M06 functionals for thermochemistry and better than B3LYP for reaction barriers, while involving only five fitting parameters. Another hybrid functional in QChem that deserves attention is the hybrid extension of the BR89B94 metaGGA functional [31, 107]. This hybrid functional yields a very good thermochemistry results, yet has only three fitting parameters.
In addition, QChem now includes the M05 and M06 suites of density functionals. These are designed to be used only with certain definite percentages of HartreeFock exchange. In particular, M06L [91] is designed to be used with no HartreeFock exchange (which reduces the cost for large molecules), and M05 [88], M052X [90], M06, and M062X [93] are designed to be used with 28%, 56%, 27%, and 54% HartreeFock exchange. M06HF [92] is designed to be used with 100% HartreeFock exchange, but it still contains some local DFT exchange because the 100% nonlocal HartreeFock exchange replaces only some of the local exchange.
Note: The hybrid functionals are not simply a pairing of an exchange and correlation functional, but are a combined exchangecorrelation functional (i.e., BLYP and B3LYP vary in the correlation contribution in addition to the exchange part).
As pointed out in Ref. Dreuw:2003 and elsewhere, the description of chargetransfer excited states within density functional theory (or more precisely, timedependent DFT, which is discussed in Section 6.3) requires full (100%) nonlocal HartreeFock exchange, at least in the limit of large donor–acceptor distance. Hybrid functionals such as B3LYP [106] and PBE0 [64] that are wellestablished and in widespread use, however, employ only 20% and 25% HartreeFock exchange, respectively. While these functionals provide excellent results for many groundstate properties, they cannot correctly describe the distance dependence of chargetransfer excitation energies, which are enormously underestimated by most common density functionals. This is a serious problem in any case, but it is a catastrophic problem in large molecules and in clusters, where TDDFT often predicts a nearcontinuum of of spurious, lowlying charge transfer states [109, 110]. The problems with TDDFT’s description of charge transfer are not limited to large donor–acceptor distances, but have been observed at 2 separation, in systems as small as uracil– [109]. Rydberg excitation energies also tend to be substantially underestimated by standard TDDFT.
One possible avenue by which to correct such problems is to parameterize functionals that contain 100% HartreeFock exchange. To date, few such functionals exist, and those that do (such as M06HF) contain a very large number of empirical adjustable parameters. An alternative option is to attempt to preserve the form of common GGAs and hybrid functionals at short range (i.e., keep the 25% HartreeFock exchange in PBE0) while incorporating 100% HartreeFock exchange at long range. Functionals along these lines are known variously as “Coulombattenuated” functionals, “rangeseparated” functionals, or (our preferred designation) “longrangecorrected” (LRC) density functionals. Whatever the nomenclature, these functionals are all based upon a partition of the electronelectron Coulomb potential into long and shortrange components, using the error function (erf):
(4.45) 
The first term on the right in Eq. () is singular but shortrange, and decays to zero on a length scale of , while the second term constitutes a nonsingular, longrange background. The basic idea of LRCDFT is to utilize the shortrange component of the Coulomb operator in conjunction with standard DFT exchange (including any component of HartreeFock exchange, if the functional is a hybrid), while at the same time incorporating full HartreeFock exchange using the longrange part of the Coulomb operator. This provides a rigorously correct description of the longrange distance dependence of chargetransfer excitation energies, but aims to avoid contaminating shortrange exchangecorrelation effects with extra HartreeFock exchange.
Consider an exchangecorrelation functional of the form
(4.46) 
in which is the correlation energy, is the (local) GGA exchange energy, and is the (nonlocal) HartreeFock exchange energy. The constant denotes the fraction of HartreeFock exchange in the functional, therefore for GGAs, for B3LYP, for PBE0, etc.. The LRC version of the generic functional in Eq. () is
(4.47) 
in which the designations “SR” and “LR” in the various exchange energies indicate that these components of the functional are evaluated using either the shortrange (SR) or the longrange (LR) component of the Coulomb operator. (The correlation energy is evaluated using the full Coulomb operator.) The LRC functional in Eq. () incorporates full HartreeFock exchange in the asymptotic limit via the final term, . To fully specify the LRC functional, one must choose a value for the range separation parameter in Eq. (); in the limit , the LRC functional in Eq. () reduces to the original functional in Eq. (), while the limit corresponds to a new functional, . It is well known that full HartreeFock exchange is inappropriate for use with most contemporary GGA correlation functionals, so the latter limit is expected to perform quite poorly. Values of are probably not worth considering [111, 84].
Evaluation of the short and longrange HartreeFock exchange energies is straightforward [112], so the crux of LRCDFT rests upon the form of the shortrange GGA exchange energy. Several different shortrange GGA exchange functionals are available in QChem, including shortrange variants of B88 and PBE exchange described by Hirao and coworkers [113, 114], an alternative formulation of shortrange PBE exchange proposed by Scuseria and coworkers [115], and several shortrange variants of B97 introduced by Chai and HeadGordon [76, 116, 117, 42]. The reader is referred to these papers for additional methodological details.
These LRCDFT functionals have been shown to remove the nearcontinuum of spurious chargetransfer excited states that appear in largescale TDDFT calculations [111]. However, certain results depend sensitively upon the rangeseparation parameter [110, 111, 84, 85], and the results of LRCDFT calculations must therefore be interpreted with caution, and probably for a range of values. In two recent benchmark studies of several LRC density functionals, Rohrdanz and Herbert [84, 85] have considered the errors engendered, as a function of , in both groundstate properties and also TDDFT vertical excitation energies. In Ref. Lange:2008, the sensitivity of valence excitations versus chargetransfer excitation energies in TDDFT was considered, again as a function of . A careful reading of these references is suggested prior to performing any LRCDFT calculations.
Within QChem 3.2, there are three ways to perform LRCDFT calculations.
The form of is different for each different GGA exchange functional, and shortrange versions of B88 and PBE exchange are available in QChem through the efforts of the Herbert group. Versions of B88 and PBE, in which the Coulomb attenuation is performed according to the procedure of Hirao [114], are denoted as B88 and PBE, respectively (since , rather than , is the Hirao group’s notation for the rangeseparation parameter). Alternatively, a shortrange version of PBE exchange called PBE is available, which is constructed according to the prescription of Scuseria and coworkers [115].
These shortrange exchange functionals can be used in the absence of longrange HartreeFock exchange, and using a combination of PBE exchange and PBE correlation, a user could, for example, employ the shortrange hybrid functional recently described by Heyd, Scuseria, and Ernzerhof [118]. Shortrange hybrids appear to be most appropriate for extended systems, however. Thus, within QChem, shortrange GGAs should be used with longrange HartreeFock exchange, as in Eq. . Longrange HartreeFock exchange is requested by setting LRC_DFT to TRUE.
LRCDFT is thus available for any functional whose exchange component consists of some combination of HartreeFock, B88, and PBE exchange (e.g., BLYP, PBE, PBE0, BOP, PBEOP, and various userspecified combinations, but not B3LYP or other functionals whose exchange components are more involved). Having specified such a functional via the EXCHANGE and CORRELATION variables, a user may request the corresponding LRC functional by setting LRC_DFT to TRUE. Longrangecorrected variants of PBE0, BOP, and PBEOP must be obtained through the appropriate userspecified combination of exchange and correlation functionals (as demonstrated in the example below). In any case, the value of must also be specified by the user. Analytic energy gradients are available but analytic Hessians are not. TDDFT vertical excitation energies are also available.
LRC_DFT
Controls the application of longrangecorrected DFT
TYPE:
LOGICAL
DEFAULT:
FALSE
OPTIONS:
FALSE
(or 0) Do not apply longrange correction.
TRUE
(or 1) Use the longrangecorrected version of the requested functional.
RECOMMENDATION:
Longrange correction is available for any combination of HartreeFock, B88, and PBE exchange (along with any standalone correlation functional).
OMEGA
Sets the Coulomb attenuation parameter .
TYPE:
INTEGER
DEFAULT:
No default
OPTIONS:
Corresponding to , in units of bohr
RECOMMENDATION:
None
Example 4.23 Application of LRCBOP to .
$comment To obtain LRCBOP, a shortrange version of BOP must be specified, using muB88 shortrange exchange plus (B88)OP correlation, which is the version of OP parameterized for use with B88. $end $molecule 1 2 O 1.347338 .017773 .071860 H 1.824285 .813088 .117645 H 1.805176 .695567 .461913 O 1.523051 .002159 .090765 H .544777 .024370 .165445 H 1.682218 .174228 .849364 $end $rem EXCHANGE GEN BASIS 631(1+,3+)G* LRC_DFT TRUE OMEGA 330 ! = 0.330 a.u. $end $xc_functional C (B88)OP 1.0 X muB88 1.0 $end
Regarding the choice of functionals and values, it has been found that the Hirao and Scuseria ansatz afford virtually identical TDDFT excitation energies, for all values of [85]. Thus, functionals based on PBE versus PBE provide the same excitation energies, as a function of . However, the PBE functional appears to be somewhat superior in the sense that it can provide accurate TDDFT excitation energies and accurate groundstate properties using the same value of [85], whereas this does not seem to be the case for functionals based on B88 or PBE [84].
Recently, Rohrdanz et al. [85] have published a thorough benchmark study of both ground and excitedstate properties, using the “LRCPBEh” functional, a hybrid (hence the “h”) that contains a fraction of shortrange HartreeFock exchange in addition to full longrange HartreeFock exchange:
(4.48) 
The statisticallyoptimal parameter set, consider both groundstate properties and TDDFT excitation energies together, was found to be and bohr [85]. With these parameters, the LRCPBEh functional outperforms the traditional hybrid functional PBE0 for groundstate atomization energies and barrier heights. For TDDFT excitation energies corresponding to localized excitations, TDPBE0 and TDLRCPBEh show similar statistical errors of 0.3 eV, but the latter functional also exhibits only 0.3 eV errors for chargetransfer excitation energies, whereas the statistical error for TDPBE0 chargetransfer excitation energies is 3.0 eV! Caution is definitely warranted in the case of chargetransfer excited states, however, as these excitation energies are very sensitive to the precise value of [110, 85]. It was later found that the parameter set (, bohr) provides similar (statistical) performance to that described above, although the predictions for specific chargetransfer excited states can be somewhat different as compared to the original parameter set [110].
Example 4.24 Application of LRCPBEh to the — heterodimer at 5 separation.
$comment This example uses the "optimal" parameter set discussed above. It can also be run by setting "EXCHANGE LRCWPBEhPBE". $end $molecule 0 1 C 0.670604 0.000000 0.000000 C 0.670604 0.000000 0.000000 H 1.249222 0.929447 0.000000 H 1.249222 0.929447 0.000000 H 1.249222 0.929447 0.000000 H 1.249222 0.929447 0.000000 C 0.669726 0.000000 5.000000 C 0.669726 0.000000 5.000000 F 1.401152 1.122634 5.000000 F 1.401152 1.122634 5.000000 F 1.401152 1.122634 5.000000 F 1.401152 1.122634 5.000000 $end $rem EXCHANGE GEN BASIS 631+G* LRC_DFT TRUE OMEGA 200 ! = 0.2 a.u. CIS_N_ROOTS 4 CIS_TRIPLETS FALSE $end $xc_functional C PBE 1.00 X wPBE 0.80 X HF 0.20 $end
Both LRC functionals and asymptotic corrections (Section 4.3.10.1) are thought to reduce selfinteraction error in approximate density functional theory. A convenient way to quantify (or at least depict) this error is by plotting the DFT energy as a function of the (fractional) number of electrons, , as the exchangecorrelation potential changes discontinuously as passes through an integer, and thus a plot of versus abruptly changes slope at integer values of [119]. Examination of such plots has been suggested as a means to “tune” the fraction of shortrange exchange in an LRC function [120], while the range separation parameter can be tuned so as to achieve the condition (exact for the HohenbergKohn functional) , where IE denotes the molecule’s lowest ionization energy [121].
Example 4.25 Example of a DFT job with a fraction number of electrons. Here, we make the anion of fluoride by subracting a fraction of an electron from the HOMO of F.
$comment Subtracting a whole electron recovers the energy of F. Adding electrons to the LUMO is possible as well. $end $rem exchangeh b3lyp basis 631+G* fractional_electron 500 ! /divide by 1000 to get the fraction, 0.5 here. $end $molecule 2 2 F $end
The BaerNeuhauserLivshits (BNL) functional [78, 79] is also based on the range separation of the Coulomb operator in Eq. . Its functional form resembles Eq. :
(4.49) 
where the recommended GGA correlation functional is LYP. The recommended GGA exchange functional is BNL, which is described by a local functional [122]. For ground state properties, the optimized value for (scaling factor for the BNL exchange functional) was found to be 0.9.
The value of in BNL calculations can be chosen in several different ways. For example, one can use the optimized value =0.5 bohr. For calculation of excited states and properties related to orbital energies, it is strongly recommend to tune as described below[121, 123].
Systemspecific optimization of is based on Koopmans conditions that would be satisfied for the exact functional[121], that is, is varied until the Koopmans IE/EA for the HOMO/LUMO is equal to IE/EA. Based on published benchmarks [79, 124], this systemspecific approach yields the most accurate values of IEs and excitation energies.
The script that optimizes is called OptOmegaIPEA.pl and is located in the $QC/bin directory. The script optimizes in the range 0.10.8 (100800). See the script for the instructions how to modify the script to optimize in a broader range. To execute the script, you need to create three inputs for a BNL job using the same geometry and basis set for a neutral molecule (N.in), anion (M.in), and cation (P.in), and then type 'OptOmegaIPEA.pl >& optomega'. The script will run creating outputs for each step (N_*, P_*, M_*) writing the optimization output into optomega.
A similar script, OptOmegaIP.pl, will optimize to satisfy the Koopmans condition for the IP only. This script minimizes , not the absolute values.
Note: (i) If the system does not have positive EA, then the tuning should be done according to the IP condition only. The IPEA script will yield a wrong value of in such cases.
(ii) In order for the scripts to work, one must specify SCF_FINAL_PRINT=1 in the inputs. The scripts look for specific regular expressions and will not work correctly without this keyword.
(iii) When tuning omega we recommend taking the amount of X BNL in the XC part as 1.0 and not 0.9.
The $xc_functional keyword for a BNL calculation reads:
$xc_functional X HF 1.0 X BNL 0.9 C LYP 1.0 $end
and the $rem keyword reads
$rem EXCHANGE GENERAL SEPARATE_JK TRUE OMEGA 500 != 0.5 Bohr$^{1}$ DERSCREEN FALSE !if performing unrestricted calcn IDERIV 0 !if performing unrestricted Hessian evaluation $end
Also available in QChem are the B97 [76], B97X [76], B97XD [116], and B97X2 [42] functionals, recently developed by Chai and HeadGordon. These authors have proposed a very simple ansatz to extend any E to E, as long as the SR operator has considerable spatial extent [76, 117]. With the use of flexible GGAs, such as Becke97 functional [70], their new LRC hybrid functionals [76, 116, 117] outperform the corresponding global hybrid functionals (i.e., B97) and popular hybrid functionals (e.g., B3LYP) in thermochemistry, kinetics, and noncovalent interactions, which has not been easily achieved by the previous LRC hybrid functionals. In addition, the qualitative failures i of the commonly used hybrid density functionals in some “difficult problems”, such as dissociation of symmetric radical cations and longrange chargetransfer excitations, are significantly reduced by these new functionals [76, 116, 117]. Analytical gradients and analytical Hessians are available for B97, B97X, and B97XD.
Example 4.26 Application of B97 functional to nitrogen dimer.
$comment Geometry optimization, followed by a TDDFT calculation. $end $molecule 0 1 N1 N2 N1 1.1 $end $rem jobtype opt exchange omegaB97 basis 631G* $end @@@ $molecule READ $end $rem jobtype sp exchange omegaB97 basis 631G* scf_guess READ cis_n_roots 10 rpa true $end
Example 4.27 Application of B97X functional to nitrogen dimer.
$comment Frequency calculation (with analytical Hessian methods). $end $molecule 0 1 N1 N2 N1 1.1 $end $rem jobtype freq exchange omegaB97X basis 631G* $end
Among these new LRC hybrid functionals, B97XD is a DFTD (density functional theory with empirical dispersion corrections) functional, where the total energy is computed as the sum of a DFT part and an empirical atomicpairwise dispersion correction. Relative to B97 and B97X, B97XD is significantly superior for nonbonded interactions, and very similar in performance for bonded interactions. However, it should be noted that the remaining shortrange selfinteraction error is somewhat larger for B97XD than for B97X than for B97. A careful reading of Refs. Chai:2008a,Chai:2008b,Chai:2008c is suggested prior to performing any DFT and TDDFT calculations based on variations of B97 functional. B97XD functional automatically involves two keywords for the dispersion correction, DFT_D and DFT_D_A, which are described in Section 4.3.6.
Example 4.28 Application of B97XD functional to methane dimer.
$comment Geometry optimization. $end $molecule 0 1 C 0.000000 0.000323 1.755803 H 0.887097 0.510784 1.390695 H 0.887097 0.510784 1.390695 H 0.000000 1.024959 1.393014 H 0.000000 0.001084 2.842908 C 0.000000 0.000323 1.755803 H 0.000000 0.001084 2.842908 H 0.887097 0.510784 1.390695 H 0.887097 0.510784 1.390695 H 0.000000 1.024959 1.393014 $end $rem jobtype opt exchange omegaB97XD basis 631G* $end
Similar to the existing doublehybrid density functional theory (DHDFT) [41, 125, 126, 127, 101], which is described in Section 4.3.9, LRCDFT can be extended to include nonlocal orbital correlation energy from secondorder MøllerPlesset perturbation theory (MP2) [128], that includes a samespin (ss) component , and an oppositespin (os) component of PT2 correlation energy. The two scaling parameters, and , are introduced to avoid doublecounting correlation with the LRC hybrid functional.
(4.50) 
Among the B97 series, B97X2 [42] is a longrange corrected doublehybrid (DH) functional, which can greatly reduce the selfinteraction errors (due to its high fraction of HartreeFock exchange), and has been shown significantly superior for systems with bonded and nonbonded interactions. Due to the sensitivity of PT2 correlation energy with respect to the choices of basis sets, B97X2 was parameterized with two different basis sets. B97X2(LP) was parameterized with the 6311++G(3df,3pd) basis set (the large Pople type basis set), while B97X2(TQZ) was parameterized with the TQ extrapolation to the basis set limit. A careful reading of Ref. Chai:2009 is thus highly advised.
B97X2(LP) and B97X2(TQZ) automatically involve three keywords for the PT2 correlation energy, DH, SSS_FACTOR and SOS_FACTOR, which are described in Section 4.3.9. The PT2 correlation energy can also be computed with the efficient resolutionofidentity (RI) methods (see Section 5.5).
Example 4.29 Application of B97X2(LP) functional to LiH molecules.
$comment Geometry optimization and frequency calculation on LiH, followed by singlepoint calculations with nonRI and RI approaches. $end $molecule 0 1 H Li H 1.6 $end $rem jobtype opt exchange omegaB97X2(LP) correlation mp2 basis 6311++G(3df,3pd) $end @@@ $molecule READ $end $rem jobtype freq exchange omegaB97X2(LP) correlation mp2 basis 6311++G(3df,3pd) $end @@@ $molecule READ $end $rem jobtype sp exchange omegaB97X2(LP) correlation mp2 basis 6311++G(3df,3pd) $end @@@ $molecule READ $end $rem jobtype sp exchange omegaB97X2(LP) correlation rimp2 basis 6311++G(3df,3pd) aux_basis rimp2augccpvtz $end
Example 4.30 Application of B97X2(TQZ) functional to LiH molecules.
$comment Singlepoint calculations on LiH. $end $molecule 0 1 H Li H 1.6 $end $rem jobtype sp exchange omegaB97X2(TQZ) correlation mp2 basis ccpvqz $end @@@ $molecule READ $end $rem jobtype sp exchange omegaB97X2(TQZ) correlation rimp2 basis ccpvqz aux_basis rimp2ccpvqz $end
M05D [99], M06D3 and B97XD3 functionals [77], developed by the Chai group, improve on the B97XD functional mentioned above. Similar to the B97XD functional, these functionals also involves emperical atomicpairwise dispersion corrections, and automatically involves the keywords for dispersion correction.
Analytical gradient is available for all three functionals in QChem, and analytical Hessian is available for B97XD3.
Example 4.31 Applications of M05D to a methane dimer.
$comment Geometry optimization, followed by singlepoint TDDFT calculations using a larger basis set. $end $molecule 0 1 C 0.000000 0.000323 1.755803 H 0.887097 0.510784 1.390695 H 0.887097 0.510784 1.390695 H 0.000000 1.024959 1.393014 H 0.000000 0.001084 2.842908 C 0.000000 0.000323 1.755803 H 0.000000 0.001084 2.842908 H 0.887097 0.510784 1.390695 H 0.887097 0.510784 1.390695 H 0.000000 1.024959 1.393014 $end $rem jobtype opt method wM05D basis 631G* $end @@@ $molecule READ $end $rem jobtype sp method wM05D basis 6311++G** cis_n_roots 30 rpa true $end
Example 4.32 Applications of M06D3 to a methane dimer.
$comment Geometry optimization, followed by singlepoint TDDFT calculations using a larger basis set. $end $molecule 0 1 C 0.000000 0.000323 1.755803 H 0.887097 0.510784 1.390695 H 0.887097 0.510784 1.390695 H 0.000000 1.024959 1.393014 H 0.000000 0.001084 2.842908 C 0.000000 0.000323 1.755803 H 0.000000 0.001084 2.842908 H 0.887097 0.510784 1.390695 H 0.887097 0.510784 1.390695 H 0.000000 1.024959 1.393014 $end $rem jobtype opt method wM06D3 basis 631G* $end @@@ $molecule READ $end $rem jobtype sp method wM06D3 basis 6311++G** cis_n_roots 30 rpa true $end
Example 4.33 Applications of B97XD3 to a methane dimer.
$comment Geometry optimization, followed by singlepoint TDDFT calculations using a larger basis set (with analytical Hessian). $end $molecule 0 1 C 0.000000 0.000323 1.755803 H 0.887097 0.510784 1.390695 H 0.887097 0.510784 1.390695 H 0.000000 1.024959 1.393014 H 0.000000 0.001084 2.842908 C 0.000000 0.000323 1.755803 H 0.000000 0.001084 2.842908 H 0.887097 0.510784 1.390695 H 0.887097 0.510784 1.390695 H 0.000000 1.024959 1.393014 $end $rem jobtype opt method wB97XD3 basis 631G* $end @@@ $molecule READ $end $rem jobtype sp method wB97XD3 basis 6311++G** cis_n_roots 30 rpa true $end
The Minnesota family of functional by Truhlar’s group has been recently updated by adding two new functionals: M11L [95] and M11 [96]. The M11 functional is a longrange corrected metaGGA, obtained by using the LRC scheme of Chai and HeadGordon (see above), with the successful parameterization of the Minnesota metaGGA functionals:
(4.51) 
with the percentage of HartreeFock exchange at short range X being 42.8. An extension of the LRC scheme to local functional (no HF exchange) was introduced in the M11L functional by means of the dualrange exchange:
(4.52) 
The correct longrange scheme is selected automatically with the input keywords. A careful reading of the references [95, 96] is suggested prior to performing any calculations with the M11 functionals.
Example 4.34 Application of M11 functional to water molecule
$comment Optimization of H2O with M11 $end $molecule 0 1 O 0.000000 0.000000 0.000000 H 0.000000 0.000000 0.956914 H 0.926363 0.000000 0.239868 $end $rem jobtype opt exchange m11 basis 631+G(d,p) $end
QChem includes four nonlocal correlation functionals that describe longrange dispersion (i.e. van der Waals) interactions:
vdWDF04, developed by Langreth, Lundqvist, and coworkers [129, 130] and implemented as described in Ref. [131];
vdWDF10 (also known as vdWDF2), which is a reparameterization [132] of vdWDF04, implemented in the same way as its precursor [131];
VV09, developed [133] and implemented [134] by Vydrov and Van Voorhis;
VV10 by Vydrov and Van Voorhis [135].
All these functionals are implemented selfconsistently and analytic gradients with respect to nuclear displacements are available [131, 134, 135]. The nonlocal correlation is governed by the $rem variable NL_CORRELATION, which can be set to one of the four values: vdWDF04, vdWDF10, VV09, or VV10. Note that vdWDF04, vdWDF10, and VV09 functionals are used in combination with LSDA correlation, which must be specified explicitly. For instance, vdWDF10 is invoked by the following keyword combination:
CORRELATION PW92 NL_CORRELATION vdWDF10
VV10 is used in combination with PBE correlation, which must be added explicitly. In addition, the values of two parameters, and must be specified for VV10. These parameters are controlled by the $rem variables NL_VV_C and NL_VV_B, respectively. For instance, to invoke VV10 with and , the following input is used:
CORRELATION PBE NL_CORRELATION VV10 NL_VV_C 93 NL_VV_B 590
The variable NL_VV_C may also be specified for VV09, where it has the same meaning. By default, is used in VV09 (i.e. NL_VV_C is set to 89). However, in VV10 neither nor are assigned a default value and must always be provided in the input.
As opposed to local (LSDA) and semilocal (GGA and metaGGA) functionals, evaluated as a single 3D integral over space [see Eq. ()], nonlocal functionals require double integration over the spatial variables:
(4.53) 
In practice, this double integration is performed numerically on a quadrature grid [131, 134, 135]. By default, the SG1 quadrature (described in Section 4.3.15 below) is used to evaluate , but a different grid can be requested via the $rem variable NL_GRID. The nonlocal energy is rather insensitive to the fineness of the grid, so that SG1 or even SG0 grids can be used in most cases. However, a finer grid may be required for the (semi)local parts of the functional, as controlled by the XC_GRID variable.
Example 4.35 Geometry optimization of the methane dimer using VV10 with rPW86 exchange.
$molecule 0 1 C 0.000000 0.000140 1.859161 H 0.888551 0.513060 1.494685 H 0.888551 0.513060 1.494685 H 0.000000 1.026339 1.494868 H 0.000000 0.000089 2.948284 C 0.000000 0.000140 1.859161 H 0.000000 0.000089 2.948284 H 0.888551 0.513060 1.494685 H 0.888551 0.513060 1.494685 H 0.000000 1.026339 1.494868 $end $rem JobType Opt BASIS augccpVTZ EXCHANGE rPW86 CORRELATION PBE XC_GRID 75000302 NL_CORRELATION VV10 NL_GRID 1 NL_VV_C 93 NL_VV_B 590 $end
In the above example, an EML(75,302) grid is used to evaluate the rPW86 exchange and PBE correlation, but a coarser SG1 grid is used for the nonlocal part of VV10.
Thanks to the efforts of the Sherrill group, the popular empirical dispersion corrections due to Grimme [73] are now available in QChem. Energies, analytic gradients, and analytic second derivatives are available. Grimme’s empirical dispersion corrections can be added to any of the density functionals available in QChem.
DFTD methods add an extra term,
(4.54)  
(4.55)  
(4.56) 
where is a global scaling parameter (near unity), is a damping parameter meant to help avoid doublecounting correlation effects at short range, is a global scaling parameter for the damping function, and is the sum of the van der Waals radii of atoms A and B.
DFTD using Grimme’s parameters may be turned on using
DFT_D EMPIRICAL_GRIMME
Grimme has suggested scaling factors of 0.75 for PBE, 1.2 for BLYP, 1.05 for BP86, and 1.05 for B3LYP; these are the default values of when those functionals are used. Otherwise, the default value of is 1.0.
It is possible to specify different values of , , the atomic coefficients, or the van der Waals radii by using the $empirical_dispersion keyword; for example:
$empirical_dispersion S6 1.1 D 10.0 C6 Ar 4.60 Ne 0.60 VDW_RADII Ar 1.60 Ne 1.20 $end
Any values not specified explicitly will default to the values in Grimme’s model.
The empirical dispersion correction from Chai and HeadGordon [116] employs a different damping function and can be activated by using
DFT_D EMPIRICAL_CHG
It uses another keyword DFT_D_A to control the strength of dispersion corrections.
DFT_D
Controls the application of DFTD or DFTD3 scheme.
TYPE:
LOGICAL
DEFAULT:
None
OPTIONS:
FALSE
(or 0) Do not apply the DFTD or DFTD3 scheme
EMPIRICAL_GRIMME
dispersion correction from Grimme
EMPIRICAL_CHG
dispersion correction from Chai and HeadGordon
EMPIRICAL_GRIMME3
dispersion correction from Grimme’s DFTD3 method
(see Section 4.3.8)
RECOMMENDATION:
NONE
DFT_D_A
Controls the strength of dispersion corrections in the ChaiHeadGordon DFTD scheme in Eq.(3) of Ref. Chai:2008b.
TYPE:
INTEGER
DEFAULT:
600
OPTIONS:
n
Corresponding to .
RECOMMENDATION:
Use default, i.e.,
While standard DFT functionals describe chemical bonds relatively well, one major deficiency is their inability to cope with dispersion interactions, i.e., van der Waals (vdW) interactions. Becke and Johnson have proposed a conceptually simple yet accurate dispersion model called the exchangedipole model (XDM) [35, 136]. In this model the dispersion attraction emerges from the interaction between the instant dipole moment of the exchange hole in one molecule and the induced dipole moment in another. It is a conceptually simple but powerful approach that has been shown to yield very accurate dispersion coefficients without fitting parameters. This allows the calculation of both intermolecular and intramolecular dispersion interactions within a single DFT framework. The implementation and validation of this method in the QChem code is described in Ref. Kong:2009.
Fundamental to the XDM model is the calculation of the norm of the dipole moment of the exchange hole at a given point:
(4.57) 
where labels the spin and is the exchangehole function. The XDM version that is implemented in QChem employs the BeckeRoussel (BR) model exchangehole function. It was not given in an analytical form and one had to determine its value at each grid point numerically. QChem has developed for the first time an analytical expression for this function based on nonlinear interpolation and spline techniques, which greatly improves efficiency as well as the numerical stability [30].
There are two different damping functions used in the XDM model of Becke and Johnson. One of them uses only the intermolecular C6 dispersion coefficient. In its QChem implementation it is denoted as "XDM6". In this version the dispersion energy is computed as
(4.58) 
where is a universal parameter, is the distance between atoms and , and is the sum of the absolute values of the correlation energy of free atoms and . The dispersion coefficients is computed as
(4.59) 
where is the exchange hole dipole moment of the atom, and is the effective polarizability of the atom in the molecule.
The XDM6 scheme is further generalized to include higherorder dispersion coefficients, which leads to the “XDM10” model in QChem implementation. The dispersion energy damping function used in XDM10 is
(4.60) 
where , and are dispersion coefficients computed at higherorder multipole (including dipole, quadrupole and octopole) moments of the exchange hole [138]. In above, is the sum of the effective vdW radii of atoms and , which is a linear function of the so called critical distance between atoms and :
(4.61) 
The critical distance, , is computed by averaging these three distances:
(4.62) 
In the XDM10 scheme there are two universal parameters, and . Their default values of 0.83 and 1.35, respectively, are due to Johnson and Becke [136], determined by least square fitting to the binding energies of a set of intermolecular complexes. Please keep in mind that these values are not the only possible optimal set to use with XDM. Becke’s group has suggested later on several other XC functional combinations with XDM that employ different and values. The user is advised to consult their recent papers for more details (e.g., Refs. Becke:2010,Kannemann:2010).
The computed vdW energy is added as a postSCF correction. In addition, QChem also has implemented the first and second nuclear derivatives of vdW energy correction in both the XDM6 and XDM10 schemes.
Listed below are a number of useful options to customize the vdW calculation based on the XDM DFT approach.
DFTVDW_JOBNUMBER
Basic vdW job control
TYPE:
INTEGER
DEFAULT:
0
OPTIONS:
0
Do not apply the XDM scheme.
1
add vdW as energy/gradient correction to SCF.
2
add VDW as a DFT functional and do full SCF (this option only works with C6 XDM formula).
RECOMMENDATION:
none
DFTVDW_METHOD
Choose the damping function used in XDM
TYPE:
INTEGER
DEFAULT:
1
OPTIONS:
1
use Becke’s damping function including C6 term only.
2
use Becke’s damping function with higherorder (C8,C10) terms.
RECOMMENDATION:
none
DFTVDW_MOL1NATOMS
The number of atoms in the first monomer in dimer calculation
TYPE:
INTEGER
DEFAULT:
0
OPTIONS:
0NATOMS
default 0
RECOMMENDATION:
none
DFTVDW_KAI
Damping factor K for C6 only damping function
TYPE:
INTEGER
DEFAULT:
800
OPTIONS:
101000
default 800
RECOMMENDATION:
none
DFTVDW_ALPHA1
Parameter in XDM calculation with higherorder terms
TYPE:
INTEGER
DEFAULT:
83
OPTIONS:
101000
RECOMMENDATION:
none
DFTVDW_ALPHA2
Parameter in XDM calculation with higherorder terms.
TYPE:
INTEGER
DEFAULT:
155
OPTIONS:
101000
RECOMMENDATION:
none
DFTVDW_USE_ELE_DRV
Specify whether to add the gradient correction to the XDM energy. only valid with Becke’s C6 damping function using the interpolated BR89 model.
TYPE:
LOGICAL
DEFAULT:
1
OPTIONS:
1
use density correction when applicable (default).
0
do not use this correction (for debugging purpose)
RECOMMENDATION:
none
DFTVDW_PRINT
Printing control for VDW code
TYPE:
INTEGER
DEFAULT:
1
OPTIONS:
0 no printing.
1
minimum printing (default)
2
debug printing
RECOMMENDATION:
none
Example 4.36 Below is a sample input illustrating a frequency calculation of a vdW complex consisted of He atom and N2 molecule.
$molecule 0 1 He .0 .0 3.8 N .000000 .000000 0.546986 N .000000 .000000 0.546986 $end $rem JOBTYPE FREQ IDERIV 2 EXCHANGE B3LYP !default SCF setting INCDFT 0 SCF_CONVERGENCE 8 BASIS 631G* XC_GRID 1 SCF_GUESS SAD !vdw parameters setting DFTVDW_JOBNUMBER 1 DFTVDW_METHOD 1 DFTVDW_PRINT 0 DFTVDW_KAI 800 DFTVDW_USE_ELE_DRV 0 $end
One should note that the XDM option can be used in conjunction with different GGA, metaGGA pure or hybrid functionals, even though the original implementation of Becke and Johnson was in combination with HartreeFock exchange, or with a specific metaGGA exchange and correlation (the BR89 exchange and the BR94 correlation described in previous sections above). For example, encouraging results were obtained using the XDM option with the popular B3LYP [137]. Becke has found more recently that this model can be efficiently combined with the old GGA exchange of Perdew 86 (the P86 exchange option in QChem), and with his hyperGGA functional B05. Using XDM together with PBE exchange plus LYP correlation, or PBE exchange plus BR94 correlation has been also found fruitful.
Recently, Grimme proposed DFTD3 method [141] to improve his previous DFTD method [73] (see Section 4.3.6). Energies and analytic gradients of DFTD3 methods are available in QChem. Grimme’s DFTD3 method can be combined with any of the density functionals available in QChem.
The total DFTD3 energy is given by
(4.63) 
where is the total energy from KSDFT and is the dispersion correction as a sum of two and threebody energies,
(4.64) 
DFTD3 method can be turned on by five keywords, DFT_D, DFT_D3_S6, DFT_D3_RS6, DFT_D3_S8 and DFT_D3_3BODY.
DFT_D
Controls the application of DFTD3 or DFTD scheme.
TYPE:
LOGICAL
DEFAULT:
None
OPTIONS:
FALSE
(or 0) Do not apply the DFTD3 or DFTD scheme
EMPIRICAL_GRIMME3
dispersion correction from Grimme’s DFTD3 method
EMPIRICAL_GRIMME
dispersion correction from Grimme (see Section 4.3.6)
EMPIRICAL_CHG
dispersion correction from Chai and HeadGordon (see Section 4.3.6)
RECOMMENDATION:
NONE
Grimme suggested four scaling factors , , and (see Equation (4) in Ref. Grimme:2010). By fixing , the other three factors were optimized for several functionals (see Table IV in Ref. Grimme:2010). For example, of 1.217 and of 0.722 for PBE, 1.094 and 1.682 for BLYP, 1.261 and 1.703 for B3LYP, 1.532 and 0.862 for PW6B95, 0.892 and 0.909 for BECKE97, and 1.287 and 0.928 for PBE0; these are the QChem default values of and . Otherwise, the default values of , , and are 1.0.
DFT_D3_S6
Controls the strength of dispersion corrections, , in Grimme’s DFTD3 method (see Table IV in Ref. Grimme:2010).
TYPE:
INTEGER
DEFAULT:
1000
OPTIONS:
n
Corresponding to .
RECOMMENDATION:
NONE
DFT_D3_RS6
Controls the strength of dispersion corrections, s, in the Grimme’s DFTD3 method (see Table IV in Ref. Grimme:2010).
TYPE:
INTEGER
DEFAULT:
1000
OPTIONS:
n
Corresponding to .
RECOMMENDATION:
NONE
DFT_D3_S8
Controls the strength of dispersion corrections, , in Grimme’s DFTD3 method (see Table IV in Ref. Grimme:2010).
TYPE:
INTEGER
DEFAULT:
1000
OPTIONS:
n
Corresponding to .
RECOMMENDATION:
NONE
DFT_D3_RS8
Controls the strength of dispersion corrections, , in Grimme’s DFTD3 method (see Equation (4) in Ref. Grimme:2010).
TYPE:
INTEGER
DEFAULT:
1000
OPTIONS:
n
Corresponding to .
RECOMMENDATION:
NONE
The threebody interaction term, mentioned in Ref. Grimme:2010, can also be turned on, if needed.
DFT_D3_3BODY
Controls whether the threebody interaction in Grimme’s DFTD3 method should be applied (see Eq. (14) in Ref. Grimme:2010).
TYPE:
LOGICAL
DEFAULT:
FALSE
OPTIONS:
FALSE
(or 0) Do not apply the threebody interaction term
TRUE
Apply the threebody interaction term
RECOMMENDATION:
NONE
Example 4.37 Applications of B3LYPD3 to a methane dimer.
$comment Geometry optimization, followed by singlepoint calculations using a larger basis set. $end $molecule 0 1 C 0.000000 0.000323 1.755803 H 0.887097 0.510784 1.390695 H 0.887097 0.510784 1.390695 H 0.000000 1.024959 1.393014 H 0.000000 0.001084 2.842908 C 0.000000 0.000323 1.755803 H 0.000000 0.001084 2.842908 H 0.887097 0.510784 1.390695 H 0.887097 0.510784 1.390695 H 0.000000 1.024959 1.393014 $end $rem jobtype opt exchange B3LYP basis 631G* DFT_D EMPIRICAL_GRIMME3 DFT_D3_S6 1000 DFT_D3_RS6 1261 DFT_D3_S8 1703 DFT_D3_3BODY FALSE $end @@@ $molecule READ $end $rem jobtype sp exchange B3LYP basis 6311++G** DFT_D EMPIRICAL_GRIMME3 DFT_D3_S6 1000 DFT_D3_RS6 1261 DFT_D3_S8 1703 DFT_D3_3BODY FALSE $end
The B97XD3 and M06D3 methods are implemented in QChem, which automatically set the empiricallyfitted parameters of the scheme, see Section 4.3.4.4.
The recent advance in doublehybrid density functional theory (DHDFT) [41, 125, 126, 127, 101], has demonstrated its great potential for approaching the chemical accuracy with a computational cost comparable to the secondorder MøllerPlesset perturbation theory (MP2). In a DHDFT, a KohnSham (KS) DFT calculation is performed first, followed by a treatment of nonlocal orbital correlation energy at the level of secondorder MøllerPlesset perturbation theory (MP2) [128]. This MP2 correlation correction includes a a samespin (ss) component, , as well as an oppositespin (os) component, , which are added to the total energy obtained from the KSDFT calculation. Two scaling parameters, and , are introduced in order to avoid doublecounting correlation:
(4.65) 
Among DH functionals, B97X2 [42], a longrange corrected DH functional, is described in Section 4.3.4.3.
There are three keywords for turning on DHDFT as below.
DH
Controls the application of DHDFT scheme.
TYPE:
LOGICAL
DEFAULT:
FALSE
OPTIONS:
FALSE
(or 0) Do not apply the DHDFT scheme
TRUE
(or 1) Apply DHDFT scheme
RECOMMENDATION:
NONE
SSS_FACTOR
Controls the strength of the samespin component of PT2 correlation energy.
TYPE:
INTEGER
DEFAULT:
0
OPTIONS:
n
Corresponding to in Eq. (4.65).
RECOMMENDATION:
NONE
SOS_FACTOR
Controls the strength of the oppositespin component of PT2 correlation energy.
TYPE:
INTEGER
DEFAULT:
0
OPTIONS:
n
Corresponding to in Eq. (4.65).
RECOMMENDATION:
NONE
For example, B2PLYP [73], which involves 53% HartreeFock exchange, 47% Becke 88 GGA exchange, 73% LYP GGA correlation and 27% PT2 orbital correlation, can be called with the following input file. The PT2 correlation energy can also be computed with the efficient resolutionofidentity (RI) methods (see Section 5.5).
Example 4.38 Applications of B2PLYP functional to LiH molecule.
$comment Geometry optimization and frequency calculation on LiH, followed by singlepoint calculations with nonRI and RI approaches. $end $molecule 0 1 H Li H 1.6 $end $rem jobtype opt exchange general correlation mp2 basis ccpvtz DH 1 SSS_FACTOR 270000 !0.27 = 270000/1000000 SOS_FACTOR 270000 !0.27 = 270000/1000000 $end $XC_Functional X HF 0.53 X B 0.47 C LYP 0.73 $end @@@ $molecule READ $end $rem jobtype freq exchange general correlation mp2 basis ccpvtz DH 1 SSS_FACTOR 270000 SOS_FACTOR 270000 $end $XC_Functional X HF 0.53 X B 0.47 C LYP 0.73 $end @@@ $molecule READ $end $rem jobtype sp exchange general correlation mp2 basis ccpvtz DH 1 SSS_FACTOR 270000 SOS_FACTOR 270000 $end $XC_Functional X HF 0.53 X B 0.47 C LYP 0.73 $end @@@ $molecule READ $end $rem jobtype sp exchange general correlation rimp2 basis ccpvtz aux_basis rimp2ccpvtz DH 1 SSS_FACTOR 270000 SOS_FACTOR 270000 $end $XC_Functional X HF 0.53 X B 0.47 C LYP 0.73 $end
A simplification of input is available for the PBE0_DH [102] and PBE0_2 [103] double hybrids from the PBE based linearlyscaled oneparameter doublehybrids, automatically setting the three keywords:
Example 4.39 Applications of PBE0_DH and PBE0_2 functionals to nitrogen molecule.
$comment Single point calculation, nonRI approach. PBE0_DH = PBE + 50% HF + 12.5% MP2 Could also do geometry optimization and frequencies. $end $molecule 0 1 N1 N2 N1 1.0 $end $rem jobtype sp exchange PBE0_DH correlation mp2 basis ccpvtz $end @@@ $comment Single point calculation, expensive nonRI approach. PBE0_2 = PBE + 79.3701% HF + 50% MP2 Could also do geometry optimization and frequencies. $end $molecule 0 1 N1 N2 N1 1.0 $end $rem jobtype sp exchange PBE0_2 correlation rimp2 basis ccpvtz aux_basis rimp2ccpvtz $end
A more detailed gist of one particular class of DH functionals, the XYG3 & XYGJOS functionals follows below thanks to Dr Yousung Jung who implemented these functionals in QChem.
A starting point of these DH functionals is the adiabatic connection formula which provides a rigorous way to define them. One considers an adiabatic path between the fictitious noninteracting KohnSham system (= 0) and the real physical system (= 1) while holding the electron density fixed at its physical state for all :
(4.66) 
where is the exchange correlation potential energy at a coupling strength . If one assumes a linear model of the latter:
(4.67) 
one obtains the popular hybrid functional that includes the HartreeFock exchange (or occupied orbitals) such as B3LYP. If one further uses the GorlingLevy’s perturbation theory (GL2) to define the initial slope at = 0, one obtains the doubly hybrid functional (see Eq. 4.65) that includes MP2 type perturbative terms (PT2) involving virtual KohnSham orbitals:
(4.68) 
In the DH functional XYG3, as implemented in QChem, the B3LYP orbitals are used to generate the density and evaluate the PT2 terms. This is different from P2PLYP, an earlier doubly hybrid functional by Grimme. P2PLYP uses truncated KohnSham orbitals while XYG3 uses converged KS orbitals to evaluate the PT2 terms. The performance of XYG3 is not only comparable to that of the G3 or G2 theory for thermochemistry, but barrier heights and noncovalent interactions, including stacking interactions, are also very well described by XYG3 [101].
The recommended basis set for XYG3 is 6311+G(3df,2p).
Due to the inclusion of PT2 terms, XYG3 or all other forms of doubly hybrid functionals formally scale as the 5th power of system size as in conventional MP2, thereby not applicable to large systems and partly losing DFT’s cost advantages. With the success of SOSMP2 and local SOSMP2 algorithms developed in QChem, the natural extension of XYG3 is to include only oppositespin correlation contributions, ignoring the samespin parts. The resulting DH functional is XYGJOS also implemented in QChem. It has 4 parameters that are optimized using thermochemistry data. This new functional is both accurate (comparable or even slightly better than XYG3) and faster. If the local algorithm is applied, the formal scaling of XYGJOS is cubic, without the locality, it has still 4th order scaling. Recently, XYGJOS becomes the only DH functional with analytical gradient [142].
Example 1: XYG3 calculation of N2. XYG3 invokes automatically the B3LYP calculation first, and use the resulting orbitals for evaluating the MP2type correction terms. One can also use XYG3 combined with RI approximation for the PT2 terms; use EXCHANGE = XYG3RI to do so, along with an appropriate choice of auxiliary basis set.
Example 4.40 XYG3 calculation of N2
$molecule 0 1 N 0.00000000 0.00000000 0.54777500 N 0.00000000 0.00000000 0.54777500 $end $rem exchange xyg3 basis 6311+G(3df,2p) $end
Example 2: XYGJOS calculation of N2. Since it uses the RI approximation by default, one must define the auxiliary basis.
Example 4.41 XYGJOS calculation of N2
$molecule 0 1 N 0.00000000 0.00000000 0.54777500 N 0.00000000 0.00000000 0.54777500 $end $rem exchange xygjos basis 6311+G(3df,2p) aux_basis rimp2ccpVtZ purecart 1111 time_mp2 true $end
Example 3: Local XYGJOS calculation of N2. The same as XYGJOS, except for the use of the attenuated Coulomb metric to solve the RI coefficients. Omega determines the locality of the metric.
Example 4.42 Local XYGJOS calculation of N2
$molecule 0 1 N 0.000 0.000 0.54777500 N 0.000 0.000 0.54777500 $end $rem exchange lxygjos omega 200 basis 6311+G(3df,2p) aux_basis rimp2ccpVtZ purecart 1111 $end
It is well known that no gradientcorrected exchange functional can simultaneously produce the correct contribution to the exchange energy density and exchange potential in the asymptotic region of molecular systems [61]. Existing GGA exchangecorrelation (xc) potentials decay much faster than the correct xc potential in the asymptotic region [143]. Highlying occupied orbitals and lowlying virtual orbitals are therefore too loosely bounded from these GGA functionals, and the minus HOMO energy becomes much less than the exact ionization potential (as required by the exact DFT) [144, 145]. Moreover, response properties could be poorly predicted from TDDFT calculations with GGA functionals [145]. Longrange corrected hybrid DFT (LRCDFT), described in Section 4.3.4, has greatly remedied this situation. However, due to the use of longrange HF exchange, LRCDFT is computationally more expensive than KSDFT with GGA functionals.
To circumvent the above problems, van Leeuwen and Baerends proposed an asymptotically corrected (AC) exchange potential [61]:
(4.69) 
that will reduce to , for an exponentially decaying density, in the asymptotic region of molecular systems, where is the reduced density gradient. The LB94 xc potential is formed by a linear combination of LDA xc potential and the LB exchange potential [61]:
(4.70) 
The parameter was determined by fitting the LB94 xc potential to the beryllium atom. As mentioned in Ref. Casida:1998,Hirata:1999b, for TDDFT and TDDFT/TDA calculations, it is sufficient to include the AC xc potential for groundstate calculations followed by TDDFT calculations with an adiabatic LDA xc kernel. The implementation of LB94 xc potential in QChem thus follows this; using LB94 xc potential for groundstate calculations, followed by TDDFT calculations with an adiabatic LDA xc kernel. This TDLDA/LB94 approach has been widely applied to study excitedstate properties of large molecules in literature.
Since the LB exchange potential does not come from the functional derivative of some exchange functional, we use the LevyPerdew virial relation [148] (implemented in QChem) to obtain its exchange energy:
(4.71) 
LB94_BETA
Set the parameter of LB94 xc potential
TYPE:
INTEGER
DEFAULT:
500
OPTIONS:
Corresponding to .
RECOMMENDATION:
Use default, i.e.,
Example 4.43 Applications of LB94 xc potential to N molecule.
$comment TDLDA/LB94 calculation is performed for excitation energies. $end $molecule 0 1 N 0.0000 0.0000 0.0000 N 1.0977 0.0000 0.0000 $end $rem jobtype = sp exchange = lb94 basis = 6311(2+,2+)G** cis_n_roots = 30 rpa = true $end
Another alternative, proposed by Pan, Fang and Chai [62], was to use a localized version of FermiAmaldi exchangecorrelation functional. The resulting exchange density functional, whose functional derivative has the correct asymptote, can be directly added to any semilocal density functional.
Three variants of this method were proposed in [62]. The simplest variant, the strictly localized FermiAmaldi (LFAs) scheme, is implemented in QChem 4.3, for molecules consisting of atoms with .
Example 4.44 LFAsPBE singlepoint TDDFT calculation with water molecule
$comment Use LFAsPBE potential for groundstate calculations, followed by TDDFT calculations with an adiabatic PBE xc kernel. $end $molecule 0 1 O H1 O oh H2 O oh H1 hoh oh = 1.0 hoh = 110.0 $end. $rem jobtype sp exchange gen basis 6311(2+,2+)G** cis_n_roots 30 rpa true $end $xc_functional X PBE 1.0 C PBE 1.0 X LFAs 1.0 $end
From the perspective of perturbation theory, Chai and Chen proposed [149] a systematic procedure for the evaluation of the derivative discontinuity (DD) of the exchangecorrelation energy functional in KohnSham density functional theory, wherein the exact DD can in principle be obtained by summing up all the perturbation corrections to infinite order. Truncation of the perturbation series at low order yields an efficient scheme for obtaining the approximate DD. In particular, the firstorder correction term is equivalent to the frozenorbital approximation method.
The firstorder scheme implemented in QChem 4.3 supports only local or GGA functionals, does not yet support metaGGA, HF, hybrids or nonlocal correlation.
FOA_FUNDGAP
Compute the frozenorbital approximation of the fundamental gap.
TYPE:
Boolean
DEFAULT:
false
OPTIONS:
false
(default) do not compute FOA DD and fundamental gap
true
compute and print FOA fundamental gap information. Implies KS_GAP_PRINT.
RECOMMENDATION:
Use in conjuction with KS_GAP_UNIT if true.
Example 4.45 frozenorbital approximation of DD with PBE and LFAsPBE functionals on carbon atom
$comment Frozenorbital derivative discontinuity, C atom, PBE $end $molecule 0 3 C $end $rem jobtype sp basis 631G* method PBE foa_fundgap true ks_gap_unit 1 ! print gap info in eV $end @@@ $comment with LFAsPBE functional instead $end $molecule READ $end $rem jobtype sp basis 631G* scf_guess READ exchange gen foa_fundgap true ks_gap_unit 1 $end $xc_functional X PBE 1.0 X LFAs 1.0 C PBE 1.0 $end
Aiming to study the groundstate properties of large strongly correlated systems with minimum computational complexity, Chai recently developed thermallyassistedoccupation density functional theory (TAODFT) [150]. Unlike conventional ab initio multireference methods, the computational complexity of TAODFT increases very insignificantly with the size of the active space (i.e., an active space restriction is not needed for TAODFT calculations), showing that TAODFT can be very promising for the study of large polyradical systems. In contrast to KSDFT, TAODFT is a density functional theory with fractional orbital occupations produced by the FermiDirac distribution (controlled by a fictitious temperature ). However, existing XC functionals (e.g., LDA and GGAs) in KSDFT may also be adopted in TAODFT [151]. TAODFT has similar computational cost as KSDFT for singlepoint energy calculations and analytical nuclear gradients, and reduces to KSDFT in the absence of strong static correlation effects.
There are several $rem variables that are used for TAODFT.
TAO_DFT
Controls whether to use TAODFT.
TYPE:
Boolean
DEFAULT:
false
OPTIONS:
false
do not use TAODFT
true
use TAODFT
RECOMMENDATION:
NONE
TAO_DFT_THETA
value of in TAODFT.
TYPE:
INTEGER
DEFAULT:
7
OPTIONS:
(hartrees), where is the value of TAO_DFT_THETA_NDP
RECOMMENDATION:
NONE
TAO_DFT_THETA_NDP
value of in TAODFT.
TYPE:
INTEGER
DEFAULT:
3
OPTIONS:
(hartrees), where is the value of TAO_DFT_THETA
RECOMMENDATION:
NONE
Note that setting TAO_DFT_THETA=0 recovers KSDFT [150].
In addition to the XC functional, is needed in TAODFT. Currently, LDA [150] and GEA (gradient expansion approximation, [151]) are available in QChem.
Example 4.46 TAOLDA calcuation on Be atom
$molecule 0 1 Be $end $rem jobtype sp basis 631G* exchange gen tao_dft true tao_dft_theta 7 ! default, theta=7 mhartree tao_dft_theta_ndp 3 ! default $end $xc_functional X S 1.0 C PW92 1.0 X ETheta_LDA 1.0 $end
Example 4.47 TAOPBE, spinrestricted calculation on stretched N
$comment Spinrestricted calculation on stretched N2 $end $molecule 0 1 N1 N2 N1 5 $end $rem jobtype sp basis 631G* exchange gen tao_dft true tao_dft_theta 40 ! theta = 40 mhartree tao_dft_theta_ndp 3 $end $xc_functional X PBE 1.0 C PBE 1.0 X ETheta_LDA 1.0 $end
Example 4.48 TAOPBE, spinunrestricted calculation on stretched N
ETheta_LDA may also be replaced with ETheta_GEA, see [151].$comment Spinunrestricted optimization calculation on stretched N2 $end $molecule 0 1 N1 N2 N1 5 $end $rem jobtype opt unrestricted true basis 631G* exchange gen tao_dft true tao_dft_theta 40 ! theta = 40 minihartrees tao_dft_theta_ndp 3 ! can omit this line scf_guess gwh SCF_GUESS_MIX 3 ! mix in 30% lumo in alpha to break symmetry $end $xc_functional X PBE 1.0 C PBE 1.0 X ETheta_LDA 1.0 $end
In practical DFT calculations, the forms of the approximate exchangecorrelation functionals used are quite complicated, such that the required integrals involving the functionals generally cannot be evaluated analytically. QChem evaluates these integrals through numerical quadrature directly applied to the exchangecorrelation integrand (i.e., no fitting of the XC potential in an auxiliary basis is done). QChem provides a standard quadrature grid by default which is sufficient for most purposes.
The quadrature approach in QChem is generally similar to that found in many DFT programs. The multicenter XC integrals are first partitioned into “atomic” contributions using a nuclear weight function. QChem uses the nuclear partitioning of Becke [152], though without the atomic size adjustments". The atomic integrals are then evaluated through standard onecenter numerical techniques.
Thus, the exchangecorrelation energy is obtained as
(4.72) 
where the first summation is over the atoms and the second is over the numerical quadrature grid points for the current atom. The function is the exchangecorrelation functional. The are the quadrature weights, and the grid points are given by
(4.73) 
where is the position of nucleus , with the defining a suitable onecenter integration grid, which is independent of the nuclear configuration.
The singlecenter integrations are further separated into radial and angular integrations. Within QChem, the radial part is usually treated by the EulerMaclaurin scheme proposed by Murry et al. [153]. This scheme maps the semiinfinite domain and applies the extended trapezoidal rule to the transformed integrand. Recently Gill and Chien [154] proposed a radial scheme based on a Gaussian quadrature on the interval with weight function . This scheme is exact for integrands that are a linear combination of a geometric sequence of exponential functions, and is therefore well suited to evaluating atomic integrals. The authors refer to this scheme as MultiExp.
Angular quadrature rules may be characterized by their degree, which is the highest degree of spherical harmonics for which the formula is exact, and their efficiency, which is the number of spherical harmonics exactly integrated per degree of freedom in the formula. QChem supports the following types of angular grids:
Lebedev These are specially constructed grids for quadrature on the surface of a sphere [155, 156, 157] based on the octahedral group. Lebedev grids of the following degrees are available:
Degree 
3rd 
5th 
7th 
9th 
11th 
15th 
17th 
19th 
23rd 
29th 
Points 
6 
18 
26 
38 
50 
86 
110 
146 
194 
302 
Additional grids with the following number of points are also available: 74, 170, 230, 266, 350, 434, 590, 770, 974, 1202, 1454, 1730, 2030, 2354, 2702, 3074, 3470, 3890, 4334, 4802, 5294. Lebedev grids typically have efficiencies near one, with efficiencies greater than one in some cases.
GaussLegendre These are spherical product rules separating the two angular dimensions and . Integration in the dimension is carried out with a Gaussian quadrature rule derived from the Legendre polynomials (orthogonal on with weight function unity), while the integration is done with equally spaced points.
A GaussLegendre grid is selected by specifying the total number of points, , to be used for the integration. This gives a grid with points, points, and a degree of .
In contrast with Lebedev grids, GaussLegendre grids have efficiency of only 2/3 (hence more GaussLegendre points are required to attain the same accuracy as Lebedev). However, since GaussLegendre grids of general degree are available, this is a convenient mechanism for achieving arbitrary accuracy in the angular integration if desired.
Combining these radial and angular schemes yields an intimidating selection of threedimensional quadratures. In practice, is it useful to standardize the grids used in order to facilitate the comparison of calculations at different levels of theory.
Both the SG0 [158] and SG1 [159] standard quadrature grids were designed to yield the performance of a large, accurate quadrature grid, but with as few points as possible for the sake of computational efficiency. This is accomplished by reducing the number of angular points in regions where sophisticated angular quadrature is not necessary, such as near the nuclei where the charge density is nearly spherically symmetric, while retaining large numbers of angular points in the valence region where angular accuracy is critical.
The SG0 grid was derived in this fashion from a MultiExpLebedev(23,170), (i.e., 23 radial points and 170 angular points per radial point). This grid was pruned whilst ensuring the error in the computed exchange energies for the atoms and a selection of small molecules was not larger than the corresponding error associated with SG1. In our evaluation, the RMS error associated with the atomization energies for the molecules in the G1 data set is 72 microhartrees. While relative energies are expected to be reproduced well by this scheme, if absolute energies are being sought, a larger grid is recommended.
The SG0 grid is implemented in QChem from H to micro Hartrees, excepted He and Na; in this scheme, each atom has around 1400point, and SG1 is used for those their SG0 grids have not been defined. It should be noted that, since the SG0 grid used for H has been reoptimized in this version of QChem (version 3.0), quantities calculated in this scheme may not reproduce those generated by the last version (version 2.1).
The SG1 grid is derived from a EulerMaclaurinLebedev(50,194) grid (i.e., 50 radial points, and 194 angular points per radial point). This grid has been found to give numerical integration errors of the order of 0.2 kcal/mol for mediumsized molecules, including particularly demanding test cases such as isomerization energies of alkanes. This error is deemed acceptable since it is significantly smaller than the accuracy typically achieved by quantum chemical methods. In SG1 the total number of points is reduced to approximately 1/4 of that of the original EML(50,194) grid, with SG1 generally giving the same total energies as EML(50,194) to within a few microhartrees (0.01 kcal/mol). Therefore, the SG1 grid is relatively efficient while still maintaining the numerical accuracy necessary for chemical reliability in the majority of applications.
Both the SG0 and SG1 grids were optimized so that the error in the energy when using the grid did not exceed a target threshold. For single point calculations this criterion is appropriate. However, derivatives of the energy can be more sensitive to the quality of the integration grid, and it is recommended that a larger grid be used when calculating these. Special care is required when performing DFT vibrational calculations as imaginary frequencies can be reported if the grid is inadequate. This is more of a problem with lowfrequency vibrations. If imaginary frequencies are found, or if there is some doubt about the frequencies reported by QChem, the recommended procedure is to perform the calculation again with a larger grid and check for convergence of the frequencies. Of course the geometry must be reoptimized, but if the existing geometry is used as an initial guess, the geometry optimization should converge in only a few cycles.
Whenever QChem calculates numerical density functional integrals, the electron density itself is also integrated numerically as a test on the quality of the quadrature formula used. The deviation of the numerical result from the number of electrons in the system is an indication of the accuracy of the other numerical integrals. If the relative error in the numerical electron count reaches 0.01%, a warning is printed; this is an indication that the numerical XC results may not be reliable. If the warning appears at the first SCF cycle, it is probably not serious, because the initialguess density matrix is sometimes not idempotent, as is the case with the SAD guess and the density matrix taken from a different geometry in a geometry optimization. If that is the case, the problem will be corrected as the idempotency is restored in later cycles. On the other hand, if the warning is persistent to the end of SCF iterations, then either a finer grid is needed, or choose an alternative method for generating the initial guess.
Users should be aware, however, of the potential flaws that have been discovered in some of the grids currently in use. Jarecki and Davidson [160], for example, have recently shown that correctly integrating the density is a necessary, but not sufficient, test of grid quality.
By default, QChem will estimate the magnitude of various XC contributions on the grid and eliminate those determined to be numerically insignificant. QChem uses specially developed cutoff procedures which permits evaluation of the XC energy and potential in only work for large molecules, where is the size of the system. This is a significant improvement over the formal scaling of the XC cost, and is critical in enabling DFT calculations to be carried out on very large systems. In very rare cases, however, the default cutoff scheme can be too aggressive, eliminating contributions that should be retained; this is almost always signaled by an inaccurate numerical density integral. An example of when this could occur is in calculating anions with multiple sets of diffuse functions in the basis. As mentioned above, when an inaccurate electron count is obtained, it maybe possible to remedy the problem by increasing the size of the quadrature grid.
Finally we note that early implementations of quadraturebased KohnSham DFT employing standard basis sets were plagued by lack of rotational invariance. That is, rotation of the system yielded a significantly energy change. Clearly, such behavior is highly undesirable. Johnson et al. rectified the problem of rotational invariance by completing the specification of the grid procedure [161] to ensure that the computed XC energy is the same for any orientation of the molecule in any Cartesian coordinate system.
For most DFT jobs two $rem variables are required: METHOD and BASIS. To request a specific pair of exchange and correlation functionals, the EXCHANGE and CORRELATION keywords can be used instead of METHOD. In addition, all of the basic input options discussed for HartreeFock calculations in Section 4.2.3, and the extended options discussed in Section 4.2.4 are all valid for DFT calculations. Below we list only the basic DFTspecific options (keywords).
METHOD
Specifies the exchangecorrelation functional.
TYPE:
STRING
DEFAULT:
No default
OPTIONS:
NAME
Use METHOD = NAME, where NAME is
one of the DFT methods listed in Table 4.2.
RECOMMENDATION:
Consult the literature to guide your selection.
METHOD = 
Exchange 
Correlation 
HF 
Hartree–Fock 

LDA 
Slater 
VoskoWilkNusair parameterization #5 (VWN) 
BLYP 
Becke 1988 (B) 
LeeYangParr (LYP) 
BP86 
B 
Perdew 1986 (P86) 
BP86VWN 
B 
P86VWN 
PBE 
PerdewBurkeErnzerhof 1996 (PBE) 

PBE0 
PBE hybrid with 25% HF exchange 

PBE50 
PBE hybrid with 50% HF exchange 

B3LYP 
B3LYP hybrid 

BHHLYP 
Modified BLYP functional with 50% HF exchange 

B5050LYP 
Modified B3LYP functional with 50% HF exchange 

B97D 
Grimme’s modified B97 with empirical dispersion 

CAMB3LYP 
Coulombattenuated B3LYP 

X3LYP 
X3LYP from Xu and Goddard 

EDF1 
EDF1 

EDF2 
EDF2 

LRCwPBE 
Longrangecorrected pure PBE 

LRCwPBEH 
Longrangecorrected hybrid PBE 

M05 
M05 hybrid 

M052X 
M052X hybrid 

M06 
M06 hybrid 

M06L 
M06L hybrid 

M06HF 
M06HF hybrid 

M062X 
M062X hybrid 

M08HX 
M08HX hybrid 

M08SO 
M08SO hybrid 

M11 
M11 longrangecorrected hybrid 

M11L 
M11L hybrid 

wB97 
B97 longrangecorrected hybrid 

wB97X 
B97X longrangecorrected hybrid 

wB97XD 
B97XD longrangecorrected hybrid with dispersion 

corrections 

wB97XD3 
B97XD3 longrangecorrected hybrid with dispersion 

corrections 

wM05D 
M05D longrangecorrected hybrid with dispersion 

corrections 

wM06D3 
M06D3 longrangecorrected hybrid with dispersion 

corrections 
EXCHANGE
Specifies the exchange functional or exchangecorrelation functional for hybrid.
TYPE:
STRING
DEFAULT:
No default exchange functional
OPTIONS:
NAME
Use EXCHANGE = NAME, where NAME is
one of the exchange functionals listed in Table 4.3.
RECOMMENDATION:
Consult the literature to guide your selection.
CORRELATION
Specifies the correlation functional.
TYPE:
STRING
DEFAULT:
None
No correlation.
OPTIONS:
None
No correlation
VWN
VoskoWilkNusair parameterization #5
LYP
LeeYangParr (LYP)
PW91, PW
GGA91 (PerdewWang)
PW92
LSDA 92 (Perdew and Wang) [46]
PK09
LSDA (ProynovKong) [47]
LYP(EDF1)
LYP(EDF1) parameterization
Perdew86, P86
Perdew 1986
PZ81, PZ
PerdewZunger 1981
PBE
PerdewBurkeErnzerhof 1996
TPSS
The correlation component of the TPSS functional
B94
Becke 1994 correlation in fully analytic form
B95
Becke 1995 correlation
B94hyb
Becke 1994 correlation as above, but readjusted for use only within
the hybrid scheme BR89B94hyb
PK06
ProynovKong 2006 correlation (known also as “tLap”)
(B88)OP
OP correlation [80], optimized for use with B88 exchange
(PBE)OP
OP correlation [80], optimized for use with PBE exchange
Wigner
Wigner
RECOMMENDATION:
Consult the literature to guide your selection.
EXCHANGE = 
Description 
HF 
Fock exchange 
Slater, S 
Slater (Dirac 1930) 
ETheta_LDA, ETheta_LSDA 
TAODFT local density approximation for [150] 
(use in conjunction with another exchange functional) 

ETheta_GEA 
TAODFT gradient expansion approximation for [151] 
(use in conjunction with another exchange functional) 

Becke86, B86 
Becke 1986 
Becke, B, B88 
Becke 1988 
muB88 
Shortrange Becke exchange, as formulated by Song et al. [114] 
Gill96, Gill 
Gill 1996 
GG99 
Gilbert and Gill, 1999 
Becke(EDF1), B(EDF1) 
Becke (uses EDF1 parameters) 
PW86, 
PerdewWang 1986 
rPW86, 
Refitted PW86 [51] 
PW91, PW 
PerdewWang 1991 
mPW1PW91 
modified PW91 
mPW1LYP 
modified PW91 exchange + LYP correlation 
mPW1PBE 
modified PW91 exchange + PBE correlation 
mPW1B95 
modified PW91 exchange + B97 correlation 
mPWB1K 
mPWB1K 
PBE 
PerdewBurkeErnzerhof 1996 
AK13 
Armiento and Kümmel, 2013 [63] 
TPSS 
The nonempirical exchangecorrelation scheme of Tao, 
Perdew, Staroverov, and Scuseria (requires also that the user 

specify “TPSS” for correlation) 

TPSSH 
The hybrid version of TPSS (with no input line for correlation) 
PBE0, PBE1PBE 
PBE hybrid with 25% HF exchange 
PBE50 
PBE hybrid with 50% HF exchange 
revPBE 
revised PBE exchange [60] 
PBEOP 
PBE exchange + oneparameter progressive correlation 
wPBE 
Shortrange PBE exchange, as formulated by Henderson et al. [115] 
muPBE 
Shortrange PBE exchange, as formulated by Song et al. [114] 
LRCwPBEPBE 
longrange corrected pure PBE 
LRCwPBEhPBE 
longrange corrected hybrid PBE 
B1B95 
Becke hybrid functional with B97 correlation 
B97 
Becke97 XC hybrid 
B971 
Becke97 reoptimized by Hamprecht et al. (1998) 
B972 
Becke971 optimized further by Wilson et al. (2001) 
B97D 
Grimme’s modified B97 with empirical dispersion 
B3PW91, Becke3PW91, B3P 
B3PW91 hybrid 
B3LYP, Becke3LYP 
B3LYP hybrid 
B3LYP5 
B3LYP based on correlation functional #5 of Vosko, Wilk, 
and Nusair (rather than their functional #3) 

CAMB3LYP 
Coulombattenuated B3LYP 
HCO3LYP 
O3LYP from Handy 
X3LYP 
X3LYP from Xu and Goddard 
B5050LYP 
modified B3LYP functional with 50% HartreeFock exchange 
BHHLYP 
modified BLYP functional with 50% HartreeFock exchange 
B3P86 
B3 exchange and Perdew 86 correlation 
B3PW91 
B3 exchange and GGA91 correlation 
B3tLAP 
Hybrid Becke exchange and PK06 correlation 
HCTH 
HCTH hybrid 
HCTH120 
HCTH120 hybrid 
HCTH147 
HCTH147 hybrid 
HCTH407 
HCTH407 hybrid 
BOP 
B88 exchange + oneparameter progressive correlation 
EDF1 
EDF1 
EDF2 
EDF2 
VSXC 
VSXC metaGGA, not a hybrid 
BMK 
BMK hybrid 
M05 
M05 hybrid 
M052X 
M052X hybrid 
M06L 
M06L hybrid 
M06HF 
M06HF hybrid 
M06 
M06 hybrid 
M062X 
M062X hybrid 
M08HX 
M08HX hybrid 
M08SO 
M08SO hybrid 
M11L 
M11L hybrid 
M11 
M11 longrange corrected hybrid 
SOGGA 
SOGGA hybrid 
SOGGA11 
SOGGA11 hybrid 
SOGGA11X 
SOGGA11X hybrid 
BR89 
BeckeRoussel 1989 represented in analytic form 
BR89B94h 
Hybrid BR89 exchange and B94hyb correlation 
omegaB97 
B97 longrange corrected hybrid 
omegaB97X 
B97X longrange corrected hybrid 
omegaB97XD 
B97XD longrange corrected hybrid with dispersion 
corrections 

omegaB97X2(LP) 
B97X2(LP) longrange corrected doublehybrid 
omegaB97X2(TQZ) 
B97X2(TQZ) longrange corrected doublehybrid 
MCY2 
The MCY2 hyperGGA exchangecorrelation (with no 
input line for correlation) 

B05 
Full exactexchange hyperGGA functional of Becke 05 with 
RI approximation for the exactexchange energy density 

BM05 
Modified B05 hyperGGA scheme with RI approximation for 
the exactexchange energy density used as a variable. 

XYG3 
XYG3 doublehybrid functional 
XYGJOS 
XYGJOS doublehybrid functional 
LXYGJOS 
LXYGJOS doublehybrid functional with localized MP2 
PBE0_DH, PBE0_2 
PBE double hybrid functionals, requires setting 
CORRELATION to an MP2 implementation 

General, Gen 
User defined combination of K, X and C (refer to the next 
section) 
NL_CORRELATION
Specifies a nonlocal correlation functional that includes nonempirical dispersion.
TYPE:
STRING
DEFAULT:
None
No nonlocal correlation.
OPTIONS:
None
No nonlocal correlation
vdWDF04
the nonlocal part of vdWDF04
vdWDF10
the nonlocal part of vdWDF10 (aka vdWDF2)
VV09
the nonlocal part of VV09
VV10
the nonlocal part of VV10
RECOMMENDATION:
Do not forget to add the LSDA correlation (PW92 is recommended) when using vdWDF04, vdWDF10, or VV09. VV10 should be used with PBE correlation. Choose exchange functionals carefully: HF, rPW86, revPBE, and some of the LRC exchange functionals are among the recommended choices.
NL_VV_C
Sets the parameter in VV09 and VV10. This parameter is fitted to asymptotic van der Waals coefficients.
TYPE:
INTEGER
DEFAULT:
89
for VV09
No default
for VV10
OPTIONS:
Corresponding to
RECOMMENDATION:
is recommended when a semilocal exchange functional is used. is recommended when a longrange corrected (LRC) hybrid functional is used. See further details in Ref. [135].
NL_VV_B
Sets the parameter in VV10. This parameter controls the short range behavior of the nonlocal correlation energy.
TYPE:
INTEGER
DEFAULT:
No default
OPTIONS:
Corresponding to
RECOMMENDATION:
The optimal value depends strongly on the exchange functional used. is recommended for rPW86. See further details in Ref. [135].
FAST_XC
Controls direct variable thresholds to accelerate exchange correlation (XC) in DFT.
TYPE:
LOGICAL
DEFAULT:
FALSE
OPTIONS:
TRUE
Turn FAST_XC on.
FALSE
Do not use FAST_XC.
RECOMMENDATION:
Caution: FAST_XC improves the speed of a DFT calculation, but may occasionally cause the SCF calculation to diverge.
XC_GRID
Specifies the type of grid to use for DFT calculations.
TYPE:
INTEGER
DEFAULT:
1
SG1 hybrid
OPTIONS:
0
Use SG0 for H, C, N, and O, SG1 for all other atoms.
1
Use SG1 for all atoms.
2
Low Quality.
The first six integers correspond to radial points and the second six
integers correspond to angular points where possible numbers of Lebedev
angular points are listed in section 4.3.13.
The first six integers correspond to radial points and the second six
integers correspond to angular points where the number of GaussLegendre
angular points .
RECOMMENDATION:
Use default unless numerical integration problems arise. Larger grids may be required for optimization and frequency calculations.
XC_SMART_GRID
Uses SG0 (where available) for early SCF cycles, and switches to the (larger) grid specified by XC_GRID (which defaults to SG1, if not otherwise specified) for final cycles of the SCF.
TYPE:
LOGICAL
DEFAULT:
FALSE
OPTIONS:
TRUE/FALSE
RECOMMENDATION:
The use of the smart grid can save some time on initial SCF cycles.
NL_GRID
Specifies the grid to use for nonlocal correlation.
TYPE:
INTEGER
DEFAULT:
1
SG1 grid
OPTIONS:
Same as for XC_GRID
RECOMMENDATION:
Use default unless computational cost becomes prohibitive, in which case SG0 may be used. XC_GRID should generally be finer than NL_GRID.
Example 4.49 QChem input for a DFT single point energy calculation on water.
$comment BLYP/STO3G water single point calculation $end $molecule 0 1 O H1 O oh H2 O oh H1 hoh oh = 1.2 hoh = 120.0 $end. $rem METHOD BLYP Becke88 exchange + LYP correlation BASIS sto3g Basis set $end
The format for entering userdefined exchangecorrelation density functionals is one line for each component of the functional. Each line requires three variables: the first defines whether the component is an exchange or correlation functional by declaring an or , respectively. The second variable is the symbolic representation of the functional as used for the EXCHANGE and CORRELATION $rem variables. The final variable is a real number corresponding to the contribution of the component to the functional. HartreeFock exchange contributions (required for hybrid density functionals) can be entered using only two variables (, for HF exchange) followed by a real number.
$xc_functional X exchange_symbol coefficient X exchange_symbol coefficient ... C correlation_symbol coefficient C correlation_symbol coefficient ... K coefficient $end
Note: (1) Coefficients are real.
(2) A userdefined functional does not require all , and components.
Examples of userdefined XCs: these are XC options that for the time being can only be invoked via the user defined XC input section:
Example 4.50 QChem input of water with B3tLap.
$comment water with B3tLap $end $molecule 0 1 O H1 O oh H2 O oh H1 hoh oh = 0.97 hoh = 120.0 $end $rem EXCHANGE gen CORRELATION none XC_GRID 000120000194 ! recommended for high accuracy BASIS G3LARGE ! recommended for high accuracy THRESH 14 ! recommended for high accuracy and better convergence $end $xc_functional X Becke 0.726 X S 0.0966 C PK06 1.0 K 0.1713 $end
Example 4.51 QChem input of water with BR89B94hyb.
$comment water with BR89B94hyb $end $molecule 0 1 O H1 O oh H2 O oh H1 hoh oh = 0.97 hoh = 120.0 $end $rem EXCHANGE gen CORRELATION none XC_GRID 000120000194 ! recommended for high accuracy BASIS G3LARGE ! recommended for high accuracy THRESH 14 ! recommended for high accuracy and better convergence $end $xc_functional X BR89 0.846 C B94hyb 1.0 K 0.154 $end
More specific is the use of the RIB05 and RIPSTS functionals. In this release we offer only singlepoint SCF calculations with these functionals. Both options require a converged LSD or HF solution as initial guess, which greatly facilitates the convergence. It also requires specifying a particular auxiliary basis set:
Example 4.52 QChem input of H using RIB05.
$comment H2, example of SP RIB05. First do a wellconverged LSD, G3LARGE is the basis of choice for good accuracy. The input lines purecar 222 SCF_GUESS CORE are obligatory for the time being here. $end $molecule 0 1 H 0. 0. 0.0 H 0. 0. 0.7414 $end $rem JOBTYPE SP SCF_GUESS CORE METHOD LDA BASIS G3LARGE purcar 222 THRESH 14 MAX_SCF_CYCLES 80 PRINT_INPUT TRUE INCDFT FALSE XC_GRID 000128000302 SYM_IGNORE TRUE SYMMETRY FALSE SCF_CONVERGENCE 9 $end @@@@ $comment For the time being the following input lines are obligatory: purcar 22222 AUX_BASIS riB05ccpvtz dft_cutoffs 0 1415 1 MAX_SCF_CYCLES 0 JOBTYPE SP $end $molecule READ $end $rem JOBTYPE SP SCF_GUESS READ EXCHANGE B05 ! EXCHANGE PSTS ! use this line for RIPSTS purcar 22222 BASIS G3LARGE AUX_BASIS riB05ccpvtz ! The aux basis for RIB05 and RIPSTS THRESH 14 PRINT_INPUT TRUE INCDFT FALSE XC_GRID 000128000302 SYM_IGNORE TRUE SYMMETRY FALSE MAX_SCF_CYCLES 0 dft_cutoffs 0 1415 1 $end
Besides postLSD, the RIB05 option can be used as postHartreeFock method as well, in which case one first does a wellconverged HF calculation and uses it as a guess read in the consecutive RIB05 run.