Q-Chem 5.0 User’s Manual

4.4 Density Functional Theory

4.4.1 Introduction

DFT [20, 21, 22, 23] has emerged as an accurate, alternative first-principles approach to quantum mechanical molecular investigations. DFT calculations account for the overwhelming majority of all quantum chemistry calculations, not only because of its proven chemical accuracy, but also because of its relatively low computational expense, comparable to Hartree-Fock theory but with treatment of electron correlation that is neglected in a HF calculation. These two features suggest that DFT is likely to remain a leading method in the quantum chemist’s toolkit well into the future. Q-Chem contains fast, efficient and accurate algorithms for all popular density functionals, making calculations on large molecules possible and practical.

DFT is primarily a theory of electronic ground state structures based on the electron density, $\rho (\ensuremath{\mathbf{r}})$, as opposed to the many-electron wave function, $\Psi (\ensuremath{\mathbf{r}}_1,\ldots ,\ensuremath{\mathbf{r}}_ N)$. (It’s excited-state extension, time-dependent DFT, is discussed in Section 6.3.) There are a number of distinct similarities and differences between traditional wave function approaches and modern DFT methodologies. First, the essential building blocks of the many-electron wave function $\Psi $ are single-electron orbitals, which are directly analogous to the Kohn-Sham orbitals in the DFT framework. Second, both the electron density and the many-electron wave function tend to be constructed via a SCF approach that requires the construction of matrix elements that are conveniently very similar.

However, traditional ab initio approaches using the many-electron wave function as a foundation must resort to a post-SCF calculation (Chapter 5) to incorporate correlation effects, whereas DFT approaches incorporate correlation at the SCF level. Post-SCF methods, such as perturbation theory or coupled-cluster theory are extremely expensive relative to the SCF procedure. On the other hand, while the DFT approach is exact in principle, in practice it relies on modeling an unknown 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 a systematically-improvable level of accuracy, but can be computationally demanding, whereas DFT approaches offer a practical route, but the theory is currently incomplete.

4.4.2 Kohn-Sham Density Functional Theory

The density functional theory by Hohenberg, Kohn, and Sham [24, 25] stems from the original work of Dirac [26], 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 where 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 Kohn-Sham formalism [25], the ground state electronic energy, $E$, can be written as

  \begin{equation} \label{eq:DFT_ energy_ components} E = E_{\mathrm{T}}^{} + E_{\mathrm{V}}^{} + E_{\mathrm{J}}^{} + E_{\mathrm{XC}}^{} \end{equation}   (4.34)

where $E_{\mathrm{T}}^{}$ is the kinetic energy, $E_{\mathrm{V}}^{}$ is the electron–nuclear interaction energy, $E_{\mathrm{J}}^{}$ is the Coulomb self-interaction of the electron density, $\rho (\ensuremath{\mathbf{r}})$ and $E_{\mathrm{XC}}^{}$ is the exchange-correlation energy. Adopting an unrestricted format, the $\alpha $ and $\beta $ total electron densities can be written as

  $\displaystyle \label{eq429} \begin{aligned} \rho _{\alpha }(\ensuremath{\mathbf{r}}) & = \sum _{i=1}^{n_\alpha } |\psi _ i^{\alpha }|^2 \\ \rho _{\beta }(\ensuremath{\mathbf{r}}) & = \sum _{i=1}^{n_\beta } |\psi _ i^{\beta } |^2 \end{aligned} $   (4.35)

where $n_{\alpha }$ and $n_{\beta }$ are the number of alpha and beta electron respectively, and $\psi _ i$ are the Kohn-Sham orbitals. Thus, the total electron density is

  \begin{equation} \label{eq430} \rho (\ensuremath{\mathbf{r}}) =\rho _\alpha (\ensuremath{\mathbf{r}}) + \rho _\beta (\ensuremath{\mathbf{r}}) \end{equation}   (4.38)

Within a finite basis set [27], the density is represented by

  \begin{equation} \label{eq431} \rho (\ensuremath{\mathbf{r}}) = \sum _{\mu \nu } P_{\mu \nu } \phi _\mu (\ensuremath{\mathbf{r}}) \phi _\nu (\ensuremath{\mathbf{r}}) \; , \end{equation}   (4.39)

where the $P_{\mu \nu }$ are the elements of the one-electron density matrix; see Eq. eq:P(mu,nu) in the discussion of Hartree-Fock theory. The various energy components in Eq. eq:DFT_energy_components can now be written

  $\displaystyle \label{eq432} E_{\mathrm{T}}^{} $ $\displaystyle = $ $\displaystyle \sum _{i=1}^{n_\alpha } \left\langle \psi _ i^\alpha \left|-\frac{1}{2}\nabla ^2\right|\psi _ i^{\alpha } \right\rangle + \sum _{i=1}^{n_\beta } \left\langle \psi _ i^\beta \left|-\frac{1}{2}\nabla ^2\right|\psi _ i^{\beta } \right\rangle \nonumber $    
  $\displaystyle $ $\displaystyle = $ $\displaystyle \sum _{\mu \nu } P_{\mu \nu } \left\langle \phi _\mu (\ensuremath{\mathbf{r}}) \left|-\frac{1}{2}\nabla ^2\right| \phi _\nu (\ensuremath{\mathbf{r}}) \right\rangle $   (4.40)
  $\displaystyle E_{\mathrm{V}}^{} $ $\displaystyle = $ $\displaystyle -\sum _{A=1}^ M Z_ A\int \frac{\rho (\ensuremath{\mathbf{r}})}{|\ensuremath{\mathbf{r}}-\ensuremath{\mathbf{R}}_ A|} d\ensuremath{\mathbf{r}} \nonumber $    
  $\displaystyle $ $\displaystyle = $ $\displaystyle -\sum _{\mu \nu } P_{\mu \nu } \sum _ A \left\langle \phi _\mu (\ensuremath{\mathbf{r}}) \left| \frac{Z_ A}{|\ensuremath{\mathbf{r}}-\ensuremath{\mathbf{R}}_ A|} \right| \phi _\nu (\ensuremath{\mathbf{r}}) \right\rangle $   (4.41)
  $\displaystyle E_{\mathrm{J}}^{} $ $\displaystyle = $ $\displaystyle \frac{1}{2}\left\langle \rho (\ensuremath{\mathbf{r}}_1) \left| \frac{1}{|\ensuremath{\mathbf{r}}_1 - \ensuremath{\mathbf{r}}_2|} \right| \rho (\ensuremath{\mathbf{r}}_2) \right\rangle \nonumber $    
  $\displaystyle $ $\displaystyle = $ $\displaystyle \frac{1}{2}\sum _{\mu \nu } \sum _{\lambda \sigma } P_{\mu \nu } P_{\lambda \sigma } \left(\mu \nu \vert \lambda \sigma \right) $   (4.42)
  $\displaystyle E_{\mathrm{XC}}^{} $ $\displaystyle = $ $\displaystyle \int f\left[\rho (\ensuremath{\mathbf{r}}),\nabla \rho (\ensuremath{\mathbf{r}}),\ldots \right] d\ensuremath{\mathbf{r}} \label{E_ XC} $   (4.43)

Minimizing $E$ with respect to the unknown Kohn-Sham orbital coefficients yields a set of matrix equations exactly analogous to Pople-Nesbet equations of the UHF case, Eq. eq:Pople-Nesbet, but with modified Fock matrix elements [cf. Eq. eq:F(mu,nu)]

  $\displaystyle \label{eq437} \begin{aligned} F_{\mu \nu }^\alpha & = H_{\mu \nu }^{\mathrm{core}} + J_{\mu \nu } - F_{\mu \nu }^{\mathrm{XC}\alpha } \\ F_{\mu \nu }^\beta & = H_{\mu \nu }^{\mathrm{core}} + J_{\mu \nu } - F_{\mu \nu }^{\mathrm{XC}\beta } \; . \end{aligned} $   (4.44)

Here, $\ensuremath{\mathbf{F}}^{\mathrm{XC}\alpha }$ and $\ensuremath{\mathbf{F}}^{\mathrm{XC}\beta }$ are the exchange-correlation parts of the Fock matrices and depend on the exchange-correlation functional used. UHF theory is recovered as a special case simply by taking $F_{\mu \nu }^{\mathrm{XC}\alpha } = K_{\mu \nu }^\alpha $, and similarly for $\beta $. Thus, the density and energy are obtained in a manner analogous to that for the HF method. Initial guesses are made for the MO coefficients and an iterative process is applied until self-consistency is achieved.

4.4.3 Overview of Available Functionals

Q-Chem currently has more than 30 exchange functionals as well as more than 30 correlation functionals, and in addition over 150 exchange-correlation (XC) functionals, which refer to functionals that are not separated into exchange and correlation parts, either because the way in which they were parameterized renders such a separation meaningless (e.g., B97-D[28] or $\omega $B97X[29]) or because they are a standard linear combination of exchange and correlation (e.g., PBE[30] or B3LYP[31, 32]). User-defined XC functionals can be created as specified linear combinations of any of the 30+ exchange functionals and/or the 30+ correlation functionals.

KS-DFT functionals can be organized onto a ladder with five rungs, in a classification scheme (“Jacob’s Ladder”) proposed by John Perdew[33] in 2001. The first rung contains a functional that only depends on the (spin-)density $\rho _{\sigma }$, namely, the local spin-density approximation (LSDA). These functionals are exact for the infinite uniform electron gas (UEG), but are highly inaccurate for molecular properties whose densities exhibit significant inhomogeneity. To improve upon the weaknesses of the LSDA, it is necessary to introduce an ingredient that can account for inhomogeneities in the density: the density gradient, $\nabla \rho _{\sigma }$. These generalized gradient approximation (GGA) functionals define the second rung of Jacob’s Ladder and tend to improve significantly upon the LSDA. Two additional ingredients that can be used to further improve the performance of GGA functionals are either the Laplacian of the density $\nabla ^{2}\rho _{\sigma }$, and/or the kinetic energy density,

  \begin{equation} \label{eq:tau} \tau _{\sigma }=\displaystyle \sum _{i}^{n_{\sigma }}\left|\nabla \psi _{i,\sigma }\right|^{2} \; . \end{equation}   (4.47)

While functionals that employ both of these options are available in Q-Chem, the kinetic energy density is by far the more popular ingredient and has been used in many modern functionals to add flexibility to the functional form with respect to both constraint satisfaction (non-empirical functionals) and least-squares fitting (semi-empirical parameterization). Functionals that depend on either of these two ingredients belong to the third rung of the Jacob’s Ladder and are called meta-GGAs. These meta-GGAs often further improve upon GGAs in areas such as thermochemistry, kinetics (reaction barrier heights), and even non-covalent interactions.

