5.6 Range-Separated Hybrid Density Functionals

5.6.3 Tuned RSH Functionals

(June 30, 2021)

Whereas the range-separation parameters for the functionals described in Section 5.6.1 are wholly empirical in nature and should not be adjusted, for the functionals described in Section 5.6.2 some adjustment was suggested, especially for TDDFT calculations and for any properties that require interpretation of orbital energies, such as HOMO/LUMO gaps. This adjustment can be performed in a non-empirical (albeit system-specific) way, 47 Baer R., Livshits E., Salzner U.
Annu. Rev. Phys. Chem.
(2010), 61, pp. 85.
Link
by “tuning” the value of ω in order to satisfy the Koopmans-like ionization energy criterion

-ϵHOMO(ω)=IE(ω) (5.18)

where IE=E(N)-E(N-1) is the ΔSCF value of the ionization energy for the N-electron system of interest. The condition ϵHOMO=-IE is a theorem in exact DFT, 842 Perdew J. P., Levy M.
Phys. Rev. B
(1997), 56, pp. 16021.
Link
but this condition is often badly violated by approximate functionals. When an RSH functional is used, both sides of Eq. (5.18) are ω-dependent and this parameter is adjusted until the condition in Eq. (5.18) is met, which requires a sequence of SCF calculations on both the neutral and ionized species, using different values of ω. The value that is obtained has come to be called the “optimally tuned” value of ω. Formally speaking, there is no guarantee that an approximate density functional can be made to satisfy Eq. (5.18) for any given molecule, thus the optimally-tuned value need not exist. In practice it is usually possible to find such a value, although it should be noted that the optimally-tuned value of ω depends on system size, and as a result this tuning procedure formally violates size-consistency. 534 Karolewski A., Kronik L., Kümmel S.
J. Chem. Phys.
(2013), 138, pp. 204115.
Link

A few variations on the simple “IE tuning” criterion in Eq. (5.18) are possible. For proper description of charge-transfer states, Baer and co-workers 47 Baer R., Livshits E., Salzner U.
Annu. Rev. Phys. Chem.
(2010), 61, pp. 85.
Link
suggest finding the value of ω that (to the extent possible) satisfies Eq. (5.18) for both the neutral donor molecule and (separately) for the anion of the acceptor species. Along similar lines, in an effort to set both the HOMO and LUMO energy levels such that the fundamental gap (IE-EA) is equal to the HOMO/LUMO gap, Kronik et al. 577 Kronik L. et al.
J. Chem. Theory Comput.
(2012), 8, pp. 1515.
Link
suggest minimizing the function

J(ω)=[ϵHOMO(ω)+IE(ω)]2+[ϵLUMO(ω)+EA(ω)]2 (5.19)

with respect to ω, where EA=E(N+1)-E(N). Minimization of J(ω) represents an attempt to satisfy the IE theorem of Eq. (5.18) for both the N-electron molecule and its (N+1)-electron anion, assuming that the latter is bound. Published benchmarks suggest that these system-specific approaches afford the most accurate values of IEs and TDDFT excitation energies. 689 Livshits E., Baer R.
Phys. Chem. Chem. Phys.
(2007), 9, pp. 2932.
Link
, 957 Salzner U., Baer R.
J. Chem. Phys.
(2009), 131, pp. 231101.
Link
, 47 Baer R., Livshits E., Salzner U.
Annu. Rev. Phys. Chem.
(2010), 61, pp. 85.
Link
, 577 Kronik L. et al.
J. Chem. Theory Comput.
(2012), 8, pp. 1515.
Link

A script that optimizes ω, called OptOmegaIPEA.pl, is located in the $QC/bin/tools directory. The script scans ω 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, the user must create three inputs for an RSH single-point energy calculation, using the same geometry and basis set: one for a neutral molecule (N.in), one for its anion (M.in), and one for the molecule’s cation (P.in). The user should then run the command

OptOmegaIPEA.pl >& optomega

This command both generates the input files (N_*, P_*, M_*) and also runs Q-Chem on these input files, writing the optimization output into optomega. This script applies the IE condition to both the neutral molecule and its anion, minimizing J(ω) in Eq. (5.19). A similar script, OptOmegaIP.pl, uses Eq. (5.18) for the neutral molecule only.

Note:  1. 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 ω in such cases. 2. 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, 689 Livshits E., Baer R.
Phys. Chem. Chem. Phys.
(2007), 9, pp. 2932.
Link
, 957 Salzner U., Baer R.
J. Chem. Phys.
(2009), 131, pp. 231101.
Link
, 47 Baer R., Livshits E., Salzner U.
Annu. Rev. Phys. Chem.
(2010), 61, pp. 85.
Link
it can equally well be applied to any RSH functional, as for example LRC-ωPBE (see, Ref. 1109). The aforementioned scripts will work with these other RSH functionals as well.

