Letter pubs.acs.org/JCTC
Cite This: J. Chem. Theory Comput. XXXX, XXX, XXX−XXX
Solvation-Mediated Tuning of the Range-Separated Hybrid Functional: Self-Sufficiency through Screened Exchange Bora Joo,† Herim Han,† and Eung-Gun Kim* Department of Polymer Science and Engineering, Dankook University, Yongin, Gyeonggi 16890, Korea S Supporting Information *
ABSTRACT: We propose a simple procedure that restores the ionization potential theorem as the sole tuning criterion for both the long- and short-range Fock exchange of the rangeseparated hybrid functional. The procedure works by screening out an opposing effect of the short-range Fock fraction at long range, through the 1/εr dielectric correction in combination with a popular continuum solvation model. Our method proves to be a consistent and accurate way of tuning for both the isolated and solvated molecules.
P
from the dielectric constant for the solvation model, was turned into a tunable parameter.9 As for range separation and screening, it is also important to mention screened hybrids, in which a constant fraction of Fock exchange is applied to the short (HSE10) or middle range (HISS11) only. HSE, in particular, has been applied with great success to metals and inorganic semiconductors. As the optimally tuned RSH finds widespread use as a new standard functional of choice, it has become evident that there are two outstanding issues to be addressed. First, unlike the range-separation parameter γ for the long range, the fractional Fock parameter α for the short range seems to require more than the IP theorem for its optimization. The use of piecewise linearity7,12 and tuning at multiple MO levels13,14 help find the optimal value, but these methods have not always worked.7,13 Data obtained at a high level of theory such as GW have been also used as a reference for tuning.15 Second, when combined with a continuum solvation model, enforcing the IP theorem on γ-tuning makes the value of γ too small. Too small a range-separation parameter means that the long-range Fock correction does not take effect until at a very large distance, leaving the base functional (GGA, for instance) intact for the most part. Although proper dielectric effects on the IP/EA (electron affinity) and exciton binding energy were produced for some representative organic π-conjugated molecules,16−18 it comes as no surprise that the solvation reintroduces the delocalization error and essentially restores the quasi-particle spectra to GGA quality.19 In this work, we address simultaneously the two issues involving the fractional Fock parameter α and solvation, simply by noticing that one holds a key to solving the other and vice versa. By effectively screening out the role of α at long range via solvation, we show that self-sufficiency is achieved with the IP
robably the most notable advances in density functional theory (DFT) in recent years have been brought by what is now collectively known as the optimally tuned rangeseparated hybrid (RSH) functional.1 The RSH functional, when its tunable parameter is determined in some non-empirical way, allows the prediction of the fundamental gaps2 and charge transfer states3 of finite systems at an unprecedented level of accuracy. The formalism combines two rather simple, yet powerful, ideas that require the following: (i) the long-range part of the exchange energy term follows the asymptotic 1/r behavior;4 and (ii) the HOMO energy equals the ionization potential (IP) in magnitude.5 The first requirement is fulfilled by introducing a switching function f(γr) that gradually converts the long-range tail of a standard density functional to the correct form (a scheme known as the range separation).4,6 In the most common implementation, the Coulomb operator is divided into the short-range (SR) and long-range (LR) terms as 1 − f (γr ) f (γr ) 1 = SR + LR = + r r r
(1)
The second requirement is used as the criterion for deciding where to start switching (the optimal tuning of the rangeseparation parameter γ).2,3 Further improvement has been made by adding the Fock exchange to the short range as well6 such that SR + LR =
1 − [α + βf (γr )] α + βf (γr ) + r r
(2) 7
which produced the quasi-particle spectra accurately. The short-range Fock exchange, introduced as a fixed fractional amount via another tunable parameter α, also allowed the application of the RSH to extended systems through an explicit dielectric constant ε so as to achieve the 1/εr decay.8 The screened RSH was later extended to a continuum solvation model, in which α was taken as a constant while ε, separate © XXXX American Chemical Society
Received: January 18, 2018 Published: May 15, 2018 A
DOI: 10.1021/acs.jctc.8b00049 J. Chem. Theory Comput. XXXX, XXX, XXX−XXX
Letter
Journal of Chemical Theory and Computation theorem as the sole criterion for optimizing the RSH. The RSH thus optimized reproduces extremely well experimental and GW data of PTCDA and NTCDA, some of the toughest molecules for DFT to treat.7,15 We used the functional form of eq 2 first proposed by Yanai et al.,6 except that an exponential function was used in place of the more common complementary error function such that f(γr) = 1 − exp(−γr).20 For a given pair of α and β, the rangeseparation parameter γ was tuned (denoting the optimal value as γ*) such that the IP theorem is achieved as closely as possible by minimizing3 J 2 (γ , α) = [HOMO(γ , α) + IP(γ , α)]2 + [LUMO(γ , α) + EA(γ , α)]2
(3)
where HOMO denotes the HOMO energy of the neutral molecule and LUMO the SOMO energy of the anion in neutral geometry. IP and EA are the vertical values of the ionization potential and electron affinity, respectively. To choose α−β pairs, we followed the prescription given by Kronik and coworkers:7
α+β=1
(4)
which ensures the correct asymptotic 1/r behavior at long range. The B88 exchange functional and the correlation functional of B3LYP (0.81 LYP + 0.19 VWN5) were used, as with the original CAM-B3LYP.6 The TZP (triple-ζ with one set of polarization functions) Slater-type orbital basis set was used. The conductor-like screening model (COSMO)21,22 was used as a continuum solvation model for dielectric screening, with a solvent radius of 3.28 Å. The geometries of the molecules were obtained using B3LYP (Grimme’s dispersion correction with the Becke−Johnson damping function23 was added for PTCDA), and used for RSH single-point calculations without further geometry optimization. All calculations were performed using the ADF code.24 We first discuss how the magnitude of J(γ*, α) is affected by the choice of α. We are reminded that γ* varies as a function of α (see Supporting Information for Figure S1). Figure 1a shows that J(γ*, α), found to be smallest at α = 0, increases monotonically with increasing α. The fractional short-range Fock exchange plays an important role in the ordering of deeplying MO levels for both PTCDA and NTCDA.7,15 Obviously, α = 0 does not produce the correct MO ordering as predicted by GW calculations7,15 (see Figure 1b for comparison). Assessing the deviation from piecewise linearity allows us to look into what α does at individual MO levels. This can be done by calculating the total energy difference ΔE between the actual and the intrinsic piecewise-linear curve:15
Figure 1. Impact of introducing the fraction α of short-range Fock exchange for PTCDA and NTCDA in the gas phase. (a) Extent of deviation from the IP theorem; (b) MO ordering of PTCDA; and (c) (de)localization error with fractional changes δ in the number of electrons of NTCDA. In panel b, the MO ordering at α = 0.28, with HOMO−2 and HOMO−3 nearly degenerate, is the same as that of GW calculations.7,15
ΔE(N + δ) = Ε(N + δ) − {[Ε(N + 1) − Ε(N )]δ + Ε(N )} (5)
instead of finding better ways to optimize α, by examining how the IP theorem responds when the molecule becomes solvated. Combining the RSH with a continuum solvation model requires special care. Standard continuum solvation models work by adding to the gas-phase Hamiltonian an external potential due to the surface charges induced by the charge distribution of the molecule.25 The resulting solvation, missing out on the correlation-induced screening, has little to no effect on the orbital energies, as noted in earlier attempts to optimize γ in combination with solvation models.19,26,27 The screening
where N is the number of electrons and δ the fraction of an electron (0 ≤ δ < 1). For NTCDA, the minimal (de)localization error ΔE is found, yet again, at α = 0 (Figure 1c). Furthermore, the errors in some ranges seem to turn around and start to subside with α approaching 0.5, in contrast to the ever increasing J(γ*, α). It is almost tempting to say, why does the IP theorem as a tuning criterion break down in the first place as soon as α is introduced? Below we seek to get at this outstanding question, B
DOI: 10.1021/acs.jctc.8b00049 J. Chem. Theory Comput. XXXX, XXX, XXX−XXX
Letter
Journal of Chemical Theory and Computation effect by the electron correlation can be added by modifying eq 4 such that α+β=
1 ε
(6)
which produces the screened asymptotic 1/εr behavior correctly. Screening the RSH this way was first proposed by Kronik and co-workers for use with periodic systems such as organic crystals,8 and here we extend their prescription to treating continuum dielectric media. In the original approach, no additional tuning of the parameters was necessary because the optimal values of both γ and α were taken from the isolated-molecule calculations and β was adjusted to accommodate the dielectric constant of the crystal according to eq 6. β may well be negative when the surroundings, unlike organic crystals, have huge dielectric constants, as recognized in a study of metal−organic interfaces.28 Our solvation calculations only require a re-evaluation of the J function on the same set of (γ*, α) pairs from the isolated molecule, this time by applying COSMO and α + β = 1/ε. Because the (γ*, α) pairs are not recalculated, the number of additional single-point calculations equals the number of α values. The entire procedure, including the γ-tuning step for the isolated molecule, is summarized in a flowchart in Scheme 1. Scheme 1. Solvation-mediated tuning procedure in a flowchart. α = 0.5 is an empirically set upper limit based on a set of molecules we have tested so far. Pushing the limit too high can also cause a compatibility issue with the semilocal correlation expression.29 In practice, setting up a dense twodimensional grid can be avoided entirely by limiting α to between 0.25 and 0.35 and using the golden proportion approach30 for γ-tuning
Figure 2. Impact of solvation on J(γ*, α) (a) and the extent of (de)localization error (b) for NTCDA. Shown in panel b are data for ε = 2, with the same coloring scheme as for the gas phase in Figure 1c. All optimal values of γ and α are given in Table S1.
with increasing dielectric constant (Table S1). The robustness of our tuning procedure is evident from the results of other molecules including pentacene and C60 (see Figure S3 and Table S2). The impact of the α-tuning is felt across individual MO levels. Well-balanced piecewise linearity among all three frontier MO levels is found at the optimal α of 0.32, and the largest deviation is notably reduced, down to about 0.01 eV, less than a third of the maximum for the isolated molecule at the same α (compare Figure 2b with Figure 1c). In Figure 3, we see that, with increasing dielectric constant, the filled MO levels rise without affecting much the gaps between adjacent levels, as they would in actual dielectric environments. This trend continues well into ε = 1 or the gas phase, for which α* from ε = 2 was used with the understanding that α* is nearly independent of the dielectric constant. Indeed, gas-phase calculations with α*(ε = 2) reproduce extremely well the MO ordering of GW (see Figures 1b and 4a for the MOs), a further discussion of which will follow below. For the isolated NTCDA, the only degeneracy down to at least HOMO−8, as predicted by GW, is found at the HOMO− 4 and HOMO−5 levels. If we consider that the HOMO−4 and HOMO−5 curves cross each other only once midway over the entire range of α, our tuning procedure being able to pinpoint the location of degeneracy is simply remarkable (the vertical
Figure 2a shows the new J(γ*, α), as obtained by applying solvation at different dielectric constants (for PTCDA, see Figure S2). A distinct minimum is visible at a nonzero α (with α* denoting the optimal value, α* ∼ 0.28 for PTCDA and α* ∼ 0.32 for NTCDA), in stark contrast to the isolated case for which J(γ*, α) increased monotonically from its minimum at α = 0 (compare Figure 2a with Figure 1a). With an increment of 0.01 for varying α, the IP theorem can be satisfied to an accuracy of