Functionals on the fourth rung of Jacob’s Ladder are called hybrid density functionals. This rung contains arguably the most popular density functional of our time, B3LYP, the first functional to see widespread application in chemistry. “Global” hybrid (GH) functionals such as B3LYP (as distinguished from the “range-separated hybrids" introduced below) add a constant fraction of “exact” (Hartree-Fock) exchange to any of the functionals from the first three rungs. Thus, hybrid LSDA, hybrid GGA, and hybrid meta-GGA functionals can be constructed, although the latter two types are much more common. As an example, the formula for the B3LYP functional, as implemented in Q-Chem, is

  \begin{equation} \label{eq:GHGGA} E_{xc}^{\textrm{B3LYP}}=c_{x}E_{x}^{\textrm{HF}}+\left(1-c_{x}-a_{x}\right)E_{x}^{\textrm{Slater}} +a_{x}E_{x}^{\textrm{B88}}+\left(1-a_{c}\right)E_{c}^{\textrm{VWN1RPA}}+a_{c}E_{c}^{\textrm{LYP}} \end{equation}   (4.48)

where $c_{x}=0.20$, $a_{x}=0.72$, and $a_{c}=0.81$.

A more recent approach to introducing exact exchange into the functional form is via range separation. Range-separated hybrid (RSH) functionals split the exact exchange contribution into a short-range (SR) component and a long-range (LR) component, often by means of the error function (erf) and complementary error function ($\rm erfc \equiv 1 - erf$):

  \begin{equation} \label{eq:erf_ partition} \frac{1}{r^{}_{12}} =\frac{\mbox{erfc}(\omega r^{}_{12})}{r^{}_{12}} + \frac{\mbox{erf}(\omega r_{12}^{})}{r^{}_{12}} \; \end{equation}   (4.49)

The first term on the right in Eq. () is singular but short-range, and decays to zero on a length scale of $\sim 1/\omega $, while the second term constitutes a non-singular, long-range background. An RSH XC functional can be expressed generically as

  \begin{equation} \label{eq:RSHGGA} E_{xc}^{\textrm{RSH}}=c^{}_{x,{\rm SR}}E_{x,{\rm SR}}^{\textrm{HF}} + c^{}_{x,\rm LR} E_{x,\rm LR}^{\textrm{HF}} +(1-c^{}_{x,\rm SR}) E_{x,\rm SR}^{\textrm{DFT}} + (1-c^{}_{x,\rm LR})E_{x,\rm LR}^{\textrm{DFT}}+E_{c}^{\textrm{DFT}} \; , \end{equation}   (4.50)

where the SR and LR parts of the Coulomb operator are used, respectively, to evaluate the HF exchange energies $E_{x,\rm SR}^{\textrm{HF}}$ and $E_{x,\rm LR}^{\textrm{HF}}$. The corresponding DFT exchange functional is partitioned in the same manner, but the correlation energy $E_ c^{\textrm{DFT}}$ is evaluated using the full Coulomb operator, $r_{12}^{-1}$. Of the two linear parameters in Eq. eq:RSHGGA, $c^{}_{x,\rm LR}$ is usually either set to 1 to define long-range corrected (LRC) RSH functionals (see Section 4.4.6) or else set to 0, which defines screened-exchange (SE) RSH functionals. On the other hand, the fraction of short-range exact exchange ($c^{}_{x,\rm SR}$) can either be determined via least-squares fitting, theoretically justified using the adiabatic connection, or simply set to zero. As with the global hybrids, RSH functionals can be fashioned using all of the ingredients from the lower three rungs. The rate at which the local DFT exchange is turned off and the non-local exact exchange is turned on is controlled by the parameter $\omega $. Large values of $\omega $ tend to lead to attenuators that are less smooth (unless the fraction of short-range exact exchange is very large), while small values of (e.g., $\omega = $0.2–0.3 bohr$^{-1}$) are the most common in semi-empirical RSH functionals.

The final rung on Jacob’s Ladder contains functionals that utilize not only occupied orbitals (via exact exchange), but virtual orbitals as well (via methods such as MP2 or the random phase approximation, RPA). These double hybrids (DH) are the most expensive density functionals available in Q-Chem, but can also be very accurate. The most basic form of a DH functional is

  \begin{equation} \label{eq:DHGGA} E_{xc}^{\textrm{DH}}=c_{x}E_{x}^{\textrm{HF}}+\left(1-c_{x}\right)E_{x}^{\textrm{DFT}} +c_{c}E_{x}^{\textrm{MP2}}+\left(1-c_{c}\right)E_{c}^{\textrm{DFT}} \; . \end{equation}   (4.51)

As with hybrids, the coefficients can either be theoretically motivated or empirically determined. In addition, double hybrids can utilize exact exchange both globally or via range-separation, and their components can be as primitive as LSDA or as advanced as in meta-GGA functionals. More information on double hybrids can be found in Section 4.4.8.

Finally, the last major advance in KS-DFT in recent years has been the development of methods that are capable of accurately describing non-covalent interactions, particularly dispersion. All of the functionals from Jacob’s Ladder can technically be combined with these dispersion corrections, although in some cases the combination is detrimental, particularly for semi-empirical functionals that were parameterized in part using data sets of non-covalent interactions, and already tend to overestimate non-covalent interaction energies. The most popular such methods available in Q-Chem are:

Below, we categorize the functionals that are available in Q-Chem, including exchange functionals (Section 4.4.3.2), correlation functionals (Section 4.4.3.3), and exchange-correlation functionals (Section  4.4.3.4). Within each category the functionals will be categorized according to Jacob’s Ladder. Exchange and correlation functionals can be invoked using the $rem variables EXCHANGE and CORRELATION, while the exchange-correlation functionals can be invoked either by setting the $rem variable METHOD or alternatively (in most cases, and for backwards compatibility with earlier versions of Q-Chem) by using the $rem variable EXCHANGE. Some caution is warranted here. While setting METHOD to PBE, for example, requests the Perdew-Burke-Ernzerhof (PBE) exchange-correlation functional,[30] which includes both PBE exchange and PBE correlation, setting EXCHANGE = PBE requests only the exchange component and setting CORRELATION = PBE requests only the correlation component. Setting both of these values is equivalent to specifying METHOD = PBE.

 

Single-Point

Optimization

Frequency

Ground State

LSDA$^{\dagger \star }$

LSDA$^{\dagger \star }$

LSDA$^{\star }$

 

GGA$^{\dagger \star }$

GGA$^{\dagger \star }$

GGA$^{\star }$

 

meta-GGA$^{\dagger \star }$

meta-GGA$^{\dagger }$

—

 

GH$^{\dagger \star }$

GH$^{\dagger \star }$

GH$^{\star }$

 

RSH$^{\dagger \star }$

RSH$^{\dagger \star }$

RSH$^{\star }$

 

NLC$^{\dagger \star }$

NLC$^{\dagger \star }$

—

 

DFT-D

DFT-D

DFT-D

 

XDM

—

—

TDDFT

LSDA$^{\dagger \star }$

LSDA$^{\dagger \star }$

LSDA

 

GGA$^{\dagger \star }$

GGA$^{\dagger \star }$

GGA

 

meta-GGA$^{\dagger \star }$

—

—

 

GH$^{\dagger \star }$

GH$^{\dagger \star }$

GH

 

RSH$^{\dagger \star }$

RSH$^{\dagger \star }$

—

 

—

—

—

 

DFT-D

DFT-D

DFT-D

 

—

—

—

$^\dagger $OpenMP parallelization available

$^\star $MPI parallelization available
Table 4.2: Available analytic properties and parallelization for SCF calculations.

Finally, Table 4.2 provides a summary, arranged according to Jacob’s Ladder, of which categories of functionals are available with analytic first derivatives (for geometry optimizations) or second derivatives (for vibrational frequency calculations). If analytic derivatives are not available for the requested job type, Q-Chem will automatically generate them via finite difference. Tests of the finite-difference procedure, in cases where analytic second derivatives are  available, suggest that finite-difference frequencies are accurate to $< 1$ cm$^{-1}$, except for very low-frequency, non-bonded modes [44]. Also listed in Table 4.2 are which functionals are available for excited-state time-dependent DFT (TDDFT) calculations, as described in Section 6.3. Lastly, Table 4.2 describes which functionals have been parallelized with OpenMP and/or MPI.

4.4.3.1 Suggested Density Functionals

Q-Chem contains over 150 exchange-correlation functionals, not counting those that can be straightforwardly appended with a dispersion correction (such as B3LYP-D3). Therefore, we suggest a few functionals from the second through fourth rungs of Jacob’s Ladder in order to guide functional selection. Most of these suggestions come from a benchmark of over 200 density functionals on a vast database of nearly 5000 data points, covering non-covalent interactions, isomerization energies, thermochemistry, and barrier heights. The single recommended method from each category is indicated in bold.

From the GGAs on Rung 2, we recommend:

From the meta-GGAs on Rung 3, we recommend:

From the hybrid GGAs on Rung 4, we recommend:

From the hybrid meta-GGAs on Rung 4, we recommend:

4.4.3.2 Exchange Functionals

Note: All exchange functionals in this section can be invoked using the $rem variable EXCHANGE. Popular and/or recommended functionals within each class are listed first and indicated in bold. The rest are in alphabetical order.

Local Spin-Density Approximation (LSDA)

Generalized Gradient Approximation (GGA)

Meta-Generalized Gradient Approximation (meta-GGA)

4.4.3.3 Correlation Functionals

Note: All correlation functionals in this section can be invoked using the $rem variable CORRELATION. Popular and/or recommended functionals within each class are listed first and indicated in bold. The rest are in alphabetical order.

Local Spin-Density Approximation (LSDA)

Generalized Gradient Approximation (GGA)

Meta-Generalized Gradient Approximation (meta-GGA)

4.4.3.4 Exchange-Correlation Functionals

Note: All exchange-correlation functionals in this section can be invoked using the $rem variable METHOD. For backwards compatibility, all of the exchange-correlation functionals except for the ones marked with an asterisk can be utilized with the $rem variable EXCHANGE. Popular and/or recommended functionals within each class are listed first and indicated in bold. The rest are in alphabetical order.

Local Spin-Density Approximation (LSDA)

Generalized Gradient Approximation (GGA)

Meta-Generalized Gradient Approximation (meta-GGA)

Global Hybrid Generalized Gradient Approximation (GH GGA)

Global Hybrid Meta-Generalized Gradient Approximation (GH meta-GGA)

Range-Separated Hybrid Generalized Gradient Approximation (RSH GGA)

Range-Separated Hybrid Meta-Generalized Gradient Approximation (RSH meta-GGA)

Double Hybrid Generalized Gradient Approximation (DH GGA)

Note: In order to use the resolution-of-the-identity approximation for the MP2 component, specify an auxiliary basis set with the $rem variable AUX_BASIS

Double Hybrid Meta-Generalized Gradient Approximation (DH MGGA)

4.4.3.5 Specialized Functionals

4.4.3.6 User-Defined Density Functionals

Users can also request a customized density functional consisting of any linear combination of exchange and/or correlation functionals available in Q-Chem. A “general” density functional of this sort is requested by setting EXCHANGE = GEN and then specifying the functional by means of an $xc_functional input section consisting of one line for each desired exchange (X) or correlation (C) component of the functional, and having the format shown below. 

$xc_functional
   X   exchange_symbol   coefficient
   X   exchange_symbol   coefficient
   ...
   C   correlation_symbol   coefficient
   C   correlation_symbol   coefficient
   ...
   K   coefficient
$end

Each line requires three variables: X or C to designate whether this is an exchange or correlation component; the symbolic representation of the functional, as would be used for the EXCHANGE or CORRELATION keywords variables as described above; and a real number coefficient for each component. Note that Hartree-Fock exchange can be designated either as “X" or as “K". Examples are shown below.

Example 4.21  Q-Chem input for H$_2$O with the B3tLap functional.

$molecule
  0  1
  O
  H1  O  oh
  H2  O  oh  H1  hoh
  oh  =   0.97
  hoh = 120.0
$end

$rem
   EXCHANGE     gen 
   CORRELATION  none 
   BASIS        g3large   ! recommended for high accuracy
   THRESH       14        ! and better convergence
$end

$xc_functional
X    Becke     0.726
X    S         0.0966
C    PK06      1.0
K    0.1713
$end

Example 4.22  Q-Chem input for H$_2$O with the BR89B94hyb functional.

$molecule
   0  1
   O
   H1  O  oh
   H2  O  oh  H1  hoh
   oh  =   0.97
   hoh = 120.0
$end

$rem
   EXCHANGE     gen
   CORRELATION  none
   BASIS        g3large   ! recommended for high accuracy
   THRESH       14        ! and better convergence
$end

$xc_functional
   X    BR89      0.846
   C    B94hyb    1.0
   K              0.154
$end

The next two examples illustrate the use of the RI-B05 and RI-PSTS functionals. These are presently available only for single-point calculations, and convergence is greatly facilitated by obtaining converged SCF orbitals from, e.g., an LDA or HF calculation first. (LDA is used in the example below but HF can be substituted.) Use of the RI approximation (Section 5.6) requires specification of an auxiliary basis set.

Example 4.23  Q-Chem input of H$_2$ using RI-B05.

$comment
   H2, example of SP RI-B05.  First do a well-converged LSD, G3LARGE is the
   basis of choice for good accuracy. The input lines
   PURECART    2222
   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
   SCF_GUESS        core
   METHOD           lda
   BASIS            g3large
   PURECART         2222
   THRESH           14
   MAX_SCF_CYCLES   80
   PRINT_INPUT      true
   INCDFT           false
   SYM_IGNORE       true
   SYMMETRY         false
   SCF_CONVERGENCE  9      
$end

@@@@

$comment
   For the time being the following input lines are obligatory:
   PURECART        2222
   AUX_BASIS       riB05-cc-pvtz
   DFT_CUTOFFS     0
   MAX_SCF_CYCLES  0
$end

$molecule
   read
$end

$rem
   SCF_GUESS       read
   EXCHANGE        b05   ! or set to psts for ri-psts
   PURECART        2222
   BASIS           g3large
   AUX_BASIS       rib05-cc-pvtz  ! the aux basis for both RI-B05 and RI-PSTS
   THRESH          4
   PRINT_INPUT     true
   INCDFT          false
   SYM_IGNORE      true
   SYMMETRY           false
   MAX_SCF_CYCLES  0
   DFT_CUTOFFS     0
$end

4.4.4 Basic DFT Job Control

Basic SCF job control was described in Section 4.3 in the context of Hartree-Fock theory and is largely the same for DFT. The keywords METHOD and BASIS are required, although for DFT the former could be substituted by specifying EXCHANGE and CORRELATION instead.

METHOD

Specifies the exchange-correlation functional.


TYPE:

STRING


DEFAULT:

No default


OPTIONS:

NAME

Use METHOD = NAME, where NAME is either HF for Hartree-Fock theory or

 

else one of the DFT methods listed in Section 4.4.3.4.


RECOMMENDATION:

In general, consult the literature to guide your selection. Our recommendations for DFT are indicated in bold in Section 4.4.3.4.


EXCHANGE

Specifies the exchange functional (or most exchange-correlation functionals for backwards compatibility).


TYPE:

STRING


DEFAULT:

No default


OPTIONS:

NAME

Use EXCHANGE = NAME, where NAME is either:

 

1) One of the exchange functionals listed in Section 4.4.3.2

 

