Accurate Ionization Potentials and Electron ... - ACS Publications

Jan 5, 2016 - Department of Chemistry and Biochemistry, Auburn University, ... Center for Computational Molecular Science and Techology, School of ...
0 downloads 0 Views 2MB Size
Article pubs.acs.org/JCTC

Accurate Ionization Potentials and Electron Affinities of Acceptor Molecules IV: Electron-Propagator Methods O. Dolgounitcheva,† Manuel Díaz-Tinoco,† V. G. Zakrzewski,† Ryan M. Richard,‡ Noa Marom,§ C. David Sherrill,‡ and J. V. Ortiz*,† †

Department of Chemistry and Biochemistry, Auburn University, Auburn, Alabama 36849-5312, United States Center for Computational Molecular Science and Techology, School of Chemistry and Biochemistry and School of Computational Science and Engineering, Georgia Institute of Technology, Atlanta, Georgia 30332-0400, United States § Department of Physics and Engineering Physics, Tulane University, New Orleans, Louisiana 70118-5645, United States ‡

S Supporting Information *

ABSTRACT: Comparison of ab initio electron-propagator predictions of vertical ionization potentials and electron affinities of organic, acceptor molecules with benchmark calculations based on the basis set-extrapolated, coupled cluster single, double, and perturbative triple substitution method has enabled identification of self-energy approximations with mean, unsigned errors between 0.1 and 0.2 eV. Among the self-energy approximations that neglect off-diagonal elements in the canonical, Hartree−Fock orbital basis, the P3 method for electron affinities, and the P3+ method for ionization potentials provide the best combination of accuracy and computational efficiency. For approximations that consider the full self-energy matrix, the NR2 methods offer the best performance. The P3+ and NR2 methods successfully identify the correct symmetry label of the lowest cationic state in two cases, naphthalenedione and benzoquinone, where some other methods fail.



INTRODUCTION

reference orbitals, GW approximations, and second-order screened exchange improvements have been compared. Ab initio electron-propagator (EP) methods for determining IPs and EAs 7−16 have been used extensively in the interpretation of photoelectron spectra of molecules and molecular ions.17−21 A given EP approximation is defined by the form of the self-energy (∑) employed in the Dyson equation. The most widely employed approximations are based on canonical Hartree−Fock (HF) orbitals and their energies, which are used to define the zeroth order EP. Perturbative treatments of the self-energy matrix yield electron binding energies and Dyson orbitals (see below) that incorporate two kinds of corrections that are neglected at the level of Koopmans’s theorem (KT): final-state orbital relaxation and electron correlation. Nondiagonal elements of the self-energy matrix in the canonical HF basis are neglected in the diagonal self-energy, or quasiparticle, approximations. These methods are applicable only for IPs and EAs where the Koopmans approximation is qualitatively valid. Greater flexibility accompanies the use of nondiagonal self-energy methods because Dyson orbitals may be expressed therein as linear combinations of HF orbitals. Nondiagonal approximations that include certain classes of self-energy terms in all orders of the

Procuring work from solar radiation is often a matter of charge separation in molecular assemblies.1−3 Design of cheap, efficacious, photovoltaic devices may be abetted by accurate predictions of the ionization potentials (IPs) and electron affinities (EAs) of organic molecular constituents. The subtleties of electron-transfer kinetics demand accurate predictions because chemical modifications may have small but technologically significant effects on electron binding energies. Systematic optimization of these modifications requires that predictions be made on many molecules with a high degree of computational efficiency. To provide reliable, well-defined standards for comparison of various approaches to calculating vertical electron binding energies, basis set extrapolated, coupled cluster single, double, and perturbative triple substitution results have been generated in part I4 of this collaborative study for 24 organic acceptor molecules that typify the search for photovoltaic materials. Because subtle nuclear relaxation effects usually are implicit in experiments on molecular IPs and EAs, such data are less suitable for judging the performance of a computational method. Part II5 considered the performance of long-range corrected hybrid functionals in generating predictions based on Kohn−Sham orbital energies or on non-self-consistent GW (see below) calculations. In Part III,6 many combinations of © 2016 American Chemical Society

Received: September 10, 2015 Published: January 5, 2016 627

DOI: 10.1021/acs.jctc.5b00872 J. Chem. Theory Comput. 2016, 12, 627−637

Article

Journal of Chemical Theory and Computation

Figure 1. Electron-accepting organic molecules.

Møller−Plesset fluctuation potential also are capable of describing correlation final states that bear little resemblance to those encompassed by KT. However, they also require greater computational effort. Because EP methods are based on the Dyson equation, they are related to the ∑ = iGW (or simply GW) family of approximations that have been employed extensively in solidstate physics.22,23 Whereas EP methods used in molecular electronic structure theory ultimately are implemented in a basis of atom-centered Gaussian functions, GW techniques have been expressed in terms of a variety of one-electron functions, including Gaussians and numerical atomic orbitals.24 In GW calculations, the intermediate step of orbital optimization, which is necessary to define a zeroth order EP or one-electron Green function, is usually a Kohn−Sham instead of a HF selfconsistent field. Several criteria of self-consistency may be used to classify the GW methods. Searches for poles of the EP (or G) may or may not be self-consistent with respect to the energy (or frequency) parameter. In addition, iterative GW approximations may seek self-consistency in G, the random-phase, screened-interaction matrix (W), or both. Recent attempts to optimize GW methodology for organic molecules have focused on propitious combinations of self-consistency criteria and density functionals used to define the zeroth order EP, G0, and screenedinteraction matrix, W0.25−32 To facilitate comparisons between related approaches that have emerged from solid-state and molecular traditions, tests of EP methods on a set organic electron-accepting molecules have been performed. The most important classes of compounds that are prominent in research on organic photovoltaic devices are included in the test set (Figure 1). Electron-withdrawing substituents such as halogens, cyano. and nitro groups are typical of these molecules, which often have lone pairs. In such cases, electron correlation may affect the order of cationic or anionic states with singly occupied π or lone-pair orbitals. The