Unfortunately, optimally-tuned values of “ωIE”, obtained using the criterion in Eq. (5.18) exhibit a troubling dependence on system size, 570 Körzdörfer T. et al.
J. Chem. Phys.
(2011), 135, pp. 204107.
Link
, 1109 Uhlig F. et al.
J. Phys. Chem. A
(2014), 118, pp. 7507.
Link
, 340 Garrett K. et al.
J. Chem. Theory Comput.
(2014), 10, pp. 3821.
Link
, 1122 Vazquez X. A. S., Isborn C. M.
J. Chem. Phys.
(2015), 143, pp. 244105.
Link
, 218 Coons M. P., You Z.-Q., Herbert J. M.
J. Am. Chem. Soc.
(2016), 138, pp. 10879.
Link
, 826 Oviedo M. B., Ilawe N. V., Wong B. M.
J. Chem. Theory Comput.
(2016), 12, pp. 3593.
Link
leading to a loss of size-extensivity. 534 Karolewski A., Kronik L., Kümmel S.
J. Chem. Phys.
(2013), 138, pp. 204115.
Link
For example, the optimally-tuned value for the cluster anion (H2O)6- is very different than the one tuned for (H2O)70-, 1109 Uhlig F. et al.
J. Phys. Chem. A
(2014), 118, pp. 7507.
Link
and the optimally-tuned value also varies with conjugation length for π-conjugated systems. 570 Körzdörfer T. et al.
J. Chem. Phys.
(2011), 135, pp. 204107.
Link
An alternative to the IE-based criterion in Eq. (5.18) is the global density-dependent (GDD) tuning procedure, 767 Modrzejewski M. et al.
J. Phys. Chem. A
(2013), 117, pp. 11580.
Link
in which the optimal value

ωGDD=Cdx2-1/2 (5.20)

is related to the average of the distance dx between an electron in the outer regions of a molecule and the exchange hole in the region of localized valence orbitals. The quantity C is an empirical parameter for a given LRC functional, which was determined for LRC-ωPBE (C=0.90) and LRC-ωPBEh (C=0.75) using the def2-TZVPP basis set. 767 Modrzejewski M. et al.
J. Phys. Chem. A
(2013), 117, pp. 11580.
Link
(A slightly different value, C=0.885, was determined for Q-Chem’s implementation of LRC-ωPBE. 617 Lao K. U., Herbert J. M.
J. Chem. Theory Comput.
(2018), 14, pp. 2955.
Link
) Since LRC-ωPBE(ωGDD) provides a better description of polarizabilities in polyacetylene as compared to ωIE, 412 Hapka M. et al.
J. Chem. Phys.
(2014), 141, pp. 134120.
Link
it is anticipated that using ωGDD in place of ωIE may afford more accurate molecular properties, especially in conjugated systems. GDD tuning of an RSH functional is involving by setting the $rem variable OMEGA_GDD = TRUE. The electron density is obviously needed to compute ωGDD in Eq. (5.20) and this is accomplished using the converged SCF density computed using the RSH functional with the value of ω given by the $rem variable OMEGA. The value of ωGDD therefore depends, in principle, upon the value of OMEGA, although in practice it is not very sensitive to this value.

OMEGA_GDD
       Controls the application of ωGDD tuning for long-range-corrected DFT
TYPE:
       LOGICAL
DEFAULT:
       FALSE
OPTIONS:
       FALSE (or 0) Do not apply ωGDD tuning. TRUE (or 1) Use ω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 ω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. 617 using LRC-ωPBE and def2-TZVPP augmented with diffuse functions on non-hydrogen atoms that are taken from Dunning’s aug-cc-pVTZ basis set.

Example 5.8  Sample input illustrating a calculation to determine the ω value for LRC-ωPBE based on the ωGDD tuning procedure.

$comment
   The initial omega value has to set.
$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


$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

View output

However the tuning is accomplished, these tuned functionals are generally thought to work by reducing 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 213 Cohen A. J., Mori-Sánchez P., Yang W.
Science
(2008), 321, pp. 792.
Link
, 770 Mori-Sánchez P., Cohen A. J.
Phys. Chem. Chem. Phys.
(2014), 16, pp. 14378.
Link
) at integer values of N, but in practice these E(N) plots bow away from straight-line segments. 213 Cohen A. J., Mori-Sánchez P., Yang W.
Science
(2008), 321, pp. 792.
Link
Examination of such plots has been suggested as a means to adjust the fraction of short-range exchange in an RSH functional, 45 Autschbach J., Srebro M.
Acc. Chem. Res.
(2014), 47, pp. 2592.
Link
while the ω parameter is set by tuning.

FRACTIONAL_ELECTRON
       Add or subtract a fraction of an electron.
TYPE:
       INTEGER
DEFAULT:
       0
OPTIONS:
       0 Use an integer number of electrons. n Add n/1000 electrons to the system.
RECOMMENDATION:
       Use only if trying to generate E(N) plots. If n<0, a fraction of an electron is removed from the system.

Example 5.9  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 F2-.

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

$molecule
   -2 2
   F
$end

$rem
   EXCHANGE              b3lyp
   BASIS                 6-31+G*
   FRACTIONAL_ELECTRON  -500   ! divide by 1000 to get the fraction, -0.5 here.
   GEN_SCFMAN            FALSE ! not yet available in new scf code
$end

View output