2) One of the XC functionals listed in Section 4.4.3.4 that is not marked with an

 

asterisk.

 

3) GEN, for a user-defined functional (see Section 4.4.3.6).


RECOMMENDATION:

In general, consult the literature to guide your selection. Our recommendations are indicated in bold in Sections 4.4.3.4 and  4.4.3.2.


CORRELATION

Specifies the correlation functional.


TYPE:

STRING


DEFAULT:

NONE


OPTIONS:

NAME

Use CORRELATION = NAME, where NAME is one of the correlation functionals

 

listed in Section 4.4.3.3.


RECOMMENDATION:

In general, consult the literature to guide your selection. Our recommendations are indicated in bold in Section 4.4.3.3.


The following $rem variables are related to the choice of the quadrature grid required to integrate the XC part of the functional, which does not appear in Hartree-Fock theory. DFT quadrature grids are described in Section 4.4.5.

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:

Functional-dependent; see Table 4.4.


OPTIONS:

0

Use SG-0 for H, C, N, and O; SG-1 for all other atoms.

$n$

Use SG-$n$ for all atoms, $n = 1, 2$, or 3

$XY$

A string of two six-digit integers $X$ and $Y$, where $X$ is the number of radial points

 

and $Y$ is the number of angular points where possible numbers of Lebedev angular

 

points, which must be an allowed value from Table 4.3 in Section 4.4.5.

$-XY$

Similar format for Gauss-Legendre grids, with the six-digit integer $X$ corresponding

 

to the number of radial points and the six-digit integer $Y$ providing the number of

 

Gauss-Legendre angular points, $Y = 2N^2$.


RECOMMENDATION:

Use the default unless numerical integration problems arise. Larger grids may be required for optimization and frequency calculations.


NL_GRID

Specifies the grid to use for non-local correlation.


TYPE:

INTEGER


DEFAULT:

1


OPTIONS:

Same as for XC_GRID


RECOMMENDATION:

Use the default unless computational cost becomes prohibitive, in which case SG-0 may be used. XC_GRID should generally be finer than NL_GRID.


XC_SMART_GRID

Uses SG-0 (where available) for early SCF cycles, and switches to the (larger) target grid specified by XC_GRID for final cycles of the SCF.


TYPE:

LOGICAL


DEFAULT:

FALSE


OPTIONS:

TRUE (or 1)

Use the smaller grid for the initial cycles.

FALSE (or 0)

Use the target grid for all SCF cycles.


RECOMMENDATION:

The use of the smart grid can save some time on initial SCF cycles.


4.4.5 DFT Numerical Quadrature

In practical DFT calculations, the forms of the approximate exchange-correlation functionals used are quite complicated, such that the required integrals involving the functionals generally cannot be evaluated analytically. Q-Chem evaluates these integrals through numerical quadrature directly applied to the exchange-correlation integrand. Several standard quadrature grids are available (“SG-$n$”, $n = 0, 1, 2, 3$), with a default value that is automatically set according to the complexity of the functional in question.

The quadrature approach in Q-Chem is generally similar to that found in many DFT programs. The multi-center XC integrals are first partitioned into “atomic” contributions using a nuclear weight function. Q-Chem uses the nuclear partitioning of Becke [176], though without the “atomic size adjustments” of Ref. Becke:1988a. The atomic integrals are then evaluated through standard one-center numerical techniques. Thus, the exchange-correlation energy is obtained as

  \begin{equation} \label{eq439} E_{\mathrm{XC}} =\sum _ A^{\rm atoms} \sum _{i\in A}^{\rm points} w_{Ai}^{} f(\mathbf{r}_{Ai}^{}) \; , \end{equation}   (4.52)

where the function $f$ is the aforementioned XC integrand and the quantities $w_{Ai}$ are the quadrature weights. The sum over $i$ runs over grid points belonging to atom $A$, which are located at positions $\mathbf{r}^{}_{Ai} = \mathbf{R}_ A^{} + \mathbf{r}_ i^{}$, so this approach requires only the choice of a suitable one-center integration grid (to define the $\mathbf{r}_ i^{}$), which is independent of nuclear configuration. These grids are implemented in Q-Chem in a way that ensures that the $E_{\rm XC}$ is rotationally-invariant, i.e., that is does not change when the molecule undergoes rigid rotation in space [177].

Quadrature grids are further separated into radial and angular parts. Within Q-Chem, the radial part is usually treated by the Euler-Maclaurin scheme proposed by Murray et al. [178], which maps the semi-infinite domain $[0,\infty )$ onto $[0,1)$ and applies the extended trapezoid rule to the transformed integrand. Alternatively, Gill and Chien [179] proposed a radial scheme based on a Gaussian quadrature on the interval $[0,1]$ with a different weight function. This “MultiExp" radial quadrature 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. However, the task of generating the MultiExp quadrature points becomes increasingly ill-conditioned as the number of radial points increases, so that a “double exponential" radial quadrature [180, 181] is used for the largest standard grids in Q-Chem, namely SG-2 and SG-3 [182]. (See Section 4.4.5.2.)

4.4.5.1 Angular Grids

For a fixed value of the radial spherical-polar coordinate $r$, a function $f(\mathbf{r}) \equiv f(r,\theta ,\phi )$ has an exact expansion in spherical harmonic functions,

  \begin{equation} f(r,\theta ,\phi ) = \sum _{\ell =0}^\infty \sum _{m=-\ell }^\ell c_{\ell m}^{} \mbox{Y}_{\ell m}(\theta ,\phi ) \; . \end{equation}   (4.53)

Angular quadrature grids are designed to integrate $f(r,\theta ,\phi )$ for fixed $r$, and are often characterized by their degree, meaning the maximum value of $\ell $ for which the quadrature is exact, as well as by their efficiency, meaning the number of spherical harmonics exactly integrated per degree of freedom in the formula. Q-Chem supports the following two types of angular grids.

Combining these radial and angular schemes yields an intimidating selection of quadratures, so it is useful to standardize the grids. This is done for the convenience of the user, to facilitate comparisons in the literature, and also for developers wishing to compared detailed results between different software programs, because the total electronic energy is sensitive to the details of the grid, just as it is sensitive to details of the basis set. Standard quadrature grids are discussed next.

4.4.5.2 Standard Quadrature Grids, SG-$\bm@general \boldmath \m@ne \mv@bold \bm@command {n}$

Four different “standard grids" are available in Q-Chem, designated SG-$n$, for $n = 0, 1, 2$, or 3; both quality and the computational cost of these grids increases with $n$. These grids are constructed starting from a “parent” grid ($N_ r,~ N_\Omega $) consisting of $N_ r$ radial spheres with $N_\Omega $ angular (Lebedev) grid points on each, then systematically pruning 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 and at long distance from the nuclei where it varies slowly. A large number of points is retained in the valence region where angular accuracy is critical. The SG-$n$ grids are summarized in Table 4.4. While many electronic structure programs use some kind of procedure to delete unnecessary grid points in the interest of computational efficiency, Q-Chem’s SG-$n$ grids are notable in that the complete grid specifications are available in the peer-reviewed literature [187, 188, 182], to facilitate reproduction of Q-Chem DFT calculations using other electronic structure programs.

Pruned

Ref.

Parent Grid

No. Grid Points

Default Grid for

Grid

 

($N_ r,N_\Omega $)

(C atom)$^ a$

Which Functionals?$^ b$

SG-0

Chien:2006

(23, 170)

1,390 (36%)

None

SG-1

Gill:1993

(50, 194)

3,816 (39%)

LDA, most GGAs and hybrids

SG-2

Dasgupta:2017

(75, 302)

7,790 (34%)

Meta-GGAs; B95- and B97-based functionals

SG-3

Dasgupta:2017

(99, 590)

17,674 (30%)

Minnesota functionals

$^ a$Number in parenthesis is the fraction of points retained from the parent grid

$^ b$Reflects Q-Chem versions since v. 4.4.2

Table 4.4: Standard quadrature grids available in Q-Chem, along with the number of grid points for a carbon atom, showing the reduction in grid points due to pruning.

The SG-0 and SG-1 grids are designed for calculations on large molecules using GGA functionals. SG-1 affords integration errors on the order of $\sim $0.2 kcal/mol for medium-sized molecules and GGA functionals, including for demanding test cases such as reaction enthalpies for isomerizations. (Integration errors in total energies are no more than a few $\mu $hartree, or $\sim $0.01 kcal/mol.) The SG-0 grid was derived in similar fashion, and affords a root-mean-square error in atomization energies of 72 $\mu $hartree with respect to SG-1, while relative energies are reproduced well [188]. In either case, errors of this magnitude are typically considerably smaller than the intrinsic errors in GGA energies, and hence acceptable. As seen in Table 4.4, SG-1 retains $<40$% of the grid points of its parent grid, which translates directly into cost savings.

Both SG-0 and SG-1 were optimized so that the integration error in the energy falls below a target threshold, but derivatives of the energy (including such properties as (hyper)polarizabilities [189]) are often more sensitive to the quality of the integration grid. Special care is required, for example, when imaginary vibrational frequencies are encountered, as low-frequency (but real) vibrational frequencies can manifest as imaginary if the grid is sparse. If imaginary frequencies are found, or if there is some doubt about the frequencies reported by Q-Chem, the recommended procedure is to perform the geometry optimization and vibrational frequency calculations again using a higher-quality grid. (The optimization should converge quite quickly if the previously-optimized geometry is used as an initial guess.)

SG-1 was the default DFT integration grid for all density functionals for Q-Chem versions 3.2–4.4. (Note that SG-0 was re-optimized for Q-Chem v. 3.0, so results may differ slightly as compared to older versions of the program.) Beginning with Q-Chem v. 4.4.2, however, the default grid is functional-dependent, as summarized in Table 4.4. This is a reflection of the fact that although SG-1 is adequate for energy calculations using most GGA and hybrid functionals (although care must be taken for some other properties, as discussed below), it is not  adequate to integrate many functionals developed since $\sim $2005. These include meta-GGAs, which are more complicated due to their dependence on the kinetic energy density ($\tau _\sigma $ in Eq. eq:tau) and/or the Laplacian of the density ($\nabla ^2\rho _\sigma $). Functionals based on B97, along with the Minnesota suite of functionals [124, 190], contain relatively complicated expressions for the exchange inhomogeneity factor, and are therefore also more sensitive to the quality of the integration grid [191, 139, 182]. To integrate these modern density functionals, the SG-2 and SG-3 grids were developed [182], which are pruned versions of the medium-quality (75, 302) and high-quality (99, 590) integration grids, respectively. Tests of properties known to be highly sensitive to the quality of the integration grid, such as vibrational frequencies, hyper-polarizabilities, and potential energy curves for non-bonded interactions, demonstrate that SG-2 is usually adequate for meta-GGAs and B97-based functionals, and in many cases is essentially converged with respect to an unpruned (250, 974) grid [182]. The Minnesota functionals are more sensitive to the grid, and while SG-3 is often adequate, it is not completely converged in the case of non-bonded interactions. [182].

4.4.5.3 Consistency Check and Cutoffs

Whenever Q-Chem calculates numerical density functional integrals, the electron density itself is also integrated numerically as a test of the quality of the numerical quadrature. The extent to which this numerical result differs from the number of electrons is an indication of the accuracy of the other numerical integrals. A warning message is printed whenever the relative error in the numerical electron count reaches 0.01%, indicating that the numerical XC results may not be reliable. If the warning appears on the first SCF cycle it is probably not serious, because the initial-guess density matrix is sometimes not idempotent. This is the case with the SAD guess discussed in Section 4.5, and also with a density matrix that is taken from a previous geometry optimization cycle, and in such cases the problem will likely correct itself in subsequent SCF iterations. If the warning persists, however, then one should consider either using a finer grid or else selecting an alternative initial guess.

By default, Q-Chem will estimate the magnitude of various XC contributions on the grid and eliminate those determined to be numerically insignificant. Q-Chem uses specially-developed cutoff procedures which permits evaluation of the XC energy and potential in only ${\cal {O}}({N})$ work for large molecules. This is a significant improvement over the formal ${\cal {O}}({N^3})$ scaling of the XC cost, and is critical in enabling DFT calculations to be carried out on very large systems. In 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. A remedy may be to increase the size of the quadrature grid.

