Article pubs.acs.org/JCTC
Accurate Ionization Potentials and Electron Affinities of Acceptor Molecules I. Reference Data at the CCSD(T) Complete Basis Set Limit Ryan M. Richard,† Michael S. Marshall,† O. Dolgounitcheva,‡ J. V. Ortiz,‡ Jean-Luc Brédas,§ Noa Marom,∥ and C. David Sherrill*,† †
Center for Computational Molecular Science and Technology, School of Chemistry and Biochemistry, and School of Computational Science and Engineering, Georgia Institute of Technology, Atlanta, Georgia 30332-0400, United States ‡ Department of Chemistry and Biochemistry, Auburn University, Auburn, Alabama 36849-5312, United States § Solar & Photovoltaics Engineering Research Center, Physical Science and Engineering Division King Abdullah University of Science and Technology, Thuwal 23955-6900, Kingdom of Saudi Arabia ∥ Department of Physics, Tulane University, New Orleans, Louisiana 70118-5645, United States S Supporting Information *
ABSTRACT: In designing organic materials for electronics applications, particularly for organic photovoltaics (OPV), the ionization potential (IP) of the donor and the electron affinity (EA) of the acceptor play key roles. This makes OPV design an appealing application for computational chemistry since IPs and EAs are readily calculable from most electronic structure methods. Unfortunately reliable, high-accuracy wave function methods, such as coupled cluster theory with single, double, and perturbative triples [CCSD(T)] in the complete basis set (CBS) limit are too expensive for routine applications to this problem for any but the smallest of systems. One solution is to calibrate approximate, less computationally expensive methods against a database of high-accuracy IP/EA values; however, to our knowledge, no such database exists for systems related to OPV design. The present work is the first of a multipart study whose overarching goal is to determine which computational methods can be used to reliably compute IPs and EAs of electron acceptors. This part introduces a database of 24 known organic electron acceptors and provides high-accuracy vertical IP and EA values expected to be within ±0.03 eV of the true non-relativistic, vertical CCSD(T)/CBS limit. Convergence of IP and EA values toward the CBS limit is studied systematically for the Hartree− Fock, MP2 correlation, and beyond-MP2 coupled cluster contributions to the focal point estimates.
1. INTRODUCTION The sun bombards the United States with about 13 PJ (P = peta) of energy every second.1 This means that about 15 seconds of sunlight could power the United States for an entire day,2 making the sun an appealing “renewable” energy source. Unfortunately, in its native form as electromagnetic radiation, the light from the sun cannot be directly used by most modern electronics and must first be interconverted into usable electricity. The large family of devices charged with this task are known as photovoltaic cells (colloquially “solar cells”), and much scientific research into creating new generations of solar cells and improving existing ones is underway.3 Presently, we limit ourselves to organic photovoltaics (OPVs, see refs 4 and 5 for recent reviews), whose composition in terms of organic molecules makes them both potentially cost-effective to manufacture and computationally tractable on a molecular basis. Although we are motivating the present study from the perspective of OPVs, we anticipate the results of the current study to be more broadly applicable to organic electronic design in general, particularly any time chargetransfer complexes play a key role in that design. The active layer in an organic solar cell is made of an electrondonor component and an electron-acceptor component. It has © 2016 American Chemical Society
been shown that the voltage afforded by the cell is directly related to the difference between the ionization potential of the donor and the electron affinity of the acceptor.4 The ionization potential (IP), EIP, of a system is given by
E IP = E+ − E0
(1)
where E+ and E0 are the energies of the cationic and neutral states, respectively. The IP may be interpreted as the amount of work required to remove (add) an electron (hole). Similarly, the electron affinity (EA), EEA, is given by E EA = E0 − E−
(2)
where E− is the energy of the anion, and the EA may be interpreted as the negative of the amount of work required to add (remove) an electron (hole). Here, we examine IP and EA values of individual molecules as a prerequisite first step, although we note that in actual solar cell design one would also need to consider the effect of the solid state on the IP and EA values. Received: September 11, 2015 Published: January 5, 2016 595
DOI: 10.1021/acs.jctc.5b00875 J. Chem. Theory Comput. 2016, 12, 595−604
Article
Journal of Chemical Theory and Computation
Figure 1. The 24 molecules in our IP/EA database, as well as the abbreviations we use throughout the manuscript, if applicable. Molecules with blue backgrounds are included in our calibration set.
def2-TZVPP basis), considers molecules of more interest in the context of organic electronics, and also considers EAs in addition to IPs. Here, we report non-relativistic vertical IPs and EAs for a series of 24 small organic molecules; the structures, names, and relevant abbreviations are given in Figure 1. We note that the current manuscript is the first part of a four part study whose overarching goal is to compare different electronic structure methods’ abilities to reproduce a common reference set of benchmark-quality IP and EA values; the remaining three parts focus on the performance of non-empirically tuned, longrange corrected DFT methods,34 GW methods,35 and electron propagator methods,36 respectively.
The molecular EAs and IPs can be readily calculated using wave function or density functional theory (DFT) methods. The question then becomes a matter of choosing a level of electronic structure theory that can accurately predict EAs and IPs for a broad class of relevant systems. The de facto standard in high-accuracy electronic structure theory is coupled cluster theory with single, double, and perturbative triples [CCSD(T)],6 extrapolated to the complete basis set (CBS) limit. Unfortunately, the computational complexity of conventional CCSD(T) scales as the seventh power of the system size, and the method is therefore only accessible to the smallest of systems. Ideally, one would like to find a computationally cheaper alternative to CCSD(T)/CBS without sacrificing accuracy; however, this requires a suitably accurate database of benchmark values to calibrate against. Surprisingly, there seems to be a lack of high-level theoretical IP and EA data for moderately sized organic molecules relevant to OPV design. For atoms or small molecules (about five heavy atoms or less), many authors have undertaken high-level computations on a few molecules,7−27 and experimental IP/EA databases do exist.28−30 However, experimental uncertainties can be significant [e.g., many of the tabulated experimental EAs in ref 31 have errors around 0.05 eV or more; also see ref 32 for a discussion on uncertainties of IPs and EAs in the National Institue of Standards and Technology’s (NIST) Webbook where uncertainties may be as high as 1.0 eV]. In addition, comparing directly to experiments introduces additional degrees of freedom such as differences in geometry, vibrational effects, etc. A systematic study comparing approximate electronic structure methods for IP/EA values is more easily performed by computing differences in electronic energies at fixed geometries (“vertical” values) and is the approach adopted here. As we were preparing the present manuscript for submission, another study presenting high-quality ionization potentials was published by Krause et al.,33 who used the CCSD(T)/def2-TZVPP level of theory. The present study is complementary in that it provides higher-quality data (CBS extrapolation instead of the fixed
2. THEORY The benchmark IP and EA values presented here are the combination of two commonly used computational chemistry techniques: focal point analysis (FPA)37 and CBS extrapolation. Both of these techniques are designed to use combinations of lower cost theories to approximate a higher level of theory. The technical details are given in Section 3, and the current section instead focuses on summaries of the techniques. In the interest of brevity, we will henceforth utilize the common abbreviation aXZ to refer to the aug-cc-pVXZ38−41 (X = D, T, Q, etc.) basis sets. 2.1. CBS Extrapolation. Quantum chemistry computations suffer from an error brought about by the basis set being insufficiently flexible, known as basis set incompleteness error (BSIE). The correlation consistent family of basis sets, of which the aXZ series are members, is specifically designed to add basis functions in a systematic way to facilitate smooth convergence to the CBS limit.38 Hence, it is possible to remove some of the BSIE by extrapolation. At some superficial level, most of the existing CBS extrapolations make recourse to work emanating from Schwartz’s42 perturbative expansions of the leading-order term to the correlation energy (CE) for a two-electron atom. Inherent in this work is 596
DOI: 10.1021/acs.jctc.5b00875 J. Chem. Theory Comput. 2016, 12, 595−604
Article
Journal of Chemical Theory and Computation
2.2. Focal Point Analysis (FPA). CBS extrapolations are a way to estimate a given quantity as if it had been computed in a larger basis set. FPAs37 are a related technique with the same end goal; however, instead of using extrapolation techniques, FPAs assume that different parts of the CE converge to the CBS limit faster than other parts. Consequently, the faster a quantity converges to the CBS limit, the smaller the basis set that is required in order to compute that term accurately. To the extent that higher-order correlation effects correspond to faster converging terms, this has the potential to save a substantial amount of computational resources. In fact several studies55,56 have shown that this is indeed the case for many of the components of the CE. Perhaps one of the more common applications of FPA, and the one used presently, is to obtain the CCSD(T) correction to MP2 more efficiently. Starting with the CCSD(T) energy in a “large” basis set, Elarge CCSD(T), we may exactly write
the assumption that the wave function is represented via a linear combination of all angular momentum eigenfunctions up to some angular momentum S . Given this, Schwartz showed that lim ES − ES − 1 ∝ S −4
S →∞
(3)
where ES is the energy computed using eigenfunctions associated with angular momenta up to S [the usually cited 1 result of (S + 2 )−4 is a more exact, intermediate step that is not explicitly shown in ref 42]. In the spirit of Schwartz’s original study, Kutzelnigg et al.43,44 demonstrated, among other things, that the same limiting behavior is encountered for the CE obtained using Møller−Plesset second-order perturbation theory (MP2) for closed-shell atoms. Unfortunately, the assumptions of Schwartz and Kutzelnigg’s studies prohibit their results from being directly applied to typical computational chemistry problems. As is commonly pointed out,45 the typical Gaussian basis sets used in computational chemistry are not saturated with respect to angular momentum, nor are molecular wave functions eigenfunctions of angular momentum. Modern endeavors46,47 to understand how these situations generalize to more typical situations have been undertaken, and the actual limiting behavior is considerably more complicated and poorly understood. Consequently, current CBS extrapolation procedures are perhaps better thought of as empirical corrections motivated by, rather than derived from, the underlying physics. This leaves us with some ambiguity in deciding which CBS extrapolation to use. A large number of CBS extrapolations have been suggested (for comprehensive analyses of the various forms, see refs 45 and 48). In thorough studies of CBS extrapolations,45,49 no given CBS extrapolation has been found to be the “best” in all cases. Fortunately, even in the worst case scenarios, Feller et al.45,49,50 conclude that an extrapolation using at most an aXZ basis set produces results that are comparable to the unextrapolated a(X+1)Z basis set results. In favorable cases, the results are typically comparable to unextrapolated a(X+2)Z basis set results. Consequentially, it appears that any CBS extrapolation can be expected to produce higher quality results than the unextrapolated results alone (note, however, that numerous studies, e.g., refs 51 and ref 52, suggest that aDZ results may be so poor as to prevent meaningful extrapolation). Given the lack of consensus on the “best” CBS extrapolation procedure, arguably the most popular CBS extrapolation procedure (and the one used presently) is the one proposed by Helgaker et al.51,53 Within this scheme the Hartee−Fock (HF) energy and the CE are separately fit to different expressions. The HF energies using the three basis sets of highest cardinal number may be fit to the three-parameter, decaying exponential of the form originally proposed by Feller:54 X CBS E HF = E HF + α e − βX
α X3
(6)
large CCSD(T),large large δMP2 ≡ ECCSD(T) − EMP2
(7)
Of course as written this represents no cost savings, but given the body of empirical evidence55−59 suggesting that δCCSD(T) MP2 CCSD(T) converges faster than Elarge in MP2, we can instead compute δMP2 a “small” basis set: large CCSD(T),small large ECCSD(T) ≈ δMP2 + EMP2
(8)
56,60,61
Previous work has shown that in some cases, especially when δCCSD(T) is small in magnitude, a double-ζ basis set can MP2 lead to the wrong sign for this correction, so that the “small” basis set should be of triple-ζ quality, or better if possible. Work by our group60 has shown that although the δCCSD(T) converges MP2 faster than either EMP2 or ECCSD(T), there can remain some small variations as the basis set is increased toward the CBS limit. Discussions in Section 4.2 of the current manuscript further support these claims.
3. COMPUTATIONAL METHODS The geometries shown in Figure 1 have been optimized within the Gaussian 09 electronic structure package,62 at the B3LYP/ 6-311G** level, while being constrained to the highest possible point-group symmetry. Cartesian coordinates are available in the Supporting Information. The benchmark quality IPs and EAs presented here are the result of CBS extrapolations and δCCSD(T) corrections. Except where noted, all additional comMP2 putations were performed with MOLPRO version 2010.1,63 and employ canonical restricted CCSD(T)64 and MP2, based on restricted (restricted open-shell) references for neutral (ionic) molecules. Energy convergence is defined as 10−9 a.u., and the electron integrals are neglected below 10−14 a.u. All of the CCSD(T) and MP2 computations neglect core correlation (i.e., we have invoked the frozen-core approximation), the errors from which are expected to be small and are discussed below. For the remainder of the manuscript, the level of theory used to compute a particular value is denoted as MP2/X+δCCSD(T) /Y, MP2 where X and Y refer to basis sets, and one or both may be a CBS extrapolation, which is denoted as a(XY)Z. The HF energy is computed in the same basis as the largest MP2 calculation; extrapolation of the HF energy produced only nominally better results, vide inf ra, and was consequently neglected in favor of pure values.
(4)
where ECBS HF , α, and β are fitting parameters, and X is the cardinal number of the basis set (e.g., X = 2 for aDZ, X = 3 for aTZ, etc.). The CE computed in the basis sets with the two highest cardinal numbers are fit using a two-parameter inverse cubic function of the form: X CBS ECE = ECE +
large CCSD(T),large large ECCSD(T) = δMP2 + EMP2
(5)
where α and ECBS CE are fitting parameters. 597
DOI: 10.1021/acs.jctc.5b00875 J. Chem. Theory Comput. 2016, 12, 595−604
Article
Journal of Chemical Theory and Computation
Table 1. Experimental (taken from ref 70, see notes in text) and Theoretical Vertical Ionization Potentials (IP) and Vertical Electron Affinities (EA), in eVa symmetry molecule name
experimental
theoretical
cation
anion
IP
EA
IP
EA
level of theory
anthracene
B2g
B3u
7.44
0.53†
0.33 0.28c
MP2/a(TQ)Z+δCCSD(T) /a(DT)Z MP2
acridine phenazine azulene
B1 B2g A2
B1 B3u B1
7.8 8.44 7.42
0.90† 1.31 0.80†
7.52 7.47b 7.01d 8.04 8.47 7.55 7.58e
0.69 1.11 0.54 0.40e
MP2/a(TQ)Z+δCCSD(T) /aDZ MP2 MP2/a(TQ)Z+δCCSD(T) /a(DT)Z MP2 MP2/a(TQ)Z+δCCSD(T) /a(DT)Z MP2
benzoquinone
B3g
B2g
10.0
1.85†
B2 B1 B1g B1g
A2 A2 B2g B2g
9.5 9.5 10.7† 9.74†
1.81 2.21 2.70 2.78
1.55 1.76g 1.71h 1.47 1.92 2.29 2.48
MP2/a(Q5)Z+δCCSD(T) /a(TQ)Z MP2
naphthalenedione dichlone F4-benzoquinone Cl4-benzoquinone
10.27 10.05f N/A 9.88 9.99 11.14 10.25
nitrobenzene
A2
B1
9.94
1.00†
F4-benzenedicarbonitrile dinitrobenzonitrile nitrobenzonitrile benzonitrile fumaronitrile mDCNB TCNE
B1g B1 B1 B1 Au A2 B3u
B3u A2 B1 B1 Bg A2 B2g
10.65* N/A 10.59 9.73 11.3* 10.20 11.79†
1.89 2.16 1.69 0.26 1.25 0.91 3.16†
TCNQ
B3u
B2g
N/A
2.80
10.19 N/A 10.76 11.15 10.62 9.93 11.48 10.45 11.99 N/A 9.57 N/A
0.54 0.32i 1.62 1.76 1.30 −0.21 0.98 0.61 3.05 2.94j 3.33 3.22h
anhydrides
maleic anhydride phthalimide phthalic anhydride Cl4-isobenzofurandione NDCA
B2 B1 A2 A2 A2
A2 A2 A2 A2 A2
11.07* 9.90* 10.1† 10.8 8.92
1.44 1.02 1.25 1.96 N/A
11.33 10.08 10.55 10.05 9.14
1.01 0.63 0.87 1.68 1.26
other
boron-dipyrromethene
A″
A′
N/A
N/A
8.07
1.67
acenes
quinones
nitro/nitriles
MP2/a(TQ)Z+δCCSD(T) /a(DT)Z MP2 MP2/a(TQ)Z+δCCSD(T) /a(DT)Z MP2 MP2/a(TQ)Z+δCCSD(T) /a(DT)Z MP2 MP2/a(TQ)Z+δCCSD(T) /a(DT)Z MP2 MP2/a(TQ)Z+δCCSD(T) /a(DT)Z MP2 MP2/a(TQ)Z+δCCSD(T) /a(DT)Z MP2 MP2/a(TQ)Z+δCCSD(T) /a(DT)Z MP2 MP2/a(TQ)Z+δCCSD(T) /a(DT)Z MP2 MP2/a(Q5)Z+δCCSD(T) /a(TQ)Z MP2 MP2/a(Q5)Z+δCCSD(T) /a(TQ)Z MP2 MP2/a(TQ)Z+δCCSD(T) /a(DT)Z MP2 MP2/a(Q5)Z+δCCSD(T) /a(TQ)Z MP2 MP2/a(TQ)Z+δCCSD(T) /a(DT)Z MP2
MP2/a(Q5)Z+δCCSD(T) /a(TQ)Z MP2 MP2/a(TQ)Z+δCCSD(T) /a(DT)Z MP2 MP2/a(TQ)Z+δCCSD(T) /a(DT)Z MP2 MP2/a(TQ)Z+δCCSD(T) /aDZ MP2 MP2/a(TQ)Z+δCCSD(T) /aDZ MP2 /aDZ MP2/a(TQ)Z+δCCSD(T) MP2 †
a
We have denoted experimental values as either adiabatic (†) or vertical (*) when such information was available. Experimental value is adiabatic. CCSD(T) Experimental value is vertical. bMP2/cc-pV(DTQ)Z + δMP4 /cc-pVDZ on Bg (sic) state from ref 11. cMP2/cc-pV(TQ5)Z + MP2/cc-pVTZ + δMP4 d CCSD(T) /cc-pVTZ + δ /cc-pVDZ on B state from ref 12. CCSD(T)/def2-TZVPP on unknown state from ref 33. eMP2/cc-pV(DTQ)Z + δMP4 MP2 MP4 3u f MP4 CCSD(T) δMP2/cc-pVTZ + δMP4 /cc-pVDZ on B2 state from ref 14. EOM-IP-CCSD/DZP+CC3 on B3g state (see ref 24 for full details). gAdiabatic CCSD(T)/aDZ on Ag (sic) state (see ref 25 for full details). hAdiabatic CCSD(T)/aDZ (no diffuse functions on hydrogen) on unknown state (see ref 23 for full details). iCCSD/D95+* on B1 state (see ref 26 for full details). jAdiabatic CCSD(T)/aDZ on B2g state (see ref 27 for full details). *
4. RESULTS AND DISCUSSION Table 1 lists the 24 molecules in our benchmark set, divided into four families (and boron-dipyrromethane). When available, we have listed the recommended experimental values from the NIST Webbook70 (if a recommended value was not available, the listed value is simply the most recent). The process used by NIST to establish recommended values is rather involved, with careful consideration of many experimental values obtained from widely differing experiments. Experimental results in Table 1 are a mix of adiabatic and vertical values; when the experimental data have been described as either adiabatic or vertical in the primary references, we have labeled them as such in the table. We suggest that the interested reader consult the literature cited for more detailed information. Lastly, we point out that vertical experimental values usually are
Ultimately, we will argue that our least accurate benchmark values lie within about ±0.03 eV of the non-relativistic, vertical CCSD(T)/CBS limit. With regard to this limit, BSIE in our estimates is responsible for about 0.02 eV of the uncertainty, the remainder comes from neglect of core correlation, which has been shown to be usually less than ±0.01 eV.10,65−67 Our own analysis, at the MP2/aug-cc-pCVDZ level of theory (presented in the Supporting Information) confirms that in all cases the corecorrelation is on the order of 10−3 eV; admittedly the double-ζ basis set is not large, but even if these estimates are off by a factor of ∼2, the core correlation contributions would still be negligible. Compared to the ab initio Born−Oppenheimer limit, our values also neglect correlation beyond perturbative triples (errors expected to be
