Article pubs.acs.org/JPCB
Calculating Electron-Transfer Coupling with Density Functional Theory: The Long-Range-Corrected Density Functionals Zhi-Qiang You,† Yi-Chen Hung,† and Chao-Ping Hsu* Institute of Chemistry, Academia Sinica, 128 Section 2 Academia Road, Nankang, Taipei 11529, Taiwan S Supporting Information *
ABSTRACT: The density functional theory (DFT) with commonly used functionals is known to be incorrect for charge-transfer problems. With long-range-corrected (LC) density functionals, the asymptotic exchange potential is gradually switched to the Hartree−Fock exchange at a long range, and the prediction for charge-transfer states is greatly improved. In this work, we test LC-DFT’s performance on charge-transfer couplings. The range-separation parameter can be tuned nonempirically for properties of a generalized DFT. We propose to minimize the difference of highest-occupied Kohn−Sham orbital energy and the ionization potential (for hole transfer) or the lowest-unoccupied orbital energy and the electron affinity (for electron transfer). For photoinduced charge transfer, the minimum in the sum of such differences for the donor and the acceptor is proposed. With the range-separation parameters optimized, we found that ET couplings derived from the LC-DFT are close to those derived from coupled cluster with singles and doubles. When compared with experimentally derived Mulliken−Hush couplings, LC-DFT couplings are greatly improved as well. We also found that the couplings from BNL and LC-BLYP functionals are generally better than those from LC-ωPBE and LC-ωPBE0. LC-DFT is suitable for calculating ET coupling, especially with this nonempirical approach for the range-separation parameter.
1. INTRODUCTION The density functional theory (DFT)1,2 has been a popular approach for solving electronic problems. Various of exchangecorrelation (XC) functionals have been developed, such as those with local density approximation (LDA),2 generalized gradient approximation (GGA),3−5 and meta-GGA functionals involving the kinetic energy density.6,7 Hybrid,8−10 and more recently, double hybrid functionals have also been developed for a wide range of applications.11−14 However, it has been known that most LDA and GGA functionals are with incorrect asymptotic potentials.15−18 Therefore, their ability to predict Rydberg states,16,19,20 charge-transfer state energy,21−23 or other important physical properties is restricted. To fix this problem, Casida proposed a shift-and-splice approach19,20 which is very simple to employ, but the XC potential has a discontinuity. A statistical-averaged overlap potential approach was proposed, which offers a continuous potential based on the orbital population.24−26 Both approaches offer improved XC potentials, not the XC energy functionals, and the potential application may be limited. Until now, the most popular approach to solve this problem has perhaps been the long-range correction in the Coulomb operator,27−29 and a frequently used partition for the range separation is10 1 − erf(ωr12) erf(ωr12) 1 = + r12 r12 r12
than 1/ω. The second term, erf(ωr12)/r12, is for the long-range potential, which becomes the Coulomb potential when r12 ≫ 1/ω. This range-separation approach is the basis of the longrange-corrected (LC) DFT. In LC-DFT, conventional DFT exchange is used with the short-range-modified Coulomb operator, and the full Hartree−Fock (HF) exchange is incorporated with the long-range modified Coulomb operator. The XC energy functionals employed are often GGA or hybrid functionals E XC = EC + (1 − C HF)E X,GGA + C HFE X,HF
where CHF denotes the coefficient of the HF exchange potential in hybrid functionals, which is zero for GGA. The corresponding LC functional is LC SR SR LR E XC = EC + (1 − C HF)E X,GGA + C HFE X,HF + E X,HF
(3)
in which the superscripts “LR” and “SR” denote long-range and short-range components, respectively, where the exchange functional is switched. Correlation potentials are expected to decay more rapidly than exchange potential,15 and they are not changed in the LC scheme. LC-DFT has become a popular alternative to traditional DFTs. It has been shown that LCDFTs yield good charge-transfer state energies,30−33 improved atomization energies, reaction barrier heights, electron affinities
(1) Special Issue: John R. Miller and Marshall D. Newton Festschrift
where “erf” denotes the error function. ω is an adjustable parameter that is used to define the switching of short-range and long-range. In eq 1, the first term is for the short-range potential, which decays to zero quickly as r12 becomes larger © XXXX American Chemical Society
(2)
Received: November 9, 2014 Revised: January 19, 2015
A
DOI: 10.1021/jp511216c J. Phys. Chem. B XXXX, XXX, XXX−XXX
Article
The Journal of Physical Chemistry B
process. Under resonance where Ei = Ef, the coupling value equals a half of the energy gap
(EAs), and ionization potentials (IPs) for atoms and molecules.34−37 The electronic coupling factor in the electron-transfer (ET) rate is an intrinsic electronic property. Because it is an important parameter determining the ET rate, the computational characterization for the coupling has gained much attention recently.38 While much published work used wavefunction-based models like semiempirical INDO or ZINDO,39 Hartree−Fock (HF), configuration interaction singles (CIS),40 or high-level coupled-cluster models,41−44 more and more work reports values of ET coupling calculated with a DFT model.45−48 ET coupling depends heavily on the behavior of the molecular orbitals in the asymptotic regions. The problem at the asymptotic region for DFT is likely associated with the ET coupling values. Because LC-DFT has a full HF exchange when the electrons are away from each other, it is a promising candidate for ET couplings. Therefore, it is necessary to systematically compare both the DFT and LC-DFT coupling values with those derived from wave function based models. In the present work, we test for the performance on ET couplings with a number of density functionals, including LDA exchange, Becke 88 exchange4 and Perdew−Burke−Ernzerhof (PBE) exchange.5 The LC functionals we test are Baer− Neuhauser−Livshits (BNL) functional,49,50 LC-BLYP developed by Hirao’s research group,28,30 as well as LC-ωPBE and LC-ωPBE0,37 which are derived from corrected exchange-hole model by Henderson et al.29 In Section 2, we review the methods for calculating ET coupling with DFT. We then introduce a nonempirical tuning procedure for the rangeseparation parameter. In Sections 3 and 4, we evaluate the performance of LC-DFT on characterizing ET coupling. A brief conclusion can be found in Section 5.
|V | = |E1 − E2| /2
Several other strategies for calculating the ET coupling are discussed in the literature.38,40,52−54 In the following sections, we review some of these approaches from the perspectives of DFT. 2.1. Charge Transfer in the Ground State. Charge transfer in the ground state is often seen in systems such as charge transfer in organic thin-film device, or most redox reactions. As an example, the hole transfer (HT) in a pair of molecules can be written as
D+A → DA+ When a negative charge transferred, also called ET, it is
⎛ −(ΔG 0 + λ)2 ⎞ 2π 2 1 |V | exp⎜ ⎟ ℏ 4λkBT ⎝ ⎠ 4πλkBT
+
) for the cationic dimer (Figure 1). first-excited state (E{DA} ex
Figure 1. Schematic energy diagram for HT (left) and ET (right) couplings. Shown are energy levels at the avoided-crossing region, where the energy gap of the lowest two cationic (anionic) states are twice HT (ET) coupling.
They can be equivalent to the first and second vertical IPs from the neutral dimer (Figure 1). Thus, the HT coupling is calculated by
(4)
|V HT| = |IP I − IP II| /2
where the ΔG is the standard free-energy difference, λ is the reorganization energy, and V is the electronic coupling, that is, the off-diagonal matrix element for the Hamiltonian of the initial and final diabatic states. V can be calculated with a number of schemes.38,52 In these calculations, the initial (final) diabatic state is modeled as states with the electron being transferred localized on donor (acceptor). In general, the diabatic states are not the eigenstates of the system. The Hamiltonian can be written in the two different representations 0
⎛ Ei V ⎞ ⎟ ⎜ ⎝ V Ef ⎠
⇔ ⎛ E1 0 ⎞ ⎜ ⎟ ⎝ 0 E2 ⎠
+
IP I ≡ Eg{DA} − Eg{DA}
E1,2
(11)
and +
IP II ≡ Eex{DA} − Eg{DA}
(12)
Similarly, for an ET process, the ET coupling under a resonance condition can be obtained from the first and second EAs of the neutral dimer (Figure 1)
|V ET| = |EAI − EAII| /2 (5)
(13)
where
with the eigenvalues E + Ef = i ± 2
(10)
where
eigenstates
diabatic states
(8)
(9) D−A → DA− Because they do not involve photoexcitation, the initial and final states are typically modeled as ground states. For the HT process, E1 and E2 in the energy-gap model (eq + 7) are adiabatic energies of the ground state (E{DA} ) and the g
2. THEORY AND METHOD ET is a process where an electron on a donor is transferred to an acceptor. Under the weak-coupling limit, the ET reaction rate is given by the Marcus’ expression51 ket =
(7)
−
⎛ E i − E f ⎞2 2 ⎜ ⎟ + V ⎝ 2 ⎠
EAI ≡ Eg{DA} − Eg{DA}
(14)
and (6)
−
EAII ≡ Eg{DA} − Eex{DA}
where we have assumed that the initial and final states are orthogonal and normalized. In eqs 5 and 6, Ei (Ef) denotes the energy of the initial (final) state and V is the coupling for an ET
(15)
Instead of directly calculating the state energies, the differences in IPs or EAs can be used for calculating B
DOI: 10.1021/jp511216c J. Phys. Chem. B XXXX, XXX, XXX−XXX
Article
The Journal of Physical Chemistry B coupling.43,44 In practice, this IP/EA scheme is easier to calculate than directly calculating the unrestricted cationic or anionic state energies. A spin-flip approach was also proposed for calculating the energy gaps.42,44 Simplified one-electron model can also be used, which is outlined later. In molecular orbital theories, such as the Hartree−Fock (HF) and the Kohn−Sham (KS) DFT, the relationship of orbital energies and IP/EA has been established. In closed-shell Hartree−Fock theory, the Koopmans’ theorem (KT)55 links the negative occupied HF orbital energies and the vertical ionization energies. It also links the negative virtual HF orbital energies as EAs. In both cases, the orbital relaxation upon ionization is neglected. With IPs and EAs in eqs 10 and 13 replaced by orbital energies, both HT and ET couplings can be calculated in terms of the HF orbitals. This is the HF-KT coupling scheme for charge-transfer couplings.42,53,56 For the exact KS-DFT, the negative KS-HOMO energy (ϵKS H) has been shown to be the first IP of a molecular system57−60 −ϵHKS(N ) = IP I(N ) ≡ Eg (N − 1) − Eg (N )
experimental parameters.70−72 With a perturbation theory to model the weak transition strength from the charge-transfer state, The MH coupling is obtained72 |V | = 2.06 × 10−2
ϵKS H
(16)
V GMH =
V FCD =
Δqmn =
(18)
where and denote the LUMO and LUMO+1 energies of KS orbitals, respectively. 2.2. ET Involving Excited States. ET may involve excited states of the system, typically accompanied by absorption of a photon. For example, the photoinduced ET is hv
hv
DA → D+A−
(Δμmn )2 + 4(μmn )2
(22)
|Δqmn|ΔEmn (Δqmm − Δqnn)2 + 4(Δqmn)2
(23)
∫r∈D ρmn (r) dr − ∫r∈A ρmn (r) dr
(24)
and ρmn(r) is the one-particle density matrix. It has been known that LDA or GGA functionals are problematic for charge-transfer excitations. GMH and FCD calculations depend on properties of the charge-transfer states, such as excitation energies, dipole moments, or charge densities. LDA or GGA typically underestimate charge-transfer excitation energies. Moreover, because charge-transfer states are not correctly described, other properties such as the dipole moment or the FCD may not be qualitatively correct. It was demonstrated that LC-DFTs yield improved charge-transfer excitations.30−32 It would be of great interest to see whether GMH or FCD calculations with LC-DFT can yield good coupling values. 2.3. Molecules Studied. In the present work, we study the ground-state charge transfer in three π-containing molecules, ethylene, pyrrole, and naphthalene, as shown in Figure 2. The first two molecules are known as basic unit of polymer semiconductors.74−76 The couplings are calculated with two identical molecules stacked in parallel arrangement, as depicted in Figure 2B. The optimized structures of neutral or ionic monomers were first calculated using DFT with B3LYP hybrid functional4,77,78 and Dunning’s double-ζ basis sets (DZP),79 except for pyrrole, where 6-311+G* basis set was used. We verified the structures
(17)
ϵKS L+1
DA → D*A → D+A− or the optical charge transfer
|μmn |ΔEmn
where Δqmn are FCD matrix elements, defined as
where and denote the HOMO and HOMO−1 energies of KS orbitals, respectively. The KT analogue for unoccupied KS orbitals has also been established by theoretical and numerical approaches.65−69 It was found that LC-DFT calculations for atoms and small molecules yields LUMO energies close to the first EA obtained from coupled-cluster calculations or experiments.65,69 Therefore, we can approximate the ET coupling in KS-DFT by ϵKS L
(21)
where μmn is the transition dipole moment, ΔEmn is the energy difference, and Δμmn is the permanent dipole moment difference, all of which are calculated between the initial and final eigenstates |m⟩ and |n⟩. A similar scheme, the fragment charge difference (FCD), was later developed.54 In FCD, the dipole moment operator in GMH is replaced by a difference charge operator, which reflects the extent of charge localization. The FCD coupling is given by
ϵKS H−1
|V ET| ≈ |ϵLKS−ϵLKS + 1| /2
RDA
where εmax is the molar extinction coefficient at maximum absorption, in the units of M−1 cm−1, and Δν1/2 is the full width at half-maximum (fwhm) of charge-transfer band. In eq 21, V, νmax, and Δν1/2 are in units of wavenumber, cm−1. Equation 21 offers a way to calculate charge-transfer coupling from characteristic values of the charge-transfer band in absorption spectra, and it is an experimental coupling value. A generalization of the MH expression (eq 21) was developed from very similar assumptions.73 With this generalized Mulliken−Hush (GMH) method, the coupling matrix element can be calculated directly from excited-state quantum calculations
where Eg(N) is the ground-state energy for an N-electron system. We note that it is a general property of DFT, and no approximation is involved, as opposed to the KT for HF. However, it has been known that commonly used LDA and GGA functionals do not yield correct orbital energies and band gaps due to the self-interaction error and incorrect XC potential in asymptotic region.15−18 When an XC functional is fixed with a proper asymptotic potential, the validity of eq 16 can be observed. With an XC potential derived from a statistical average of model orbital potentials, Chong et al. showed that KS eigenvalues are very close to vertical ionization energies.61 More theoretical and numerical studies confirmed the physical implication of the occupied KS orbital energies.62,63 LC-DFT provides a useful route to correct asymptotic potentials. Recent works64,65 showed that LC-DFT can restore the correct KS orbital energies close to the corresponding IPs for several small molecules. Therefore, in KS-DFT with correct asymptotic XC potential, the HT coupling can be calculated as |V HT| ≈ |ϵHKS−ϵHKS− 1| /2
νmaxεmax Δν1/2
(19)
(20)
For excited-state ET, the Mulliken−Hush (MH) expression provides an estimation of charge transfer coupling from C
DOI: 10.1021/jp511216c J. Phys. Chem. B XXXX, XXX, XXX−XXX
Article
The Journal of Physical Chemistry B
charge-transfer bands in the absorption spectra.83,84 In one series, molecules are with (1,4)-dimethoxybenzene as electron donor and (1,4)-benzoquinone as acceptor (1). Another pair of molecules (2) has the same electron donor but with (7,7,8,8)tetracyanoquinodimethane as the electron acceptor. In each series, we studied both pseudogeminal and pseudo-ortho structures, which are named as a and b, respectively. The ground-state structures were optimized using B3LYP/6311G*. The couplings and charge-transfer excitations were calculated using the single excitation theory, HF-CIS, and TDDFT, with the basis set cc-pVDZ. Because the experimental charge-transfer spectra were obtained in chloroform, we performed all calculations with the polarizable continuum model (PCM). All computation was performed with a developmental version of Q-Chem.85 2.4. Tuning the Range-Separation Parameter. The optimal value of the range-separation parameter ω in eq 1 varies with different XC functionals and different molecular systems studied.49,50,86−89 Using values suggested in the literatures may not be sufficient Because they were from benchmark results across many different molecular systems, and in most of the cases the best-fit property was excitation energy. In general, tuning ω value in LC-DFT is still required. One general property we aimed is that the negative KS-HOMO energy equals the vertical ionization energy.57−60 Thus, we propose to choose to the parameter that minimizes the KS . Following a simple difference between IP and −ϵHOMO 32,50 we propose to tune ω by minimizing optimization scheme, the error in IP or EA defined as
Figure 2. (A) Molecules tested for the present work. From left to right: ethylene, pyrrole, and naphthalene. (B) Stacked dimer, with naphthalene as an example.
by a positive definite Hessian matrix. For simple test systems composed of small molecules, an approximated reaction coordinate, R, is often used in the literature42,80,81 Q(R ) = (1 − R )Qr + R Qp, 0 ≤ R ≤ 1
(25)
where Q represents the nuclear coordinates of the system, Qr and Qp are the initial (reactant) and final (product) state nuclear coordinates that are composed of corresponding optimized neutral and ionic molecules. R is the reaction coordinate, with R = 0 for the initial state and R = 1 for the final state. For symmetric dimers in the present study, R = 0.5 was used as an approximated transition state. We also tested the performance of LC-DFT for excited-state ET with two series of molecular complexes. The first series is composed of a tetracyanoethylene (TCNE) electron acceptor and an aromatic donor, which includes benzene, toluene, oxylene, or naphthalene (Figure 3A) in the present work. The
KS |ΔEIP(ω)| ≡ |− ϵHOMO,neutral (ω) − IP(ω)| KS KS KS = |− ϵHOMO,neutral (ω) − [Ecation (ω) − Eneutral (ω)]|
(26) KS |ΔE EA (ω)| ≡ |− ϵLUMO,neutral (ω) − EA(ω)| KS KS KS = |− ϵLUMO,neutral (ω) − [Eneutral (ω) − Eanion (ω)]|
(27) KS
ϵKS HOMO
where E (ω), the KS ground-state energy, and (ω), the KS-HOMO eigenvalue, are calculated using the LC functional with the parameter ω. The subscripts “neutral”, “cation”, and “anion” indicate that the value was calculated in a neutral, cationic, and anionic monomer, respectively. For EA in eq 27, we also tested with the energy of the singly occupied MO energy of the radical anion, instead of the LUMO of the neutral molecule, and the optimized parameters were the same. Minimizing the errors in IP (eq 26) or in EA (eq 27) is likely to yield an optimal condition for the ground-state HT and ET coupling calculations. For excited-state ET, the charge transfer involves an electron removal from donor and an electron attachment to acceptor. Thus, we propose to find the optimal parameter ω for excited-state ET problems with
Figure 3. Systems studied for excited electron transfer: (A) chargetransfer complexes formed by a tetracyanoethylene (TCNE) acceptor and an aromatic donor, benzene, toluene, o-xylene, and naphthalene. (B) [3,3]Paracyclophane molecules.
min|ΔE DIP(ω)| + |ΔEAEA (ω)| ω
structures were obtained from ref 32. The charge-transfer excitations and the FED couplings were calculated using TDDFT. With the MH expression (eq 21), it is possible to derive empirical coupling values and compare them with our calculated results.82 The second series we studied are four different [3,3]paracyclophanes83,84 (Figure 3B). These are molecules with
(28)
The parameters that produce the best IP (eq 26) and EA (eq 27) are listed in Table 1. It is seen that the parameter systematically decreases as the size of the molecule increases. We note that the parameters for good IP and EA are close, and this result implies a consistency between the two conditions eqs 26 and 27. D
DOI: 10.1021/jp511216c J. Phys. Chem. B XXXX, XXX, XXX−XXX
Article
The Journal of Physical Chemistry B Table 1. Optimal Range-Separation Parameter Values IPa
EAb
functional
ethylene
pyrrole
naphthalene
ethylene
pyrrole
naphthalene
BNL LC-BLYP LC-ωPBE LC-ωPBE0
0.41c 0.39 0.36 0.31
0.33 0.31 0.30 0.25
0.28 0.26 0.26 0.21
0.43 0.41 0.37 0.32
0.35 0.33 0.31 0.26
0.29 0.28 0.27 0.22
a Parameters that minimize ΔE as in eq 26, calculated using aug-cc-pVDZ basis set. bParameters that minimize ΔE as in eq 27, calculated using DZP basis set. cIn the units of bohr−1.
Table 2. Range-Separation Parameters for Excited-State Electron-Transfer Problems TCNE complexes
a
cyclophanes
functional
benzene
toluene
o-xylene
naphthalene
1a
1b
2a
2b
BNL LC-BLYP LC-ωPBE
0.33a 0.32 0.30
0.32 0.30 0.29
0.31 0.28 0.27
0.32 0.31 0.30
0.30 0.27 0.26
0.31 0.28 0.26
0.27 0.26 0.25
0.27 0.26 0.25
Parameters obtained following eq 28, calculated with cc-PVDZ basis set, in units of bohr−1.
derived from EOM-CCSD. The mean absolute errors in IP and EA, across the molecules studied, are also calculated
In the tests for excited-state charge-transfer couplings, the parameters were obtained with the expression in eq 28, and the results are listed in Table 2. 2.5. Other Computational Details. The ground-state ET coupling values were calculated using the KT-IP(EA) schemes, eq 17and 18, with different density functional approximations. The distance dependence of the couplings is defined as V = V0 exp( −βd)
MAE (%) of IP =
(29)
100 × N
DEV (%) of β = 100 ×
N
∑ i=1
MAE (%) of EA =
β DFT − β CCSD β CCSD
m=1
100 × NM
NM
KS (m) − EACCSD(m)| ∑ |−ϵLUMO m=1
(33)
where m is the index for different molecular systems tested, and NM is the total number of the systems studied (which is 3 in this case). For excited-state charge transfer, an error value was derived for every molecule in each DFT or LC-DFT calculation error (%) = 100 ×
|VFCD| − |VMH| |VMH|
(34)
where VMH is the MH empirical coupling and VFCD is the FCD coupling value.
3. GROUND-STATE ELECTRON TRANSFER 3.1. HT and ET Coupling Values. The ground HT and ET couplings were calculated for the XC functionals tested. The difference in coupling values, as compared with CCSD results, were quantified as MAEs (eqs 30−33), and they are included in Figures 4 and 5. In Figure 4, it is seen that the MAE for IPs are reduced significantly in LC-DFTS, from well above 20% in traditional DFTs to ∼1% in LC-DFTs. The errors for the HT couplings (V) and the exponential decay rates (β) are also improved with LC-DFTs. The overall errors from LDA and GGA functionals are ∼35%, and the errors from hybrid-GGA functionals, B3LYP and PBE0, are slightly better, ∼30%. With LC functionals, the resulting MAEs for coupling are significantly reduced. The mean errors from BNL and LC-BLYP functionals are