4.4.6 Range-Separated Hybrid Density Functionals

Whereas RSH functionals such as LRC-$\omega $PBE are attempts to add 100% LR Hartree-Fock exchange with minimal perturbation to the original functional (PBE, in this example), other RSH functionals are of a more empirical nature and their range-separation parameters have been carefully parameterized along with all of the other parameters in the functional. These cases are functionals are discussed first, in Section 4.4.6.1, because their range-separation parameters should be taken as fixed. User-defined values of the range-separation parameter are discussed in Section 4.4.6.2, and Section 4.4.6.3 discusses a procedure for which an optimal, system-specific value of this parameter ($\omega $ or $\mu $) can be chosen for functionals such as LRC-$\omega $PBE or LRC-$\mu $PBE.

4.4.6.1 Semi-Empirical RSH Functionals

Semi-empirical RSH functionals for which the range-separation parameter should be considered fixed include the $\omega $B97, $\omega $B97X, and $\omega $B97X-D functionals developed by Chai and Head-Gordon [29, 39]; $\omega $B97X-V and $\omega $B97M-V from Mardirossian and Head-Gordon [139, 151]; M11 from Peverati and Truhlar [152]; $\omega $B97X-D3, $\omega $M05-D, and $\omega $M06-D3 from Chai and coworkers [153, 140]; and the screened exchange functionals N12-SX and MN12-SX from Truhlar and co-workers [149]. More recently, Mardirossian and Head-Gordon developed two RSH functionals, $\omega $B97X-V and $\omega $B97M-V, via a combinatorial approach by screening over 100,000 possible functionals in the first case and over 10 billion possible functionals in the second case. Both of the latter functionals use the VV10 non-local correlation functional in order to improve the description of non-covalent interactions and isomerization energies. $\omega $B97M-V is a 12-parameter meta-GGA with 15% short-range exact exchange and 100% long-range exact exchange and is one of the most accurate functionals available through rung 4 of Jacob’s Ladder, across a wide variety of applications. This has been verified by benchmarking the functional on nearly 5000 data points against over 200 alternative functionals available in Q-Chem [151].

4.4.6.2 User-Defined RSH Functionals

As pointed out in Ref. Dreuw:2003 and elsewhere, the description of charge-transfer excited states within density functional theory (or more precisely, time-dependent DFT, which is discussed in Section 6.3) requires full (100%) non-local HF exchange, at least in the limit of large donor–acceptor distance. Hybrid functionals such as B3LYP [31, 32] and PBE0 [193] that are well-established and in widespread use, however, employ only 20% and 25% HF exchange, respectively. While these functionals provide excellent results for many ground-state properties, they cannot correctly describe the distance dependence of charge-transfer 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 non-covalent clusters, where TDDFT often predicts a near-continuum of spurious, low-lying charge transfer states [194, 195]. The problems with TDDFT’s description of charge transfer are not limited to large donor–acceptor distances, but have been observed at $\sim $2  separation, in systems as small as uracil–(H$_2$O)$_4$ [194]. 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% HF exchange, though few such functionals exist to date. An alternative option is to attempt to preserve the form of common GGAs and hybrid functionals at short range (i.e., keep the 25% HF exchange in PBE0) while incorporating 100% HF exchange at long range, which provides a rigorously correct description of the long-range distance dependence of charge-transfer excitation energies, but aims to avoid contaminating short-range exchange-correlation effects with additional HF exchange. The separation is accomplished using the range-separation ansatz that was introduced in Section 4.4.3. In particular, functionals that use 100% HF exchange at long range ($c^{}_{x,\rm LR}=1$ in Eq. eq:RSHGGA) are known as “long-range-corrected” (LRC) functionals. An LRC version of PBE0 would, for example, have $c^{}_{x,\rm SR} = 0.25$.

To fully specify an LRC functional, one must choose a value for the range separation parameter $\omega $ in Eq. (). In the limit $\omega \rightarrow 0$, the LRC functional in Eq. eq:RSHGGA reduces to a non-RSH functional where there is no “SR” or “LR”, because all exchange and correlation energies are evaluated using the full Coulomb operator, $r_{12}^{-1}$. Meanwhile the $\omega \rightarrow \infty $ limit corresponds to a new functional, $E_{xc}^{\rm RSH} = E_{c}^{} + E_ x^{\rm HF}$. Full HF exchange is inappropriate for use with most contemporary GGA correlation functionals, so the latter limit is expected to perform quite poorly. Values of $\omega > 1.0$ bohr$^{-1}$ are likely not worth considering, according to benchmark tests [196, 147].

Evaluation of the short- and long-range HF exchange energies is straightforward [197], so the crux of any RSH functional is the form of the short-range GGA exchange functional, and several such functionals are available in Q-Chem. These include short-range variants of the B88 and PBE exchange described by Hirao and co-workers [53, 146], called $\mu $B88 and $\mu $PBE in Q-Chem [198], and an alternative formulation of short-range PBE exchange proposed by Scuseria and co-workers [62], which is known as $\omega $PBE. These functionals are available in Q-Chem thanks to the efforts of the Herbert group [147, 148]. By way of notation, the terms “$\mu $PBE”, “$\omega $PBE”, etc., refer only to the short-range exchange functional, $E_{x,\rm SR}^{\rm DFT}$ in Eq. eq:RSHGGA. These functionals could be used in “screened exchange” mode, as described in Section 4.4.3, as for example in the HSE03 functional [199], therefore the designation “LRC-$\omega $PBE”, for example, should only be used when the short-range exchange functional $\omega $PBE is combined with 100% Hartree-Fock exchange in the long range.

In general, LRC-DFT functionals have been shown to remove the near-continuum of spurious charge-transfer excited states that appear in large-scale TDDFT calculations [196]. However, certain results depend sensitively upon the value of the range-separation parameter $\omega $ [196, 147, 148, 195, 200], especially in TDDFT calculations (Section 6.3) and therefore the results of LRC-DFT calculations must therefore be interpreted with caution, and probably for a range of $\omega $ values. This can be accomplished by requesting a functional that contains some short-range GGA exchange functional ($\omega $PBE or $\mu $PBE, in the examples mentioned above), in combination with setting the $rem variable LRC_DFT = TRUE, which requests the addition of 100% Hartree-Fock exchange in the long-range. Basic job-control variables and an example can be found below. The value of the range-separation parameter is then controlled by the variable OMEGA, as shown in the examples below.

LRC_DFT

Controls the application of long-range-corrected DFT


TYPE:

LOGICAL


DEFAULT:

FALSE


OPTIONS:

FALSE (or 0)

Do not apply long-range correction.

TRUE (or 1)

Add 100% long-range Hartree-Fock exchange to the requested functional.


RECOMMENDATION:

The $rem variable OMEGA must also be specified, in order to set the range-separation parameter.


OMEGA

Sets the range-separation parameter, $\omega $, also known as $\mu $, in functionals based on Hirao’s RSH scheme.


TYPE:

INTEGER


DEFAULT:

No default


OPTIONS:

$n$

Corresponding to $\omega = n/1000$, in units of bohr$^{-1}$


RECOMMENDATION:

None


COMBINE_K

Controls separate or combined builds for short-range and long-range K


TYPE:

LOGICAL


DEFAULT:

FALSE


OPTIONS:

FALSE (or 0)

Build short-range and long-range K separately (twice as expensive as a global hybrid)

TRUE (or 1)

Build short-range and long-range K together ($\approx $ as expensive as a global hybrid)


RECOMMENDATION:

Most pre-defined range-separated hybrid functionals in Q-Chem use this feature by default. However, if a user-specified RSH is desired, it is necessary to manually turn this feature on.


Example 4.24  Application of LRC-BOP to $\rm (H_2O)_2^-$.

$comment
The value of omega is 0.47 by default but can
be overwritten by specifying OMEGA.
$end

$molecule
-1 2
O           1.347338    -0.017773    -0.071860
H           1.824285     0.813088     0.117645
H           1.805176    -0.695567     0.461913
O          -1.523051    -0.002159    -0.090765
H          -0.544777    -0.024370    -0.165445
H          -1.682218     0.174228     0.849364
$end

$rem
   EXCHANGE      LRC-BOP    
   BASIS         6-31(1+,3+)G*
   LRC_DFT       TRUE
   OMEGA         300      ! = 0.300 bohr**(-1)
$end

Recently, Rohrdanz et al. [148] have published a thorough benchmark study of both ground- and excited-state properties, using the “LRC-$\omega $PBEh” functional, in which the “h” indicates a short-range hybrid (i.e., the presence of some short-range HF exchange). Empirically-optimized parameters of $c^{}_{x,\rm SR} =0.2$ (see Eq. eq:RSHGGA) and $\omega = 0.2$ bohr$^{-1}$ were obtained [148], and these parameters are taken as the defaults for LRC-$\omega $PBEh. Caution is warranted, however, especially in TDDFT calculations for large systems, as excitation energies for states that exhibit charge-transfer character can be rather sensitive to the precise value of $\omega $.[195, 148] In such cases (and maybe in general), the “tuning” procedure described in Section 4.4.6.3 is recommended.

Example 4.25  Application of LRC-$\omega $PBEh to the $\rm C_2H_4$–$\rm C_2F_4$ dimer at 5  separation.

$comment
    This example uses the "optimal" parameter set discussed above.  
    It can also be run by setting METHOD = LRC-wPBEh.
$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         6-31+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 also the asymptotic corrections that will be discussed in Section 4.4.9.1 are thought to reduce self-interaction error in approximate DFT. 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, $N$, because $E(N)$ should in principle consist of a sequence of line segments with abrupt changes in slope (the so-called derivative discontinuity [201, 202]) at integer values of $N$, but in practice these $E(N)$ plots bow away from straight-line segments [201]. Examination of such plots has been suggested as a means to adjust the fraction of short-range exchange in an LRC functional [203], while the range-separation parameter is tuned as described in Section 4.4.6.3.

Example 4.26  Example of a DFT job with a fractional number of electrons. Here, we make the $-1.x$ anion of fluoride by subtracting a fraction of an electron from the HOMO of F$^{2-}$.

$comment
   Subtracting a whole electron recovers the energy of F-.
   Adding electrons to the LUMO is possible as well.
$end

$rem
   EXCHANGE              b3lyp
   BASIS                 6-31+G*
   FRACTIONAL_ELECTRON  -500   ! divide by 1000 to get the fraction, -0.5 here.
$end

$molecule
-2 2
F
$end

4.4.6.3 Tuned RSH Functionals

Whereas the range-separation parameters for the functionals described in Section 4.4.6.1 are wholly empirical in nature and should not be adjusted, for the functionals described in Section 4.4.6.2 some adjustment was suggested, especially for TDDFT calculations and for any properties that require interpretation of orbital energies, e.g., HOMO/LUMO gaps. This adjustment can be performed in a non-empirical, albeit system-specific way, by “tuning” the value of $\omega $ in order to satisfy certain criteria that ought rigorously to be satisfied by an exact density functional.

System-specific optimization of $\omega $ is based on the Koopmans condition that would be satisfied for the exact density functional, namely, that for an $N$-electron molecule [204]

  \begin{equation} \label{eq:IP-tuning} -\varepsilon ^{}_{\rm HOMO}(N) = \mbox{IE}(N) \equiv E(N-1) - E(N) \; . \end{equation}   (4.54)

In other words, the HOMO eigenvalue should be equal to minus the ionization energy (IE), where the latter is defined by a $\Delta $SCF procedure [205, 206]. When an RSH functional is used, all of the quantities in Eq. eq:IP-tuning are $\omega $-dependent, so this parameter is adjusted until the condition in Eq. eq:IP-tuning is met, which requires a sequence of SCF calculations on both the neutral and ionized species, using different values of $\omega $. For proper description of charge-transfer states, Baer and co-workers [204] suggest finding the value of $\omega $ that (to the extent possible) satisfies both Eq. eq:IP-tuning for the neutral donor molecule and its analogue for the anion of the acceptor species. Note that for a given approximate density functional, there is no guarantee that the IE condition can actually be satisfied for any value of $\omega $, but in practice it usually can be, and published benchmarks suggest that this system-specific approach affords the most accurate values of IEs and TDDFT excitation energies [207, 208, 204]. It should be noted, however, that the optimal value of $\omega $ can very dramatically with system size, e.g., it is very different for the cluster anion (H$_2$O)$^-_{6}$ than it is for cluster (H$_2$O)$^-_{70}$ [200].