test set therefore constitutes a challenge for methods that aim to describe such effects. To identify the most efficient EP methods, the arithmetic and memory demands of several approximations will be compared. This information suffices to identify the EP methods that are most useful in calculating the electron binding energies of molecular components of photovoltaic devices. The Conclusions include a discussion of additional methodological improvements pertaining to basis-set saturation that may enable greater computational efficiency to be realized.



ELECTRON-PROPAGATOR METHODS The p-q element of the EP matrix in its spectral form7,15 reads ⎧ ⟨Ψ0N |ap†|ΨrN − 1⟩⟨ΨrN − 1|aq|Ψ0N ⟩ ⎪ ∑ Gpq(E) = lim ⎨ η→ 0⎪ E − E0N + ErN − 1 − iη ⎩ r +

∑ r

⟨Ψ0N |aq|ΨrN + 1⟩⟨ΨrN + 1|ap†|Ψ0N ⟩ ⎫ ⎪ ⎬ N+1 N E − Er + E0 + iη ⎪ ⎭

(1)

where p and q are general spin-orbital indices; |ΨN0 ⟩ is the exact nondegenerate ground state of an N-electron system with energy EN0 ; ENr ± 1 and |ΨNr ± 1⟩ respectively denote energies and states with N ± 1 electrons; η is a positive infinitesimal constant that enables convergence of the Fourier transform from the time-dependent expression of the EP. In eq 1, the first term describes electron detachments from the ground state in which poles are negatives of IPs, whereas the second term describes electron attachments where poles are negatives of EAs. Poles of the usual, unperturbed EP, G0(E), where [G0(E)]pq = δpq(E − ϵp)−1

(2)

equal canonical HF orbital energies, ϵp, obtained by solving the HF equations, 628

DOI: 10.1021/acs.jctc.5b00872 J. Chem. Theory Comput. 2016, 12, 627−637

Article

Journal of Chemical Theory and Computation

FϕpHF = ϵpϕpHF

⎛ ⎞−1 d ∑(E) πr = ⎜1 − C†r (E) |E = Er Cr (E)⎟ ⎝ ⎠ dE

(3)

where F is the Fock operator with the familiar Coulomb exchange potential. In the frozen-orbital, single-determinant approximation employed in KT, IPs and EAs are negatives of occupied and virtual orbital energies. KT results usually are quantitatively inadequate and often produce erroneous assignments of ionic states. An improved EP can be obtained from G0(E) by employing the Dyson equation, which in its inverse matrix form reads G−1(E) = G−0 1(E) − ∑(E)

In the diagonal approximation, Dyson orbitals become proportional to normalized canonical HF orbitals such that ϕpDyson =

=N

(5)

Table 1. Scaling of Arithmetic and Storage Requirementsa

∫ dx2dx3 ... dxN ΨN(x1, x2 , x3 , ..., xN )

Ψ*N − 1, r(x 2 , x3 , ..., xN )

(6)

where xs is the space-spin coordinate for electron s. For EAs, where an electron is added to the reference state, the Dyson orbital is given by ϕrDyson(x1)

1/2

= (N + 1)

a

iterative

D2 P3 (P3+) IPs P3 (P3+) EAs D3 (OVGF) NR2 IPs NR2 EAs 2p-h TDA 3+ ADC(3)

OV2 O3V2 OV4 OV4 O2V3 OV4 OV4 OV4 O2V4

simple 2 3

OV O2V3 O2V3 O3V3 O2V4 O2V4 O2V4

storage OV2 OV3 V4 V4 OV3 V4 V4 V4 V4

O = number of occupied orbitals, V = number of virtual orbitals.

integrals) required by each EP method. Iterative, arithmetic bottlenecks arise from the search for poles. Simple, noniterative, arithmetic bottlenecks stem from couplings between simple (a) and triple (aaa†) ionization operators that need be evaluated only once. Storage requirements pertain to ERIs. Because the number of virtual orbitals, V, is generally much larger than the number of occupied orbitals, O, methods that employ ERIs with four virtual-orbital indices require more memory. Semidirect algorithms that regenerate these integrals also may be used.33,34 The D2 approximation requires a small set of four-index ERIs in the canonical HF basis where one index corresponds to a canonical orbital of interest and the other three indices correspond to 2h1p or 2p1h operators. The lengthiest calculation has an arithmetic scaling factor of OV2. ERI transformations to the HF basis constitute the bottleneck for D2 calculations. Outer valence Green function (OVGF) methods also assume a diagonal self-energy matrix and describe low-energy transitions.9,10,35 All versions of OVGF are based on the diagonal third-order (D3) terms in the self-energy and include estimates of higher-order contributions. Three frequently applied procedures have been denominated as versions A, B, and C.9,10,36 An algorithm that introduces several numerical criteria for the choosing among these alternatives has been developed9,10,36 and implemented.37 P3, the partial third-order quasiparticle method,38,39 omits several terms in the D3 self-energy. For IPs and EAs, this approximation retains all second- order terms and has,

(7)

The pole strength for a given Dyson orbital,

∫ |ϕrDyson(x)|2 dx

(8)

approaches unity when correlation effects are weak and approaches zero when such effects are strong. For IPs and EAs where KT provides a reasonable qualitative description of initial and final states, pole strengths are generally above 0.85. For correlation final states with shake-up (two-hole−oneparticle, or 2h1p) character for IPs or with core-excited (twoparticle−one-hole, or 2p1h) character for EAs, pole strengths below 0.8 are typical. Poles of the electron propagator correspond to eigenvalues obtained from [F + Σ(E)]C(E) = EC(E)