A script that optimizes $\omega $, called OptOmegaIPEA.pl, is located in the $QC/bin directory. The script scans $\omega $ over the range 0.1–0.8 bohr$^{-1}$, corresponding to values of the $rem variable OMEGA in the range 100–800. See the script for the instructions how to modify the script to scan over a wider 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), its anion (M.in), and its cation (P.in), and then run the command 

OptOmegaIPEA.pl >& optomega

which both generates the input files (N_*, P_*, M_*) and runs Q-Chem on them, writing the optimization output into optomega. This script applies the IE condition to both the neutral molecule and its anion, minimizing the sum of $(\mbox{IE}+\varepsilon ^{}_{\rm HOMO})^2$ for these two species. A similar script, OptOmegaIP.pl, uses Eq. eq:IP-tuning for the neutral molecule only.

Note:  (i) If the system does not have positive EA, then the tuning should be done according to the IP condition only. The IP/EA script will yield an incorrect value of $\omega $ 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.

Although the tuning procedure was originally developed by Baer and co-workers using the BNL functional [207, 208, 204], it has more recently been applied using functionals such as LRC-$\omega $PBE (see, e.g., Ref. Uhlig:2014), and the scripts will work with functionals other than BNL.

4.4.6.4 Tuned RSH Functionals Based on the Global Density-Dependent Condition

The value of range-separation parameter based on IP tuning procedure ($\omega ^{}_{\rm IP}$) exhibits a troublesome dependence on system size [200, 209, 210, 211, 212]. An alternative method to select $\omega $ is the global density-dependent (GDD) tuning procedure[213], in which the optimal value

  \begin{equation} \label{eq:wgdd} \omega ^{}_{\rm GDD} = C \langle d^2_{\rm x} \rangle ^{-1/2} \end{equation}   (4.55)

is related to the average distance $d_{\rm x}$ between an electron in the outer regions of a molecule and the exchange hole in the region of localized orbitals. The quantity $C$ is an empirical constant for a given LRC functional, which was determined in Ref. Modrzejewski:2013 for LRC-$\omega $PBE ($C = 0.90$) and LRC-$\omega $PBEh ($C = 0.75$), using the def2-TZVPP basis set. Since LRC-$\omega $PBE($\omega ^{}_{\rm GDD}$) provides a better description of polarizabilities in polyacetylene as compared to $\omega ^{}_{\rm IP}$[214], it is anticipated that using $\omega ^{}_{\rm GDD}$ in place of $\omega ^{}_{\rm IP}$ may afford more accurate molecular properties, especially in conjugated systems. The $\omega $ value based on the GDD tuning procedure for RSH functionals can be invoked by setting the $rem variable OMEGA_GDD = TRUE. Basic job-control variables and an example can be found below. Since $\omega ^{}_{\rm GDD}$ depends on the electron density, the converged density obtained by RSH functionals is used to compute $\omega ^{}_{\rm GDD}$. The range-separation parameter for RSH functionals is still necessary to set by the variable OMEGA, for example, OMEGA = 300 for LRC-$\omega $PBE and OMEGA = 200 for LRC-$\omega $PBEh although the final $\omega ^{}_{\rm GDD}$ value is not quite sensitive to the initial OMEGA input.

OMEGA_GDD

Controls the application of $\omega ^{}_{\rm GDD}$ tuning for long-range-corrected DFT


TYPE:

LOGICAL


DEFAULT:

FALSE


OPTIONS:

FALSE (or 0)

Do not apply $\omega ^{}_{\rm GDD}$ tuning.

TRUE (or 1)

Use $\omega ^{}_{\rm GDD}$ tuning.


RECOMMENDATION:

The $rem variable OMEGA must also be specified, in order to set the initial range-separation parameter.


OMEGA_GDD_SCALING

Sets the empirical constant $C$ in $\omega ^{}_{\rm GDD}$ tuning procedure.


TYPE:

INTEGER


DEFAULT:

885


OPTIONS:

$n$

Corresponding to $C = n/1000$.


RECOMMENDATION:

The quantity $n$ = 885 was determined by Lao and Herbert in Ref. Lao:2017 using LRC-$\omega $PBE and def2-TZVPP augmented with diffuse functions on non-hydrogen atoms that are taken from Dunning’s aug-cc-pVTZ basis set.


Example 4.27  Sample input illustrating a calculation to determine the $\omega $ value for LRC-$\omega $PBE based on the $\omega ^{}_{\rm GDD}$ tuning procedure.

$comment
    The initial omega value has to set.
$end

$rem
exchange          gen
basis             aug-cc-pvdz
lrc_dft           true
omega             300
omega_gdd         true
$end

$xc_functional
x   wPBE   1.0
c    PBE   1.0
$end

$molecule
0 1
O -0.042500 0.091700 0.110000
H  0.749000 0.556800 0.438700
H -0.825800 0.574700 0.432500
$end

4.4.7 DFT Methods for van der Waals Interactions

This section describes three different procedures for obtaining a better description of dispersion (i.e., van der Waals) interactions DFT calculations: non-local correlation functionals (Section 4.4.7.1), empirical atom–atom dispersion potentials (“DFT-D”, Section 4.4.7.2), and finally the Becke-Johnson exchange-dipole model (XDM, Section 4.4.7.3).

4.4.7.1 Non-Local Correlation (NLC) Functionals

Q-Chem includes four non-local correlation functionals that describe long-range dispersion:

All of these functionals are implemented self-consistently and analytic gradients with respect to nuclear displacements are available [216, 218, 35]. The non-local correlation is governed by the $rem variable NL_CORRELATION, which can be set to one of the four values: vdW-DF-04, vdW-DF-10, VV09, or VV10. Note that the vdW-DF-04, vdW-DF-10, and VV09 functionals are used in combination with LSDA correlation, which must be specified explicitly. For instance, vdW-DF-10 is invoked by the following keyword combination: 

CORRELATION        PW92
NL_CORRELATION     vdW-DF-10

VV10 is used in combination with PBE correlation, which must be added explicitly. In addition, the values of two parameters, $C$ and $b$ 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 $C = 0.0093$ and $b = 5.9$, 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, $C = 0.0089$ is used in VV09 (i.e. NL_VV_C is set to 89). However, in VV10 neither $C$ nor $b$ are assigned a default value and must always be provided in the input.