(9)

where the argument of the self-energy matrix is equal to an eigenvalue. The corresponding eigenvector, C(E), provides a linear combination of orbitals that is proportional to the Dyson orbital. The diagonal self-energy, or quasiparticle, approximation leads to a simplified form of the Dyson equation that reads E = ϵp + Σpp(E)

method

∫ dx2dx3 ... dxN +1ΨN +1,r

(x1 , x 2 , x3 , ..., xN + 1)Ψ*N (x 2 , x3 , ..., xN + 1)

πr =

(13)

Diagonal self-energy approximations have been applied widely.17−21 Employment of eq 10 with the second-order self-energy defines the diagonal second order, or D2, approximation. Table 1 lists the scaling of arithmetic operations and storage of transformed electron repulsion integrals (ERIs, or Coulomb

where F is now the Fock operator generated by the referencestate density matrix. Er is the eigenvalue of the self-consistent, energy-dependent operator F + ∑(E) and ϕDyson is a Dyson r orbital. For a vertical IP from an N-electron reference state that is associated with a final state r, the Dyson orbital reads 1/2

(12)

⎛ ⎞−1 dΣpp(E) ⎜ πp = ⎜1 − |E = Ep ⎟⎟ dE ⎝ ⎠

where ∑(E), the self-energy, describes final-state orbital relaxation and electron correlation. IPs and EAs may be obtained from the equivalent (quasiparticle) expression15

ϕrDyson(x1)

πp ϕpHF

where pole strengths are evaluated by the formula

(4)

[F + ∑(Er )]ϕrDyson = Er ϕrDyson

(11)

(10)

Correlation contributions to matrix elements of F in eq 5 are absorbed into ∑(E). Pole strengths are obtained from 629

DOI: 10.1021/acs.jctc.5b00872 J. Chem. Theory Comput. 2016, 12, 627−637

Article

Journal of Chemical Theory and Computation

Figure 2. Ionization potential error distribution (eV): Diagonal EP methods.

respectively, some 2h1p and 2p1h terms (e.g., ring and ladder diagrams) in third order. (Note that separate self-energy approximations are employed for IPs and EAs.) For IPs (EAs), the third-order 2ph (2hp) terms are omitted. Constant thirdorder terms are eliminated in both versions. The renormalized partial third-order, or P3+, method is obtained by introducing a renormalization factor in some terms.40

Nondiagonal self-energy methods also are examined presently. The two-particle−one-hole, Tamm−Dancoff approximation9 (2ph-TDA) includes all first-order couplings among h, p, 2h1p and 2p1h operators and retains ring, ladder, and mixed ring−ladder terms in all orders. In the 3+ method, higher-order couplings that yield all thirdorder terms in the self-energy are calculated.10 Some higher630

DOI: 10.1021/acs.jctc.5b00872 J. Chem. Theory Comput. 2016, 12, 627−637

Article

Journal of Chemical Theory and Computation

Figure 3. Ionization potential error distribution (eV): Nondiagonal EP methods.

TDA, 3+, ADC(3), and NR2 methods. Complete basis set (CBS) effects were estimated by extrapolating these results with an inverse cubic function.48 All pole strengths for the IPs and EAs under consideration exceeded 0.85 and confirm the qualitative validity of perturbative improvements to KT. A complete set of calculated results is provided in the Supporting Information.

order terms also are included in 3+. The third-order algebraic diagrammatic construction, or ADC(3), method surpasses 3+ by adding fourth-order and some higher-order terms in the energy-independent part of the self-energy matrix.9 The nondiagonal, renormalized, second-order (NR2) methods for IPs and EAs retain some third-order and higherorder terms (including those that describe final-state relaxation of orbitals) and are complete only in second order in the selfenergy.41 For IPs, NR2 is the least demanding of the nondiagonal methods examined below. For EAs, 3+ and NR2 have the smallest arithmetic and scaling requirements.



RESULTS AND DISCUSSION Histograms in Figures 2−5 show the errors in IPs and EAs made by each combination of self-energy approximation and basis set. In the top left corner of each histogram, the mean error (μ), standard deviation (σ), and mean unsigned error (MUE) also are displayed. Lack of convergence in HF calculations with the aug-cc-pVTZ basis on anthracene and phenazine and in constant self-energy algorithms in ADC(3)/ aug-cc-pVDZ calculations on acridine require that averages be based on less than 24 data in some cases. In all cases, IPs and



COMPUTATIONAL METHODS The present calculations were executed with a version of Gaussian 0937 that was modified to perform P3+ calculations or with the development version of Gaussian.42 Augmented, correlation-consistent double and triple ζ basis sets43−47 were used in combination with the D2, D3, OVGF, P3, P3+, 2ph631

DOI: 10.1021/acs.jctc.5b00872 J. Chem. Theory Comput. 2016, 12, 627−637

Article

Journal of Chemical Theory and Computation

Figure 4. Electron affinity error distribution (eV): Diagonal EP methods.

EAs increase as the basis set grows. Basis set extrapolations produce the largest predictions. Ionization Potentials. For IPs, the D2 method produces underestimates in most cases (Figure 2a−c). Even after basis set extrapolation, D2 has a MUE of 0.3 eV. The standard deviation of 0.26 eV indicates relatively high variability in the quality of the results. IPs for benzoquinone, maleic anhydride,

and naphthalenedione lie outside two standard deviations (2σ) from the mean error. Because D2 is the most efficient method, it can provide very approximate predictions when the alternatives are infeasible. The next most easily executed diagonal methods, P3 and P3+, have a tendency to overestimate IPs (Figure 2j−o), with the latter method showing a MUE of 0.2 eV and a small σ of 0.1 eV. (Note that P3+/aug-cc-pVTZ, 632

DOI: 10.1021/acs.jctc.5b00872 J. Chem. Theory Comput. 2016, 12, 627−637

Article

Journal of Chemical Theory and Computation

Figure 5. Electron affinity error distribution (eV): Nondiagonal EP methods.

As with P3 and P3+, there are no predicted IPs that lie outside 2σ. NR2 produces more overestimates than underestimates. This combination of precision, accuracy, and computational efficiency makes NR2 the preferred method for IP calculations. Larger measures of error and slower execution (see arithmetic scaling factors in Table 1) characterize the 2p-h TDA (Figure 3a−c). Nitrobenzonitrile is the outlier in Figure 3. This method is not competitive for IP calculations. The 3+ method’s MUE of 0.25 eV (Figure 3j−l) and greater computational requirements render it less promising than NR2. (The limitations of the perturbative order concept are evident again.49) Benzoquinone, maleic anhydride, and naphthalenedione have IPs that are outside 2σ (about 0.44 eV) of standard values. ADC(3), the only method with iterative sixth-power arithmetic scaling, produces error measures that are competitive with those of NR2 (Figure 3g−i). In fact, its mean error and MUE are the lowest of any method for IPs. However, NR2 has a smaller

P3+/aug-cc-pVDZ and P3/aug-cc-pVDZ results produce fortuitously good agreement with standard values.) None of the molecules have errors that are greater than 2σ for P3 and P3+. In general, the higher-order formula of P3+ mitigates the overestimates of P3 by about 0.1 eV. This combination of accuracy and efficiency makes P3+ an attractive method for IPs. D3 and OVGF, the most computationally demanding of the diagonal methods, produce MUEs near 0.2 eV (Figure 2d−i). Their standard deviations are larger than those of P3 and P3+, especially those of D3. D3 errors for benzoquinone, maleic anhydride, and tetrachloroisobenzofurandione are greater than 2σ. Only the OVGF result for maleic anhydride lies outside 2σ (about 0.34 eV in this case) of standard values. The oscillations of perturbative total energy methods are similar to the present D2 and D3 results and illustrate the limitations of the order concept in designing effective self-energy approximations.49 The fastest of the nondiagonal methods, NR2, yields a MUE of 0.14 eV with a standard deviation of 0.12 eV (Figure 3d−f). 633

DOI: 10.1021/acs.jctc.5b00872 J. Chem. Theory Comput. 2016, 12, 627−637

Article

Journal of Chemical Theory and Computation standard deviation. Only naphthalenedione has an IP that is off by more than 2σ. Electron Affinities. Although D2 results with the aug-ccpVDZ basis are in good agreement with standard values, the effects of basis improvement make this method of relatively little predictive value (Figure 4a−c). However, the standard deviation of 0.16 eV is lower for EAs than for IPs. Only the case of benzonitrile produces a result that is off by more than 2σ. Because of its modest computational requirements, D2 may be of use when other methods are infeasible. The best predictions with the diagonal self-energy approximation belong to the P3 method (Figure 4j−l). The average error, MUE, and standard deviation are all near 0.1 eV. All P3 results are within 2σ of reference values. Higher-order estimates provided by P3+ lead to unacceptable disagreements with standard values (Figure 4m−o). The benzonitrile result is off by more than 2σ. The D3 method has a marked tendency to underestimate EAs and has error indices that are unacceptably large (Figure 4d−f). These tendencies are ameliorated to some extent (Figure 4g−i) by the OVGF method. (The worst outlier for D3 and OVGF, with errors of more than 2σ, is tetrachloroisobenzofurandione.) However, the two latter methods appear to have no systematic advantage over P3. For the nondiagonal methods, the most accurate alternative, NR2, also is the most computationally efficient. The average error, standard deviation, and MUE are all near 0.1 eV (Figure 5d−f). Only TCNE’s result is off by barely more than 2σ (i.e., 0.20 eV). This method has a slight tendency to overestimate EAs. Overestimates produced by 2p-h TDA (Figure 5a−c) grow worse with basis set improvements and become unacceptably large once basis set extrapolations are made. (The 2σ outlier is nitrobenzonitrile.) This method is not competitive with P3 or NR2. The 3+ method (Figure 5j−l) tends to underestimate EAs. (Here, bodipy is the 2σ outlier.) Its large errors and disadvantageous computational characteristics also make it uncompetitive. ADC(3)’s predictive accuracy is close to that of NR2 (Figure 5g−i). Its standard deviation is somewhat lower than that of NR2. ADC(3) produces no 2σ outliers. ADC(3) has a tendency to underestimate EAs that contrasts with NR2’s overestimates. Higher Ionization Potentials. Calculations that produce erroneous orderings of cationic states may lead to fortuitous agreement or anomalous disagreement with experiment. Differences between ionic state orderings from KT and EP methods are commonplace in heterocyclic molecules. Therefore, several IPs and EAs were calculated for molecules with small energy differences between the highest occupied or lowest virtual canonical HF orbital energies. For naphthalenedione, three cationic states listed in Table 2 are close in energy. (EP results demonstrate that a fourth state predicted to be nearby at the KT level is less stable than the other three.) Methods are listed in the order of their computational difficulty. Only the 2p-h TDA, P3+, and NR2 methods agree with the benchmark assignment of the 2B2 label to the ground state of the cation. D2 produces fortuitously good agreement with the benchmark but predicts that another cationic state is 1 eV lower in energy. P3, D3, OVGF, 3+, and ADC(3) predict IPs that are close to the standard value but with the wrong state assignment. Close agreement with the benchmark is obtained for the P3+ and NR2 methods. Values from G0W0+SOSEX@PBE calculations are lower than NR2 results by 0.2−0.3 eV. The experimental value for the lowest IP is lower than the benchmark by almost 0.3 eV. The three

Table 2. Extrapolated Naphthalenedione IPs 2

method

2

B2 (eV)

KT D2 P3 P3+ D3 OVGF NR2 2p-h TDA 3+ ADC(3) benchmark G0W0+SOSEX@PBE experiment51

B1 (eV)

11.86 9.96 10.33 10.03 10.90 10.09 9.80 8.77 10.38 10.37 9.88 9.54 9.60

2

A2 (eV)

8.98 9.73 10.25 10.14 9.89 9.84 10.06 9.46 9.74 9.93

11.86 8.98 10.07 10.03 9.85 9.79 9.95 9.38 9.72 9.90

9.75 9.90

9.73 9.77

experimental results are exceeded by NR2 predictions by about 0.2 eV. Four cationic states are close in energy for benzoquinone. Table 3 shows that the orderings of final states from KT, D2, Table 3. Extrapolated Benzoquinone IPs method KT D2 P3 P3+ D3 OVGF NR2 2p-h TDA 3+ ADC(3) benchmark G0W0+SOSEX@PBE experiment52

2

B2u (eV) 12.70 9.42 10.89 10.56 11.77 10.90 10.35 9.88 11.02 11.02 10.30 10.41

2

B3g (eV) 12.06 9.46 10.65 10.40 11.13 10.43 10.13 9.17 10.58 10.67 10.27 9.95 10.11

2

B3u (eV)

2

B1g (eV)

11.38 10.92 11.42 11.34 11.24 11.00 11.17 10.66 11.07 11.12

11.23 11.26 11.46 11.39 11.07 11.07 11.34 10.86 10.98 11.20

10.76 11.06

10.89 11.23

and D3 are incorrect. NR2, P3+, and OVGF are in close agreement with the standard value, with NR2 giving the best result. G0W0+SOSEX@PBE results are lower than NR2 IPs by 0.1−0.4 eV. The experimental value for the lowest IP is lower than the standard calculation by 0.16 eV. Discrepancies between NR2 predictions and the four experimental results are 0.11 eV or less.



CONCLUSIONS Figure 6 summarizes the mean unsigned errors (MUEs) obtained with various ab initio electron-propagator (EP) methods for vertical ionization potentials (IPs) and electron affinities (EAs) of the test set of electron-accepting molecules. For the methods that make the diagonal self-energy approximation in the canonical Hartree−Fock (HF) basis, MUEs of 0.2 eV may be obtained for IPs calculated with P3+, D3, and OVGF. (The D2 method may provide semiquantitative predictions with much lower computational effort.) P3+ requires the least arithmetic of these fifth-power-scaling alternatives. OVGF, which is a computationally trivial extension of D3, provides a good compromise of accuracy and efficiency. Methods with nondiagonal self-energy approximations entail the introduction of sixth-power steps. Of these methods, the most efficient is NR2, which also produces a MUE of 0.15 eV. 634

DOI: 10.1021/acs.jctc.5b00872 J. Chem. Theory Comput. 2016, 12, 627−637

Article

Journal of Chemical Theory and Computation

Figure 6. Mean unsigned errors (eV).



The more arithmetically demanding ADC(3) method produces a slightly lower MUE for IPs. For EA calculations, the only diagonal self-energy approximation with a MUE below 0.2 eV is P3. Errors with the other diagonal methods are unacceptably large. The most accurate of the nondiagonal methods, NR2, with a MUE of 0.1 eV, also requires the least arithmetic and has a slight tendency to overestimate EAs. Predictions with ADC(3) have comparable accuracy and a contrasting trend of underestimating EAs. For IP calculations, P3+ provides a good combination of accuracy and efficiency. P3 is preferable for EA calculations. Although the A, B, and C versions of OVGF introduce no numerical parameters, the OVGF selection procedure does. When longer calculations are practical, the NR2 method offers another favorable combination of accuracy and efficiency. Its requirements for transformed integrals are no greater than those of P3 or P3+. NR2, with its treatment of couplings between 2hp operators that is complete through first order, is more able than P3+ to correctly describe higher IPs with lower pole strengths. P3+ is an expedient alternative for larger molecules. P3, P3+, and NR2 have the following feature in common, final-state orbital relaxation effects are described through third order. NR2 retains ring and ladder relaxation terms in all orders. Efficient estimates of NR2, P3, and P3+ calculations with large basis sets may be made by taking advantage of composite methods and improved virtual orbitals. The former methods postulate parallel behavior between D2 and higher-order predictions with respect to basis set size. By performing D2 calculations with large basis sets and higher-order calculations with relatively small basis sets, accurate estimates may be made efficiently. Quasiparticle virtual orbitals50 that are adapted to specific IPs may be used to increase the efficiency of calculations with large basis sets. Implementation and characterization of these techniques are currently in progress.

ASSOCIATED CONTENT

* Supporting Information S

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.jctc.5b00872. Data on EP calculations and accompanying statistical analysis. (PDF)



AUTHOR INFORMATION

Corresponding Author

*Phone: 334-844-6998. E-mail: [email protected]. Notes

The authors declare no competing financial interest.



ACKNOWLEDGMENTS Work at Tulane University was supported by the Louisiana Alliance for Simulation-Guided Materials Applications (LASiGMA), funded by National Science Foundation award EPS1003897. Computer time was provided by the National Energy Research Scientific Computing Center (NERSC), which is supported by the Office of Science of the U.S. Department of Energy under contract DE-AC02-05CH11231. We thank Professor Patrick Rinke of Aalto University of Science, Aalto, Finland, for useful discussions.



REFERENCES

(1) Grätzel, M. Photoelectrochemical Cells. Nature 2001, 414, 338− 344. (2) Brédas, J.-L. Organic Photovoltaics. Energy Environ. Sci. 2009, 2, 251−261. (3) Delgado, J. L.; Bouit, P.-A.; Filippone, S.; Herranz, M. A.; Martín, N. Organic Photovoltaics: A Chemical Approach. Chem. Commun. 2010, 46, 4853−4865. (4) Richard, R. M.; Marshall, M. S.; Dolgounitcheva, O.; Ortiz, J. V.; Brédas, J. L.; Marom, N.; Sherrill, C. D. Accurate Ionization Potentials and Electron Affinities of Acceptor Molecules. I. Reference Data at the CCSD(T) Complete Basis Set Limit. J. Chem. Theory Comput. 2016, DOI: 10.1021/acs.jctc.5b00875. 635

DOI: 10.1021/acs.jctc.5b00872 J. Chem. Theory Comput. 2016, 12, 627−637

Article

Journal of Chemical Theory and Computation (5) Gallandi, L.; Marom, N.; Rinke, P.; Körzdörfer, T. Accurate Ionization Potentials and Electron Affinities of Acceptor Molecules II: Non-empirically Tuned Long-range Corrected Hybrid Functionals. J. Chem. Theory Comput. 2016, DOI: 10.1021/acs.jctc.5b00873. (6) Knight, J.; Wang, X.; Gallandi, L.; Ren, X.; Rinke, P.; Körzdörfer, T.; Marom, N.; Dolgounitcheva, O.; Ortiz, J. V. Accurate Ionization Potentials and Electron Affinities of Acceptor Molecules for Organic Photovoltaics III: Benchmarking the Accuracy of GW Methods. J. Chem. Theory Comput. 2016, DOI: 10.1021/acs.jctc.5b00871. (7) Linderberg, J.; Ö hrn, Y. Propagators in Quantum Chemistry; Wiley-Interscience: Hoboken NJ, 2004. (8) Herman, M. F.; Freed, K. F.; Yeager, D. L. Analysis And Evaluation Of Ionization Potentials, Electron Affinities, And Excitation Energies By The Equations Of Motion-Green’s Function Method. Adv. Chem. Phys. 1981, 48, 1−69. (9) von Niessen, W.; Schirmer, J.; Cederbaum, L. S. Computational Methods For The One-Particle Green’s Function. Comput. Phys. Rep. 1984, 1, 57−125. (10) Ortiz, J. V. The Electron Propagator Picture Of Molecular Electronic Structure. In Computational Chemistry: Reviews of Current Trends; Leszcynski, J., Ed.; World Scientific: Singapore, 1997; Vol. 2, pp 1−61. (11) Ortiz, J. V. Toward an Exact One-Electron Picture of Chemical Bonding. In A Tribute to the Worksof Yngve Ö hrn; Löhrn, A., Ed.; Adv. Quantum Chem.; Academic Press, 1999; Vol. 35, pp 33−52. (12) Simons, J. Equation of Motion Methods for Computing Electron Affinities and Ionization Potentials. In Theory and Applications of Computational Chemistry: The First Forty Years; Dykstra, C. E., Ed.; Elsevier: Amsterdam, 2005; pp 443−464. (13) Flores-Moreno, R.; Melin, J.; Dolgounitcheva, O.; Zakrzewski, V. G.; Ortiz, J. V. Three Approximations To The Nonlocal And Energy-Dependent Correlation Potential In Electron Propagator Theory. Int. J. Quantum Chem. 2010, 110, 706−715. (14) Flores-Moreno, R.; Ortiz, J. V. Efficient and Accurate Electron Propagator Methods and Algorithms. In Practical Aspects of Computational Chemistry; Leszczynski, J., Shukla, M. K., Eds.; Springer, 2010; pp 1−17. (15) Ortiz, J. V. Electron Propagator Theory: An Approach To Prediction And Interpretation In Quantum Chemistry. Wiley Interdisciplinary Reviews: Computational Molecular Science 2013, 3, 123−142. (16) Danovich, D. Green’s Function Methods For Calculating Ionization Potentials, Electron Affinities, And Excitation Energies. Wiley Interdisciplinary Reviews: Computational Molecular Science 2011, 1, 377−387. (17) Ortiz, J. V.; Zakrzewski, V. G.; Dolgounitcheva, O. One-Electron Pictures of Electronic Structure: Propagator Calculations on Photoelectron Spectra of Aromatic Molecules. In Conceptual Perspectives in Quantum Chemistry; Calais, J.-L., Kryachko, E., Eds.; Springer: The Netherlands, 1997; Vol. 3, pp 465−517. (18) Dolgounitcheva, O.; Zakrzewski, V.; Ortiz, Z. V. Electron Propagator Calculations on the Ionization Energies of Nucleic Acid Bases, Base-Water Complexes and Base Dimers. In Fundamental World of Quantum Chemistry: A Tribute to the Memory of Per-Olov Löwdin; Brändas, E. J., Kryachko, E. S., Eds.; Kluwer: Dordrecht, 2003; Vol. 2, pp 525−55. (19) Zakrzewski, V. G.; Dolgounitcheva, O.; Zakjevskii, A. V.; Ortiz, J. V. Ab Initio Electron Propagator Methods: Applications To Nucleic Acids Fragments And Metallophthalocyanines. Int. J. Quantum Chem. 2010, 110, 2918−2930. (20) Zakrzewski, V. G.; Dolgounitcheva, O.; Zakjevskii, A. V.; Ortiz, J. Ab Initio Electron Propagator Methods: Applications to Fullerenes and Nucleic Acid Fragments. In Annu. Rep. Comput. Chem.; Wheeler, R. A., Ed.; Elsevier, 2010; Vol. 6, Chapter 6, pp 79−94. (21) Zakrzewski, V. G.; Dolgounitcheva, O.; Zakjevskii, A. V.; Ortiz, J. V. Ab Initio Electron Propagator Calculations On Electron Detachment Energies Of Fullerenes, Macrocyclic Molecules, And Nucleotide Fragments. Adv. Quantum Chem. 2011, 62, 105−136.

(22) Hedin, L. New Method For Calculating The One-Particle Green’s Function With Applications To The Electron-Gas Problem. Phys. Rev. 1965, 139, A796−A823. (23) Hybertsen, M. S.; Louie, S. G. Electron Correlation In Semiconductors And Insulators: Band Gaps And Quasiparticle Energies. Phys. Rev. B: Condens. Matter Mater. Phys. 1986, 34, 5390− 5413. (24) Ren, X.; Sanfilippo, A.; Rinke, P.; Wieferink, J.; Tkatchenko, A.; Reuter, K.; Blum, V.; Scheffler, M. An Accurate Resolution Of The Identity Approach To Hartree-Fock, Hybrid Functionals, MP2, RPA And GW With Numeric Atom-Centered Basis Functions. New J. Phys. 2012, 14, 053020. (25) Rinke, P.; Qteish, A.; Neugebauer, J.; Freysoldt, C.; Scheffler, M. Combining GW Calculations With Exact-Exchange Density-Functional Theory: An Analysis Of Valence-Band Photoemission For Compound Semiconductors. New J. Phys. 2005, 7, 126. (26) Fuchs, F.; Furthmüller, J.; Bechstedt, F.; Shishkin, M.; Kresse, G. Quasiparticle Band Structure Based On A Generalized Kohn-Sham Scheme. Phys. Rev. B: Condens. Matter Mater. Phys. 2007, 76, 11509. (27) Marom, N.; Moussa, J. E.; Ren, X.; Tkatchenko, A.; Chelikowsky, J. R. Electronic Structure Of Dye-Sensitized TiO2 Clusters From Many-Body Perturbation Theory. Phys. Rev. B: Condens. Matter Mater. Phys. 2011, 84, 245115. (28) Faber, C.; Attaccalite, C.; Olevano, V.; Runge, E.; Blase, X. FirstPrinciples GW Calculations For DNA And RNA Nucleobases. Phys. Rev. B: Condens. Matter Mater. Phys. 2011, 83, 115123. (29) Marom, N.; Caruso, F.; Ren, X.; Hofmann, O. T.; Körzdörfer, T.; Chelikowsky, J. R.; Rubio, A.; Scheffler, M.; Rinke, P. Benchmark Of GW methods For Azabenzenes. Phys. Rev. B: Condens. Matter Mater. Phys. 2012, 86, 247127. (30) Caruso, F.; Rinke, P.; Ren, X.; Scheffler, M.; Rubio, A. Unified Description Of Ground And Excited States Of Finite Systems: The Self-Consistent GW Approach. Phys. Rev. B: Condens. Matter Mater. Phys. 2012, 86, 081102. (31) Körzdörfer, T.; Marom, N. Strategy For Finding A Reliable Starting Point For G0W0 Demonstrated For Molecules. Phys. Rev. B: Condens. Matter Mater. Phys. 2012, 86, 041110. (32) Atalla, V.; Yoon, M.; Caruso, F.; Rinke, P.; Scheffler, M. Hybrid Density Functional Theory Meets Quasiparticle Calculations: A Consistent Electronic Structure Approach. Phys. Rev. B: Condens. Matter Mater. Phys. 2013, 88, 165122. (33) Zakrzewski, V. G.; Ortiz, J. V. Semidirect Algorithms In Electron Propagator Calculations. Int. J. Quantum Chem. 1994, 52, 23−27. (34) Zakrzewski, V. G.; Ortiz, J. V. Semidirect Algorithms For ThirdOrder Electron Propagator Calculations. Int. J. Quantum Chem. 1995, 53, 583−590. (35) Cederbaum, L. S. One-Body Green’s Function For Atoms And Molecules: Theory And Application. J. Phys. B: At. Mol. Phys. 1975, 8, 290−303. (36) Zakrzewski, V. G.; Ortiz, J. V.; Nichols, J. A.; Heryadi, D.; Yeager, D. L.; Golab, J. T. Comparison Of Perturbative And Multiconfigurational Electron Propagator Methods. Int. J. Quantum Chem. 1996, 60, 29−36. (37) Frisch, M. J.; Trucks, G. W.; Schlegel, H. B.; Scuseria, G. E.; Robb, M. A.; Cheeseman, J. R.; Scalmani, G.; Barone, V.; Mennucci, B.; Petersson, G. A.; Nakatsuji, H.; Caricato, M.; Li, X.; Hratchian, H. P.; Izmaylov, A. F.; Bloino, J.; Zheng, G.; Sonnenberg, J. L.; Hada, M.; Ehara, M.; Toyota, K.; Fukuda, R.; Hasegawa, J.; Ishida, M.; Nakajima, T.; Honda, Y.; Kitao, O.; Nakai, H.; T. Vreven, J.; Montgomery, J. A.; Peralta, J. E.; Ogliaro, F.; Bearpark, M.; Heyd, J. J.; Brothers, E.; Kudin, K. N.; Staroverov, V. N.; Kobayashi, R.; Normand, J.; Raghavachari, K.; Rendell, A.; Burant, J. C.; Iyengar, S. S.; Tomasi, J.; Cossi, M.; Rega, N.; Millam, J. M.; Klene, M.; Knox, J. E.; Cross, J. B.; Bakken, V.; Adamo, C.; Jaramillo, J.; Gomperts, R.; Stratmann, R. E.; Yazyev, O.; Austin, A. J.; Cammi, R.; Pomelli, C.; Ochterski, J. W.; Martin, R. L.; Morokuma, K.; Zakrzewski, V. G.; Voth, G. A.; Salvador, P.; Dannenberg, J. J.; Dapprich, S.; Daniels, A. D.; Farkas, O.; Foresman, J. B.; Ortiz, J. V.; Cioslowski, J.; Fox, D. J. Gaussian 09, Revision D.01; Gaussian, Inc.: Wallingford, CT, 2009. 636

DOI: 10.1021/acs.jctc.5b00872 J. Chem. Theory Comput. 2016, 12, 627−637

Article

Journal of Chemical Theory and Computation (38) Ortiz, J. V. Partial Third-Order Quasiparticle Theory: Comparisons For Closed-shell Ionization Energies And An Application To The Borazine Photoelectron Spectrum. J. Chem. Phys. 1996, 104, 7599−7605. (39) Ferreira, A. M.; Seabra, G.; Dolgounitcheva, O.; Zakrzewski, V. G.; Ortiz, J. V. Application and Testing of Diagonal, Partial ThirdOrder Electron Propagator Approximations. In Quantum-Mechanical Prediction of Thermochemical Data; Cioslowski, J., Ed.; Understanding Chemical Reactivity; Springer: The Netherlands, 2001; Vol. 22, pp 131−160. (40) Ortiz, J. V. An Efficient, Renormalized Self-energy For Calculating The Electron Binding Energies Of Closed-shell Molecules And Anions. Int. J. Quantum Chem. 2005, 105, 803−808. (41) Ortiz, J. V. A Nondiagonal, Renormalized Extension Of Partial Third-order Quasiparticle Theory: Comparisons For Closed-shell Ionization Energies. J. Chem. Phys. 1998, 108, 1008−1014. (42) Frisch, M. J.; Trucks, G. W.; Schlegel, H. B.; Scuseria, G. E.; Robb, M. A.; Cheeseman, J. R.; Scalmani, G.; Barone, V.; Mennucci, B.; Petersson, G. A.; Nakatsuji, H.; Caricato, M.; Li, X.; Hratchian, H. P.; Bloino, J.; Janesko, B. G.; Izmaylov, A. F.; Lipparini, F.; Zheng, G.; Sonnenberg, J. L.; Liang, W.; Hada, M.; Ehara, M.; Toyota, K.; Fukuda, R.; Hasegawa, J.; Ishida, M.; Nakajima, T.; Honda, Y.; Kitao, O.; Nakai, H.; Vreven, T.; Throssell, K.; Montgomery, J. A.Jr.; Peralta, J. E.; Ogliaro, F.; Bearpark, M.; Heyd, J. J.; Brothers, E.; Kudin, K. N.; Staroverov, V. N.; Keith, T.; Kobayashi, R.; Normand, J.; Raghavachari, K.; Rendell, A.; Burant, J. C.; Iyengar, S. S.; Tomasi, J.; Cossi, M.; Rega, N.; Millam, J. M.; Klene, M.; Knox, J. E.; Cross, J. B.; Bakken, V.; Adamo, C.; Jaramillo, J.; Gomperts, R.; Stratmann, R. E.; Yazyev, O.; Austin, A. J.; Cammi, R.; Pomelli, C.; Ochterski, J. W.; Martin, R. L.; Morokuma, K.; Zakrzewski, V. G.; Voth, G. A.; Salvador, P.; Dannenberg, J. J.; Dapprich, S.; Parandekar, P. V.; Mayhall, N. J.; Daniels, A. D.; Farkas, O.; Foresman, J. B.; Ortiz, J. V.; Cioslowski, J.; Fox, D. J. Gaussian Development Version Revision I.04; Gaussian, Inc.: Wallingford, CT, 2009. (43) Dunning, T. H. Gaussian Basis Sets For Use In Correlated Molecular Calculations. I. The Atoms Boron Through Neon And Hydrogen. J. Chem. Phys. 1989, 90, 1007−1023. (44) Kendall, R. A.; Dunning, T. H.; Harrison, R. J. Electron Affinities Of The First-Row Atoms Revisited. Systematic Basis Sets And Wave Functions. J. Chem. Phys. 1992, 96, 6796−6806. (45) Woon, D. E.; Dunning, T. H. Gaussian Basis Sets For Use In Correlated Molecular Calculations. III. The Atoms Aluminum Through Argon. J. Chem. Phys. 1993, 98, 1358−1371. (46) Davidson, E. R. Comment On Dunning’s CorrelationConsistent Basis Sets. Chem. Phys. Lett. 1996, 260, 514−518. (47) Dunning, T. H.; Peterson, K. A.; Wilson, A. K. Gaussian Basis Sets For Use In Correlated Molecular Calculations. X. The Atoms Aluminum Through Argon Revisited. J. Chem. Phys. 2001, 114, 9244− 9253. (48) Helgaker, T.; Klopper, W.; Koch, H.; Noga, J. Basis-Set Convergence Of Correlated Calculations On Water. J. Chem. Phys. 1997, 106, 9639−9646. (49) Hirata, S.; Hermes, M. R.; Simons, J.; Ortiz, J. V. General-order many-body Green’s function method. J. Chem. Theory Comput. 2015, 11, 1595−1606. (50) Flores-Moreno, R.; Zakrzewski, V. G.; Ortiz, J. V. Quasi-Particle Virtual Orbitals In Electron Propagator Calculations. J. Chem. Phys. 2008, 128, 164105. (51) Millefiori, S.; Gulino, A.; Casarin, M. UV photoelectron spectra, reducion potentials and MO calculations of intramolecularly hydrogen-bonded naphtoquinones. J. Chim. Phys. 1990, 87, 317−330. (52) Brundle, C. R.; Robin, M. B.; Kuebler, N. A. Perfluoro effect in photoelectron spectroscopy. II. Aromatic molecules. J. Am. Chem. Soc. 1972, 94, 1466−1475.

637

DOI: 10.1021/acs.jctc.5b00872 J. Chem. Theory Comput. 2016, 12, 627−637