Unlike local (LSDA) and semi-local (GGA and meta-GGA) functionals, for non-local functionals evaluation of the correlation energy requires a double integral over the spatial variables [cf. Eq. ()]:

  \begin{equation} E_{\rm c}^{\rm nl} = \int f(\ensuremath{\mathbf{r}},\ensuremath{\mathbf{r'}}) \, d\ensuremath{\mathbf{r}} \; d\ensuremath{\mathbf{r'}}. \end{equation}   (4.56)

In practice, this double integration is performed numerically on a quadrature grid [216, 218, 35]. By default, the SG-1 quadrature (described in Section 4.4.5.2 below) is used to evaluate $E_{\rm c}^{\rm nl}$, but a different grid can be requested via the $rem variable NL_GRID. The non-local energy is rather insensitive to the fineness of the grid, so that SG-1 or even SG-0 grids can be used in most cases. However, a finer grid may be required to integrate other components of the functional, as controlled by the XC_GRID variable (see Section 4.4.5.2).

Several exchange-correlation functionals in Q-Chem are pre-defined to use the VV10 non-local correlation functional. The two functionals from the original paper[35] can be used by specifying METHOD = VV10 or METHOD LC-VV10. In addition, the combinatorially-optimized functionals of Mardirossian and Head-Gordon ($\omega $B97X-V, B97M-V, and $\omega $B97M-V) use non-local correlation and can be invoked by setting METHOD to wB97X-V, B97M-V, or wB97M-V.

Example 4.28  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              aug-cc-pVTZ
   EXCHANGE           rPW86
   CORRELATION        PBE
   XC_GRID            2
   NL_CORRELATION     VV10
   NL_GRID            1
   NL_VV_C            93
   NL_VV_B            590
$end

In the above example, the SG-2 grid is used to evaluate the rPW86 exchange and PBE correlation, but a coarser SG-1 grid is used for the non-local part of VV10. Furthermore, the above example is identical to specifying METHOD = VV10.

NL_CORRELATION

Specifies a non-local correlation functional that includes non-empirical dispersion.


TYPE:

STRING


DEFAULT:

None

No non-local correlation.


OPTIONS:

None

No non-local correlation

vdW-DF-04

the non-local part of vdW-DF-04

vdW-DF-10

the non-local part of vdW-DF-10 (also known as vdW-DF2)

VV09

the non-local part of VV09

VV10

the non-local part of VV10


RECOMMENDATION:

Do not forget to add the LSDA correlation (PW92 is recommended) when using vdW-DF-04, vdW-DF-10, 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 $C$ in VV09 and VV10. This parameter is fitted to asymptotic van der Waals $C_6$ coefficients.


TYPE:

INTEGER


DEFAULT:

89

for VV09

No default

for VV10


OPTIONS:

$n$

Corresponding to $C = n/10000$


RECOMMENDATION:

$C = 0.0093$ is recommended when a semi-local exchange functional is used. $C = 0.0089$ is recommended when a long-range corrected (LRC) hybrid functional is used. For further details see Ref. Vydrov:2010b.


NL_VV_B

Sets the parameter $b$ in VV10. This parameter controls the short range behavior of the non-local correlation energy.


TYPE:

INTEGER


DEFAULT:

No default


OPTIONS:

$n$

Corresponding to $b = n/100$


RECOMMENDATION:

The optimal value depends strongly on the exchange functional used. $b = 5.9$ is recommended for rPW86. For further details see Ref. Vydrov:2010b.


USE_RVV10

Used to turn on the rVV10 NLC functional


TYPE:

LOGICAL


DEFAULT:

FALSE


OPTIONS:

FALSE

VV10 NLC is utilized (the default for NL_CORRELATION)

TRUE

rVV10 NLC is utilized


RECOMMENDATION:

Set to TRUE if the rVV10 NLC is desired.


4.4.7.2 Empirical Dispersion Corrections: DFT-D

There are currently three unique DFT-D methods in Q-Chem, which are requested via the $rem variable DFT_D, and which are discussed below.

DFT_D

Controls the empirical dispersion correction to be added to a DFT calculation.


TYPE:

LOGICAL


DEFAULT:

None


OPTIONS:

FALSE

(or 0) Do not apply the DFT-D2, DFT-CHG, or DFT-D3 scheme

EMPIRICAL_GRIMME

DFT-D2 dispersion correction from Grimme

EMPIRICAL_CHG

DFT-CHG dispersion correction from Chai and Head-Gordon

EMPIRICAL_GRIMME3

DFT-D3(0) dispersion correction from Grimme (deprecated as of Q-Chem 5.0)

D3_ZERO

DFT-D3(0) dispersion correction from Grimme et al.

D3_BJ

DFT-D3(BJ) dispersion correction from Grimme et al.

D3_CSO

DFT-D3(CSO) dispersion correction from Schröder et al.

D3_ZEROM

DFT-D3M(0) dispersion correction from Smith et al.

D3_BJM

DFT-D3M(BJ) dispersion correction from Smith et al.

D3_OP

DFT-D3(op) dispersion correction from Witte et al.

D3

Recommended; chooses the "best" available D3 dispersion correction


RECOMMENDATION:

NONE


Thanks to the efforts of the Sherrill group, Grimm’s popular DFT-D2 empirical dispersion correction [28] is available in Q-Chem, along with analytic gradients and frequencies. This correction can be added to any density functional that is available in Q-Chem, although its use with the non-local correlation functionals discussed in Section 4.4.7.1 is obviously not recommended.

DFT-D2 adds an extra term,

  $\displaystyle \label{eqn:C6} E_{\rm disp} $ $\displaystyle = -s_6 \sum _ A \sum _{B<A} \frac{C_{6}^{AB}}{R_{AB}^{6}}f_\ensuremath{\mathrm{}}{damp}(R_{AB}) $   (4.57)
  $\displaystyle C_{6}^{AB} $ $\displaystyle = \bigl (C_{6}^{A} C_{6}^{B}\bigr )^{1/2} $   (4.58)
  $\displaystyle f_\ensuremath{\mathrm{}}{damp}(R_{AB}) $ $\displaystyle = \left[1+e^{-d(R_{AB}/R_{AB}^0 - 1)}\right]^{-1} $   (4.59)

where $s_6$ is a global scaling parameter (near unity) and $f_\ensuremath{\mathrm{}}{damp}$ is a damping function meant to help avoid double-counting correlation effects at short range (and to prevent the divergence of $R_{AB}^{-6}$), which depends on another parameter, $d$. The quantity $R_{AB}^0$ is the sum of the van der Waals radii of atoms A and B. Grimme has suggested using $s_6 = 0.75$ for PBE, $s_6 = 1.2$ for BLYP, $s_6 = 1.05$ for BP86, and $s_6 = 1.05$ for B3LYP, and these are the default values when those functionals are used. Otherwise, the default value is $s_6 = 1$.

DFT-D2 parameters, including the atomic van der Waals radii, can be altered using a $empirical_dispersion input section. 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 optimized by Grimme.

Note:  (i) DFT-D2 is only defined for elements up to Xe.
(ii) B97-D is an exchange-correlation functional that automatically employs the DFT-D2 dispersion correction when used via METHOD B97-D.

An alternative to Grimme’s DFT-D2 is the empirical dispersion correction of Chai and Head-Gordon [39], which employs a slightly different damping function

  $\displaystyle \label{eqn:CHG} f_{\rm damp}(R_{AB}) = \left[1+a(R_{AB}/R_{AB}^0)^{-12}\right]^{-1} $   (4.60)

The Chai–Head-Gordon dispersion correction is activated by setting DFT_D = EMPIRICAL_CHG, and the damping parameter $a$ is controlled by the keyword DFT_D_A.

DFT_D_A

Controls the strength of dispersion corrections in the Chai–Head-Gordon DFT-D scheme, Eq. eqn:CHG.


TYPE:

INTEGER


DEFAULT:

600


OPTIONS:

$n$

Corresponding to $a = n/100$.


RECOMMENDATION:

Use the default.


Note:  (i) DFT-CHG is only defined for elements up to Xe.
(ii) The $\omega $B97X-D and $\omega $M05-D functionals automatically employ the DFT-CHG dispersion correction when used via METHOD = wB97X-D or wM05-D.

Finally, Grimme’s more recent DFT-D3 method [40], and improvement on his DFT-D2 approach, [28] is available along with analytic first and second derivatives, and can be combined with any of the density functionals available in Q-Chem. The total DFT-D3 energy is given by $E_{\text {DFT-D3}} = E_{\text {KS-DFT}} + E_{\rm disp}$, where $E_{\text {KS-DFT}}$ is the total energy from KS-DFT and the dispersion correction $E_{\rm disp} = E^{(2)}+E^{(3)}$ is a sum of two- and three-body energies. As of Q-Chem 5.0, several versions of DFT-D3 are implemented.

All versions of DFT-D3 are based on the addition of an empirical dispersion term to the energy. This term is given by

  \begin{equation} E_{\text {D3}} = - \sum \limits _{i<j}\sum \limits _{n=6,8} s_ n \frac{C_{n,ij}}{r_{ij}^ n}f_{\text {damp},n}(r_{ij}), \label{D3eqn} \end{equation}   (4.61)

where the sum runs over all unique pairs of atoms $i$ and $j$; the $C_{n,ij}$ are $n$th-order isotropic dispersion coefficients for atom pair $ij$; $r_{ij}$ is the distance between atoms $i$ and $j$; $s_ n$ are global, functional-dependent scaling parameters; and $f_{\text {damp},n}$ are damping functions. The original zero-damping formulation, DFT-D3(0), of Grimme et al. [40] utilizes a damping function of the form

  \begin{equation} f_{\text {damp},n}^{\text {D3(0)}}(r_{ij}) = \left[1+6\left(\frac{r_{ij}}{s_{r,n}r_{0,ij}}\right)^{-\beta _ n}\right]^{-1}, \label{zerodamping} \end{equation}   (4.62)

where $r_{0,ij}$ are pre-tabulated van der Waals radii, $s_{r,6}$ is a functional-dependent parameter, $s_{r,8}=1$, $\beta _6=12$, and $\beta _8=14$. In this approach, $s_6$ is generally fixed to unity, and $s_8$ is optimized for each density functional.

The more recent BJ-damping model, [41] DFT-D3(BJ), adds in some finite amount of correlation even at short distances. Its damping function is given by

  \begin{equation} f_{\text {damp},n}^{\text {D3(BJ)}}(r_{ij}) = \frac{r_{ij}^ n}{r_{ij}^ n + \left(\alpha _1 r_{0,ij} + \alpha _2\right)^ n}, \label{bjdamping} \end{equation}   (4.63)

where $\alpha _1$ and $\alpha _2$ are adjustable nonlinear parameters fit to each density functional. As in DFT-D3(0), $s_6$ is generally fixed to unity and $s_8$ is optimized for each functional. The DFT-D3(BJ) approach generally outperforms the original DFT-D3(0) formulation. [41]

The C-six-only (CSO) approach of Schröder et al. [42] utilizes a damping function with one nonlinear parameter, $\alpha _1$:

  \begin{equation} f_{\text {damp},6}^{\text {D3(CSO)}}(r_{ij}) = \frac{r_{ij}^6}{r_{ij}^6+(2.5)^6}\left[1+\frac{\alpha _1}{s_6 \left[1 + \exp \left(r_{ij}-2.5 r_{0,ij}\right)\right]}\right]. \label{csodamping} \end{equation}   (4.64)

The DFT-D3(BJ) approach was re-parameterized by Smith et al. [43] to yield the DFT-D3M(BJ) approach. This approach was parameterized heavily on non-equilibrium systems. Those same authors altered and modified the original DFT-D3(0) to obtain DFT-D3M(0), which has one additional parameter on top of the original DFT-D3(0) set, $\alpha _1$:

  \begin{equation} f_{\text {damp},n}^{\text {D3(0)}}(r_{ij}) = \left[1+6\left(\frac{r_{ij}}{s_{r,n}r_{0,ij}}+\alpha _1 r_{0,ij}\right)^{-\beta _ n}\right]^{-1}. \label{zeromdamping} \end{equation}   (4.65)

The optimized power approach of Witte et al. [219] treats the exponent, $\beta _6$, as an optimizable parameter (Note: $\beta _8=\beta _6+2$), and is given by

  \begin{equation} f_{\text {damp},n}^{\text {D3(op)}}(r_{ij}) = \frac{r_{ij}^{\beta _ n}}{r_{ij}^{\beta _ n} + (\alpha _1 r_{0,ij} + \alpha _2)^{\beta _ n}}. \label{opdamping} \end{equation}   (4.66)

DFT-D3(0) is requested by setting DFT_D = D3_ZERO. The model depends on four scaling factors [$s_6$, $s_{r,6}$, $s_8$, and $s_{r,8}$, as defined in eq. (4.62), which can be adjusted using the $rem variables described below.

DFT-D3(BJ) is requested by setting DFT_D = D3_BJ. The model depends on four scaling factors [$s_6$, $s_8$, $\alpha _1$, and $\alpha _2$, as defined in eq. (4.63), which can be adjusted using the $rem variables described below.

DFT-D3(CSO) is requested by setting DFT_D = D3_CSO. The model depends on two scaling factors [$s_6$ and $\alpha _1$, as defined in eq. (4.64), which can be adjusted using the $rem variables described below.

DFT-D3M(0) is requested by setting DFT_D = D3_ZEROM. The model depends on five scaling factors [$s_6$, $s_8$, $s_{r,6}$, $s_{r,8}$, and $\alpha _1$, as defined in eq. (4.65), which can be adjusted using the $rem variables described below.

DFT-D3M(BJ) is requested by setting DFT_D = D3_BJM. The model depends on four scaling factors [$s_6$, $s_8$, $\alpha _1$, and $\alpha _2$, as defined in eq. (4.63), which can be adjusted using the $rem variables described below.

DFT-D3(op) is requested by setting DFT_D = D3_OP. The model depends on four scaling factors [$s_6$, $s_8$, $\alpha _1$, $\alpha _2$, and $\beta _6$ as defined in eq. (4.63), which can be adjusted using the $rem variables described below.

Or, one may simply set DFT_D = D3, and a D3 tail will be chosen automatically, if one is available for the selected functional.

Note:  (i) DFT-D3(0) is defined for elements up to Pu ($Z=94$).
(ii) The B97-D3(0), $\omega $B97X-D3, $\omega $M06-D3 functionals automatically employ the DFT-D3(0) dispersion correction when invoked by setting METHOD equal to B97-D3, wB97X-D3, or wM06-D3.
(iii) The old way of invoking DFT-D3, namely through the use of EMPIRICAL_GRIMME3, is still supported, though its use is discouraged, since D3_ZERO accomplishes the same thing, but with additional precision for the relevant parameters.
(iv) When DFT_D = D3, parameters may not be overwritten, with the exception of DFT_D3_3BODY; this is intended as a user-friendly option. This is also the case when EMPIRICAL_GRIMME3 is employed for a functional parameterized in Q-Chem. When any of D3_ZERO, D3_BJ, etc. are chosen, Q-Chem will automatically populate the parameters with values if available for the given functional, but will still allow these defaults to be overwritten.

DFT_D3_S6

The linear parameter $s_6$ in eq. (4.61). Utilized in all forms of DFT-D3.


TYPE:

INTEGER


DEFAULT:

100000


OPTIONS:

$n$

Corresponding to $s_6 = n/100000$.


RECOMMENDATION:

NONE


DFT_D3_RS6

The nonlinear parameter $s_{r,6}$ in eqs. (4.62) and (4.65). Utilized in DFT-D3(0) and DFT-D3M(0).


TYPE:

INTEGER


DEFAULT:

100000


OPTIONS:

$n$

Corresponding to $s_{r,6} = n/100000$.


RECOMMENDATION:

NONE


DFT_D3_S8

The linear parameter $s_8$ in eq. (4.61). Utilized in DFT-D3(0), DFT-D3(BJ), DFT-D3M(0), DFT-D3M(BJ), and DFT-D3(op).


TYPE:

INTEGER


DEFAULT:

100000


OPTIONS:

$n$

Corresponding to $s_8 = n/100000$.


RECOMMENDATION:

NONE


DFT_D3_RS8

The nonlinear parameter $s_{r,8}$ in eqs. (4.62) and (4.65). Utilized in DFT-D3(0) and DFT-D3M(0).


TYPE:

INTEGER


DEFAULT:

100000


OPTIONS:

$n$

Corresponding to $s_{r,8} = n/100000$.


RECOMMENDATION:

NONE


DFT_D3_A1

The nonlinear parameter $\alpha _1$ in eqs. (4.63), (4.64), (4.65), and (4.66). Utilized in DFT-D3(BJ), DFT-D3(CSO), DFT-D3M(0), DFT-D3M(BJ), and DFT-D3(op).


TYPE:

INTEGER


DEFAULT:

100000


OPTIONS:

$n$

Corresponding to $\alpha _1 = n/100000$.


RECOMMENDATION:

NONE


DFT_D3_A2

The nonlinear parameter $\alpha _2$ in eqs. (4.63) and (4.66). Utilized in DFT-D3(BJ), DFT-D3M(BJ), and DFT-D3(op).


TYPE:

INTEGER


DEFAULT:

100000


OPTIONS:

$n$

Corresponding to $\alpha _2 = n/100000$.


RECOMMENDATION:

NONE


DFT_D3_POWER

The nonlinear parameter $\beta _6$ in eq. (4.66). Utilized in DFT-D3(op). Must be greater than or equal to 6 to avoid divergence.


TYPE:

INTEGER


DEFAULT:

600000


OPTIONS:

$n$

Corresponding to $\beta _6 = n/100000$.


RECOMMENDATION:

NONE


The three-body interaction term [40], $E^{(3)}$, must be explicitly turned on, if desired.

DFT_D3_3BODY

Controls whether the three-body interaction in Grimme’s DFT-D3 method should be applied (see Eq. (14) in Ref. Grimme:2010).


TYPE:

LOGICAL


DEFAULT:

FALSE


OPTIONS:

FALSE (or 0)

Do not apply the three-body interaction term

TRUE

Apply the three-body interaction term


RECOMMENDATION:

NONE


Example 4.29  Applications of B3LYP-D3(0) with custom parameters to a methane dimer.

$comment
   Geometry optimization, followed by single-point 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           6-31G*
   DFT_D           D3_ZERO
   DFT_D3_S6       100000
   DFT_D3_RS6      126100
   DFT_D3_S8       170300
   DFT_D3_3BODY    FALSE
$end

@@@

$molecule
   read
$end

$rem
   JOBTYPE         sp
   EXCHANGE        B3LYP
   BASIS           6-311++G**
   DFT_D           D3_ZERO
   DFT_D3_S6       100000
   DFT_D3_RS6      126100
   DFT_D3_S8       170300
   DFT_D3_3BODY    FALSE
$end

4.4.7.3 Exchange-Dipole Model (XDM)

Becke and Johnson have proposed a conceptually simple yet accurate dispersion model called the exchange-dipole model (XDM) [168, 220]. In this model the dispersion attraction emerges from the interaction between the instantaneous 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 Q-Chem 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:

  \begin{equation} d_{\sigma }(\mathbf{r})=-\int h_{\sigma }(\mathbf{r},\mathbf{r}’) \; \mathbf{r}’ \; d\mathbf{r}’ - \mathbf{r} \end{equation}   (4.67)

where $\sigma $ labels the spin and $h_{\sigma }(\mathbf{r},\mathbf{r}’)$ is the exchange-hole function. The XDM version that is implemented in Q-Chem employs the Becke-Roussel (BR) model exchange-hole function. In previous work on the BR model, $h_\sigma $ was not available in analytical form and one had to determine its value numerically at each grid point. Q-Chem has developed for the first time an analytical expression for this function based on non-linear interpolation and spline techniques, which greatly improves efficiency as well as the numerical stability [221].

There are two different damping functions used in the XDM model of Becke and Johnson. One of them uses only the intermolecular $C_6$ dispersion coefficient, and its implementation in Q-Chem is denoted as "XDM6". In this version the dispersion energy is computed as

  \begin{equation} E_{\rm vdW}=\sum _{i>j} E_{{\rm vdW},ij} =-\sum _{i>j}\frac{C_{6,ij}}{R_{ij}^{6} +kC_{6,ij} (E_{i}^{C}+E_{j}^{C})^{-1}} \end{equation}   (4.68)

where $k$ is a universal parameter, $R_{ij}^{{}}$ is the distance between atoms $i$ and $j$, and $E_{ij}^{C}$ is the sum of the absolute values of the correlation energies of the free atoms $i$ and $j$. The dispersion coefficients $C_{6,ij}$ is computed as

  \begin{equation} C_{6,ij}=\frac{\langle d_{\rm X}^{2}\rangle ^{}_ i \langle d_{\rm X}^{2}\rangle ^{}_{j} \; \alpha ^{}_{i} \, \alpha ^{}_{j}}{\langle d_{\rm X}^{2}\rangle ^{}_{i}\alpha _{j}+\langle d_{\rm X}^{2}\rangle ^{}_{j} \alpha ^{}_{i}} \end{equation}   (4.69)

where $\langle d_{\rm X}^{2}\rangle ^{}_{i}$ is the exchange-hole dipole moment of the atom, and $\alpha ^{}_{i}$ is the effective polarizability of the atom $i$ in the molecule.

The XDM6 scheme can be further generalized to include higher-order dispersion coefficients, which leads to the “XDM10” model in Q-Chem, whose dispersion energy is

  \begin{equation} E_{\rm vdW}=-\sum _{i>j}\left(\frac{C_{6,ij}}{R_{{\rm vdW},ij}^{6}+R_{ij}^{6}} +\frac{C_{8,ij}}{R_{{\rm vdW},ij}^{8}+R_{ij}^{8}} +\frac{C_{10,ij}}{R_{{\rm vdW},ij}^{10}+R_{ij}^{10}}\right) \end{equation}   (4.70)

where $C_{6,ij}$, $C_{8,ij}$ and $C_{10,ij}$ are dispersion coefficients computed using higher-order multipole moments of the exchange hole, including dipole, quadrupole and octupole) [222]. The quantity $R_{{\rm vdW},ij}^{{}}$ is the sum of the effective vdW radii of atoms $i$ and $j$, which is a linear function

  \begin{equation} R_{{\rm vdW},ij}^{{}}=a_{1}R_{{\rm C},ij}^{{}}+a_{2} \end{equation}   (4.71)

of the so called critical distance $R_{{\rm C},ij}^{{}}$ between atoms $i$ and $j$, computed according to

  \begin{equation} R_{{\rm C},ij}^{{}} = \frac{1}{3}\left[ \left(\frac{C_{8,ij}}{C_{6,ij}}\right)^{1/2} +\left(\frac{C_{10,ij}}{C_{6,ij}}\right)^{1/4} +\left(\frac{C_{10,ij}}{C_{8,ij}}\right)^{1/2} \right] \; . \end{equation}   (4.72)

In the XDM10 scheme there are two universal parameters, $a_{1}$ and $a_{2}$. Their default values of 0.83 and 1.35, respectively, are due to Johnson and Becke [220], determined by least-squares fitting to the binding energies of a set of intermolecular complexes. These values are not the only possible optimal set to use with XDM. Becke’s group later suggested several other XC functional combinations with XDM that employ different values of $a_{1}$ and $a_{2}$. The user is advised to consult the recent literature for details [223, 224].

The computed vdW energy $E_{\rm vdW}$ is added as a post-SCF correction. Analytic gradients and Hessians are available for both XDM6 and XDM10. Additional job control and customization options are listed below.

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 XDM6).


RECOMMENDATION:

None


DFTVDW_METHOD

Choose the damping function used in XDM


TYPE:

INTEGER


DEFAULT:

1


OPTIONS:

1

Use Becke’s damping function including $C_6$ term only.

2

Use Becke’s damping function with higher-order ($C_8$ and $C_{10}$) terms.


RECOMMENDATION:

None


DFTVDW_MOL1NATOMS

The number of atoms in the first monomer in dimer calculation


TYPE:

INTEGER


DEFAULT:

0


OPTIONS:

0–$N_{\rm atoms}$


RECOMMENDATION:

None


DFTVDW_KAI

Damping factor $k$ for $C_6$-only damping function


TYPE:

INTEGER


DEFAULT:

800


OPTIONS:

10–1000


RECOMMENDATION:

None


DFTVDW_ALPHA1

Parameter in XDM calculation with higher-order terms


TYPE:

INTEGER


DEFAULT:

83


OPTIONS:

10-1000


RECOMMENDATION:

None


DFTVDW_ALPHA2

Parameter in XDM calculation with higher-order terms.


TYPE:

INTEGER


DEFAULT:

155


OPTIONS:

10-1000


RECOMMENDATION:

None


DFTVDW_USE_ELE_DRV

Specify whether to add the gradient correction to the XDM energy. only valid with Becke’s $C_6$ damping function using the interpolated BR89 model.


TYPE:

LOGICAL


DEFAULT:

1


OPTIONS:

1

Use density correction when applicable.

0

Do not use this correction (for debugging purposes).


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.30  Sample input illustrating a frequency calculation of a vdW complex consisted of He atom and N$_2$ molecule.

$molecule
   0 1
   He 0.000000  0.00000   3.800000
   N  0.000000  0.000000  0.546986
   N  0.000000  0.000000 -0.546986
$end

$rem
   JOBTYPE             FREQ
   IDERIV              2
   EXCHANGE            B3LYP
   INCDFT              0
   SCF_CONVERGENCE     8
   BASIS               6-31G*
   !vdw parameters settings
   DFTVDW_JOBNUMBER    1
   DFTVDW_METHOD       1
   DFTVDW_PRINT        0
   DFTVDW_KAI          800
   DFTVDW_USE_ELE_DRV  0
$end

XDM can be used in conjunction with different GGA, meta-GGA pure or hybrid functionals, even though the original implementation of Becke and Johnson was in combination with HF exchange, or with a specific meta-GGA exchange and correlation (the BR89 exchange and BR94 correlation functionals). For example, encouraging results have been obtained using XDM with B3LYP [221]. Becke has found more recently that this model can be efficiently combined with the P86 exchange functional, with the hyper-GGA functional B05. Using XDM together with PBE exchange plus LYP correlation, or PBE exchange plus BR94 correlation, has been also found fruitful.

4.4.8 Double-Hybrid Density Functional Theory

The recent advance in double-hybrid density functional theory (DH-DFT) [158, 225, 226, 227, 156], has demonstrated its great potential for approaching the chemical accuracy with a computational cost comparable to the second-order Møller-Plesset perturbation theory (MP2). In a DH-DFT, a Kohn-Sham (KS) DFT calculation is performed first, followed by a treatment of non-local orbital correlation energy at the level of second-order Møller-Plesset perturbation theory (MP2) [228]. This MP2 correlation correction includes a a same-spin (ss) component, $E_{c}^\ensuremath{\mathrm{}}{ss}$, as well as an opposite-spin (os) component, $E_{c}^\ensuremath{\mathrm{}}{os}$, which are added to the total energy obtained from the KS-DFT calculation. Two scaling parameters, $c_\ensuremath{\mathrm{}}{ss}$ and $c_\ensuremath{\mathrm{}}{os}$, are introduced in order to avoid double-counting correlation:

  \begin{equation} \label{eq:dft2a} E_{\text {DH-DFT}} = E_{\text {KS-DFT}} + c_\ensuremath{\mathrm{}}{ss} E_{c}^\ensuremath{\mathrm{}}{ss} + c_\ensuremath{\mathrm{}}{os} E_{c}^\ensuremath{\mathrm{}}{os} \end{equation}   (4.73)

A starting point for understanding where a functional form like Eq. eq:dft2a might come from is the adiabatic connection formula that provides a rigorous way to define double-hybrid functionals. One considers an adiabatic path between the fictitious non-interacting Kohn-Sham reference system ($\lambda = 0$) and the real physical system ($\lambda = 1$), while holding the electron density fixed at its value for the real system, for all $\lambda $. Then

  \begin{equation} \label{sung1} E_{\rm XC} \left[\rho \right]=\int _{0}^{1}U_{{\rm XC},\lambda }[\rho ] \; d\lambda \, , \end{equation}   (4.74)

where $U_{{\rm XC}, \lambda }$ is the exchange-correlation energy at a coupling strength $\lambda $, meaning that the exchange-correlation energy if the electron–electron terms in the Hamiltonian had the form $\lambda /r^{}_{ij}$ rather than $1/r^{}_{ij}$. Using a linear model of this quantity,

  \begin{equation} \label{sung2} U_{{\rm XC},\lambda } =a+b\lambda \; , \end{equation}   (4.75)

one obtains the popular hybrid functional that includes the HF exchange (or occupied orbitals) such as B3LYP [31, 32]. If one further uses the Gorling-Levy’s perturbation theory (GL2) to define the initial slope at $\lambda $= 0, one obtains the doubly hybrid functional form in Eq. eq:dft2a that includes MP2 type perturbative terms (PT2) involving virtual Kohn-Sham orbitals [229]:

  \begin{equation} \label{sung3} U_{{\rm XC},\lambda } =\left. \frac{\partial U_{{\rm XC},\lambda } }{\lambda } \right|_{\lambda =0} =2E_{\rm C}^{\rm GL2} \, . \end{equation}   (4.76)

The adiabatic connection formula has been used to develop double hybrid functionals such as XYG3 [156]. Note that XYG3, as implemented in Q-Chem, uses B3LYP orbitals to generate the density and evaluate the PT2 terms. This is different from B2PLYP, an earlier doubly hybrid functional from Grimme [158]. The latter uses truncated Kohn-Sham orbitals while XYG3 uses converged KS orbitals to evaluate the PT2 terms. The performance of XYG3 is not only comparable to that of the G2 or G3 theory for thermochemistry, but barrier heights and non-covalent interactions, including stacking interactions, are also very well described by XYG3 [156]. The recommended basis set for XYG3 is 6-311+G(3df,2p).

Due to the inclusion of PT2 terms, the cost of double-hybrid calculations is formally ${\cal {O}}({N^5})$, as in conventional MP2, thereby not applicable to large systems and partly losing DFT’s cost advantages. However, the highly successful SOS-MP2 and local SOS-MP2 algorithms in Q-Chem can be leveraged to develop double-hybrid functionals based on these ${\cal {O}}({N^4})$ methods. A version of XYG3 that uses SOS-MP2 is the XYGJ-OS functional [157]. This functional has 4 parameters that are optimized using thermochemical data. It is not only faster than XYG3, but comparable to XYG3 (or perhaps slightly better) in accuracy. If the local SOS-MP2 algorithm is applied, the scaling of XYGJ-OS is further reduced to ${\cal {O}}({N^3})$. Recently, XYGJ-OS became the first double-hybrid functional whose analytic energy gradient has been derived and implemented [230].

Other more empirical double-hybrid functionals have been implemented in Q-Chem. Among the $\omega $B97 series of functionals, $\omega $B97X-2 [155] is a long-range corrected double hybrid that can greatly reduce the self-interaction errors (due to its high fraction of HF exchange), and has been shown significantly superior to other functionals for systems with both bonded and non-bonded interactions. Due to the sensitivity of PT2 correlation energy with respect to the choices of basis sets, $\omega $B97X-2 was parameterized with two different basis sets: the $\omega $B97X-2(LP) was parameterized for use with 6-311++G(3df,3pd), while $\omega $B97X-2(TQZ) was parameterized with the T/Q (triple-$\zeta $/quadruple-$\zeta $) extrapolation to the basis set limit. A careful reading of Ref. Chai:2009 is highly advised before using either of these two functionals.

Job control variables for double-hybrid DFT are described below. Note that the PT2 correlation energy can also be computed with the efficient resolution-of-identity (RI) methods; see Section 5.6.

DH

Controls the application of DH-DFT scheme.


TYPE:

LOGICAL


DEFAULT:

FALSE


OPTIONS:

FALSE (or 0)

Do not apply the DH-DFT scheme

TRUE (or 1)

Apply DH-DFT scheme


RECOMMENDATION:

NONE


The following to $rem variables pertain to the $\omega $B97X-2(LP) and $\omega $B97X-2(TQZ) functionals.

SSS_FACTOR

Controls the strength of the same-spin component of PT2 correlation energy.


TYPE:

INTEGER


DEFAULT:

0


OPTIONS:

$n$

Corresponding to $c_\ensuremath{\mathrm{}}{ss} = n/10^6$ in Eq. (4.73).


RECOMMENDATION:

NONE


SOS_FACTOR

Controls the strength of the opposite-spin component of PT2 correlation energy.


TYPE:

INTEGER


DEFAULT:

0


OPTIONS:

$n$

Corresponding to $c_\ensuremath{\mathrm{}}{os} = n/10^6$ in Eq. (4.73).


RECOMMENDATION:

NONE


Example 4.31  Applications of B2PLYP functional to LiH.

$comment
   Single-point calculation on LiH without RI, followed by a single-point
   calculation with RI. 
$end

$molecule 
   0 1 
   H  
   Li H 1.6 
$end 

$rem
   JOBTYPE     sp
   EXCHANGE    B2PLYP
   BASIS       cc-pvtz
$end

@@@

$molecule
   read
$end

$rem
   JOBTYPE     sp
   EXCHANGE    B2PLYP
   BASIS       cc-pvtz
   AUX_BASIS   rimp2-cc-pvtz
$end

Example 4.32  Application of $\omega $B97X-2(TQZ) functional (with and without RI) to LiH molecules.

$comment
   Single-point calculation on LiH (with RI).
$end

$molecule
   0 1
   H
   Li H 1.6
$end

$rem
   JOBTYPE     sp
   EXCHANGE    omegaB97X-2(TQZ)
   BASIS       cc-pvqz
   AUX_BASIS   rimp2-cc-pvqz
$end

@@@

$comment
   Single-point calculation on LiH (without RI).
$end

$molecule
   read
$end

$rem
   JOBTYPE      sp
   EXCHANGE      omegaB97X-2(TQZ)
   BASIS      cc-pvqz
$end

In the following example of XYG3, Q-Chem automatically performs a B3LYP calculation first, then uses the resulting orbitals to evaluate the PT2 correlation terms. Once can also use XYG3 combined with the RI approximation; use EXCHANGE = XYG3RI to do so, along with an appropriate choice of auxiliary basis set.

Example 4.33  XYG3 calculation of N$_2$

$molecule
   0 1
   N      0.000000   0.000000   0.547775
   N      0.000000   0.000000  -0.547775
$end

$rem
   EXCHANGE  xyg3
   BASIS     6-311+G(3df,2p)
$end

The next example illustrates XYGJ-OS. This functional uses the RI approximation by default, so it is necessary to specify an auxiliary basis set.

Example 4.34  XYGJ-OS calculation of N$_2$

$molecule
  0 1
  N      0.000000   0.000000   0.547775
  N      0.000000   0.000000  -0.547775
$end

$rem
   EXCHANGE    xygjos
   BASIS       6-311+G(3df,2p)
   AUX_BASIS   rimp2-cc-pVtZ
   PURECART    1111
   TIME_MP2    true
$end

The final example uses the local version of XYGJ-OS, which is the same as the original XYGJ-OS but with the use of the attenuated Coulomb metric to solve the RI coefficients. Here, the keyword omega determines the locality of the metric.

Example 4.35  Local XYGJ-OS calculation of N$_2$

$molecule
   0 1
   N    0.000    0.000   0.54777500
   N    0.000    0.000  -0.54777500
$end

$rem
   EXCHANGE   lxygjos
   OMEGA      200
   BASIS      6-311+G(3df,2p)
   AUX_BASIS  rimp2-cc-pVtZ
   PURECART   1111
$end

4.4.9 Asymptotically Corrected Exchange-Correlation Potentials

No GGA exchange functional can simultaneously produce the correct contribution to the exchange energy density and exchange potential in the asymptotic region of molecular systems [231], and existing GGA exchange-correlation (XC) potentials decay much faster than the correct $-1/r$ XC potential in the asymptotic region [232]. High-lying occupied orbitals and low-lying virtual orbitals are therefore too loosely bound, and $-\varepsilon ^{}_{\rm HOMO}$ becomes far smaller than the ionization energy, despite the exact condition that these should be the same for the exact functional [205, 206]. Moreover, response properties may be poorly predicted from TDDFT calculations with GGA functionals [205]. Long-range corrected hybrid DFT (LRC-DFT), described in Section 4.4.6, has greatly remedied this situation, but is more expensive that KS-DFT with GGA functionals due to the use of Hatree-Fock exchange. The asymptotic corrections described in this section are designed to remedy the same problems but within the GGA framework.

4.4.9.1 LB94 Scheme

An asymptotically corrected (AC) exchange potential proposed by van Leeuwen and Baerends is [231]

  \begin{equation} \label{eq:vx_ LB} v_{x}^{\text {LB}} = -\beta \left( \frac{x^2}{1+3 \beta \mbox{sinh}^{-1}(x)}\right) \end{equation}   (4.77)

where $x = \frac{ \vert \nabla \rho \vert }{\rho ^{4/3}}$ is the reduced density gradient. For an exponentially-decaying density, this potential reduces to $-1/r$ in the asymptotic region of molecular systems. The LB94 xc potential is formed by a linear combination of LDA XC potential and the LB exchange potential [231]:

  \begin{equation} v_{xc}^{\text {LB94}} = v_{xc}^{\text {LDA}} + v_{x}^{\text {LB}} \; . \end{equation}   (4.78)

The parameter $\beta $ in Eq. eq:vx_LB was determined by fitting to the exact XC potential for Be atom. As mentioned in Refs. Casida:1998 and Hirata:1999b, for TDDFT calculations, it is sufficient to include the AC XC potential for ground-state calculations followed by TDDFT calculations with an adiabatic LDA XC kernel. The implementation of the LB94 XC potential in Q-Chem takes this approach, using the LB94 XC potential for the ground state calculations, followed by a TDDFT calculation with an adiabatic LDA XC kernel. This TDLDA/LB94 approach has been widely applied to study excited-state properties of large molecules.

Since the LB exchange potential in Eq. eq:vx_LB does not come from the functional derivative of an exchange energy functional, the Levy-Perdew virial relation [235] is used instead to obtain the exchange energy:

  \begin{equation} E_ x^{\text {LB}} = -\int v_ x^{\text {LB}}([\rho ],\textbf{r})[3\rho (\textbf{r})+\textbf{r}\nabla \rho (\textbf{r})]d\textbf{r} \end{equation}   (4.79)

An LB94 calculation is requested by setting EXCHANGE = LB94 in the $rem section. Additional job control and examples appear below.

LB94_BETA

Sets the $\beta $ parameter for the LB94 XC potential


TYPE:

INTEGER


DEFAULT:

500


OPTIONS:

$n$

Corresponding to $\beta = n/10000$.


RECOMMENDATION:

Use the default.


Example 4.36  Applications of LB94 XC potential to N$_2$ 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       = 6-311(2+,2+)G**
   CIS_N_ROOTS = 30
   RPA         = true
$end

4.4.9.2 Localized Fermi-Amaldi (LFA) Schemes

Another alternative, proposed by Pan, Fang and Chai [236], is to use a localized version of Fermi-Amaldi exchange-correlation functional. The resulting exchange density functional, whose functional derivative has the correct $-1/r$ asymptotic behavior, can be directly added to any semi-local density functional. Three variants of this method were proposed in Ref. Chai:2013b. The simplest of these, the strictly-localized Fermi-Amaldi (LFAs) scheme, is implemented in Q-Chem, for molecules consisting of atoms with $Z\le 55$.

Example 4.37  LFAs-PBE single-point TD-DFT calculation with water molecule

$comment
   Use LFAs-PBE potential for ground-state 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         6-311(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

4.4.10 Derivative Discontinuity Restoration

From the perspective of perturbation theory, Chai and Chen [237] proposed a systematic procedure for the evaluation of the derivative discontinuity of the exchange-correlation energy functional in KS-DFT, wherein the exact derivative discontinuity 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 derivative discontinuity. In particular, the first-order correction term is equivalent to the frozen-orbital approximation method. Its implementation in Q-Chem supports only local and GGA functionals at present, not meta-GGA, hybrid, or non-local functionals. Job control variables and examples appear below.

FOA_FUNDGAP

Compute the frozen-orbital approximation of the fundamental gap.


TYPE:

Boolean


DEFAULT:

FALSE


OPTIONS:

FALSE

Do not compute FOA derivative discontinuity and fundamental gap.

TRUE

Compute and print FOA fundamental gap information. Implies KS_GAP_PRINT.


RECOMMENDATION:

Use in conjunction with KS_GAP_UNIT if true.


Example 4.38  frozen-orbital approximation of derivative discontinuity with PBE and LFAs-PBE functionals on carbon atom

$comment
    Frozen-orbital derivative discontinuity, C atom, PBE
$end

$molecule
   0 3
   C
$end

$rem
   JOBTYPE       sp
   BASIS         6-31G*
   METHOD        PBE
   FOA_FUNDGAP   true
   KS_GAP_UNIT   1    ! print gap info in eV
$end

@@@

$comment
   with LFAs-PBE functional instead
$end

$molecule
   READ
$end

$rem
   JOBTYPE       sp
   BASIS         6-31G*
   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

4.4.11 Thermally-Assisted-Occupation Density Functional Theory (TAO-DFT)

Aiming to study the ground-state properties of large, strongly correlated systems with minimum computational complexity, Prof. Jeng-Da Chai recently developed thermally-assisted-occupation density functional theory (TAO-DFT) [238]. Unlike conventional multi-reference methods, the computational complexity of TAO-DFT increases very insignificantly with the size of the active space (i.e., an active space restriction is not needed for TAO-DFT calculations), and TAO-DFT appears to be very promising for the study of large poly-radical systems. TAO-DFT is a DFT scheme with fractional orbital occupations produced by the Fermi-Dirac distribution, controlled by a fictitious temperature $\theta $, and existing XC functionals (e.g., LDA or GGAs) can be used in TAO-DFT [239]. The computational cost of the method is similar to that of KS-DFT for single-point energy calculations and analytical nuclear gradients, and reduces to the cost of KS-DFT in the absence of strong static correlation effects.

There are several $rem variables that are used for TAO-DFT.

TAO_DFT

Controls whether to use TAO-DFT.


TYPE:

Boolean


DEFAULT:

false


OPTIONS:

false

Do not use TAO-DFT

true

Use TAO-DFT


RECOMMENDATION:

NONE


TAO_DFT_THETA

Controls the value of the fictitious temperature $\theta $ in TAO-DFT.


TYPE:

INTEGER


DEFAULT:

7


OPTIONS:

$m$

$\theta =m\times 10^{-n}$ (hartrees), where $n$ is the value of TAO_DFT_THETA_NDP


RECOMMENDATION:

NONE


TAO_DFT_THETA_NDP

Controls the value of the fictitious temperature $\theta $ in TAO-DFT.


TYPE:

INTEGER


DEFAULT:

3


OPTIONS:

$n$

$\theta =m\times 10^{-n}$ (hartrees), where $m$ is the value of TAO_DFT_THETA


RECOMMENDATION:

NONE


Note that setting TAO_DFT_THETA = 0 recovers ordinary KS-DFT [238]. In addition to the XC functional, a functional $E_{\theta }[\rho ]$ is needed in TAO-DFT. Currently available in Q-Chem are an LDA version [238] (the ETheta_LDA functional) as well as a version based on the gradient expansion approximation (GEA) [239] (ETheta_GEA functional), and the latter may be substituted for the former in the sample jobs below.

Example 4.39  TAO-LDA calculation on Be atom

$molecule
   0 1
   Be
$end

$rem
   JOBTYPE             sp
   BASIS               6-31G*
   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.40  TAO-PBE, spin-restricted calculation on stretched N$_2$

$molecule
   0 1
   N1
   N2 N1 5
$end

$rem
   JOBTYPE             sp
   BASIS               6-31G*
   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.41  TAO-PBE, spin-unrestricted calculation on stretched N$_2$

$molecule
   0 1
   N1
   N2 N1 5
$end

$rem
   JOBTYPE             opt
   UNRESTRICTED        true
   BASIS               6-31G*
   EXCHANGE            gen
   TAO_DFT             true
   TAO_DFT_THETA       40  ! theta = 40 mhartrees
   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