Size-dependence of non-empirically tuned DFT starting points for $G_0W_0$ applied to ${\pi}$-conjugated molecular chains

TL;DR

This paper investigates how the size of π-conjugated molecular chains influences the effectiveness of non-empirically tuned DFT starting points in $G_0W_0$ calculations for predicting ionization potentials, finding that IP-tuned hybrids perform well across sizes.

Contribution

It demonstrates that size-dependent non-empirical tuning of DFT functionals leads to accurate $G_0W_0$ IP predictions for long molecular chains.

Findings

Size dependence affects the tuning of hybrid functionals.

IP-tuned hybrids yield accurate $G_0W_0$ results across chain lengths.

Comparison with coupled-cluster calculations confirms accuracy.

Abstract

$G_0W_0$ calculations for predicting vertical ionization potentials (IPs) and electron affinities of molecules and clusters are known to show a significant dependence on the density functional theory (DFT) starting point. A number of non-empirical procedures to find an optimal starting point have been proposed, typically based on tuning the amount of HF exchange in the underlying hybrid functional specifically for the system at hand. For the case of $\pi$-conjugated molecular chains, these approaches lead to a significantly different amount of HF exchange for different oligomer sizes. In this study, we analyze if and how strongly this size dependence affects the ability of non-empirical tuning approaches to predict accurate IPs for $\pi$-conjugated molecular chains of increasing chain length. To this end, we employ three different non-empirical tuning procedures for the $G_0W_0$ starting point to calculate the IP of polyene oligomers up to 22 repeat units and compare the results to highly accurate coupled-cluster calculations. We find that, despite its size dependence, using an IP-tuned hybrid functional as a starting point for $G_0W_0$ yields excellent agreement with the reference data for all chain lengths.

Equations

$$E_i^{\mathrm{QP}}={ϵ}_i^{\mathrm{DFT}}+⟨{\varphi }_i^{\mathrm{DFT}}|\hat{\sum }(E_i^{\mathrm{QP}})-{\hat{v}}_{xc}|{\varphi }_i^{\mathrm{DFT}}⟩$$

$$\begin{array}{rl}\mathrm{DSLE}& =\mathrm{EA}(M-1)-\mathrm{IP}(M)\\ & ={ϵ}_{\mathrm{LUMO},\mathrm{cat}}-{ϵ}_{\mathrm{HOMO},\mathrm{neut}}.\end{array}$$

$$\frac{\delta E}{\delta f_{\mathrm{HOMO}}}={ϵ}_{\mathrm{HOMO}}.$$

$$\mathrm{IP}=E_M-E_{M-1}=-{ϵ}_{\mathrm{HOMO}}.$$

$$E_{xc}^{\mathrm{PBEh},\alpha }=\alpha E_x^{\mathrm{HF}}+(1-\alpha )E_x^{\mathrm{PBE}}+E_c^{\mathrm{PBE}}.$$

$${\Delta }_{\mathrm{IP}}=|-{ϵ}_{\mathrm{HOMO},\mathrm{neut}}^{\alpha }-(E_{\mathrm{neut}}^{\alpha }-E_{\mathrm{cat}}^{\alpha })|$$

$$\begin{array}{rl}{\Delta }_{\mathrm{DSLE}}& =\mathrm{EA}(M-1)-\mathrm{IP}(M)\\ & ={ϵ}_{\mathrm{LUMO},\mathrm{cat}}^{\alpha ,\mathrm{QP}}-{ϵ}_{\mathrm{HOMO},\mathrm{neut}}^{\alpha ,\mathrm{QP}},\end{array}$$

$${\Delta }_{G_0{W}_0}=|{ϵ}_{\mathrm{HOMO}}^{\alpha ,\mathrm{DFT}}-{ϵ}_{\mathrm{HOMO}}^{\alpha ,\mathrm{QP}}|$$

Full text

1 Introduction

An accurate prediction of the ionization potential (IP) of $\pi$-conjugated polymers is of utmost importance when considering them for applications in organic electronics. Due to their huge system size, accurate theoretical methods can not be used to investigate $\pi$-conjugated polymers directly. To circumvent this problem, relevant properties are often calculated for oligomers of increasing length and extrapolated to the infinite chain limit, therefore allowing the use of electronic structure methods such as density functional theory (DFT). This is known as the oligomer approach, which is widely followed for the prediction of polymer properties.1, 2, 3

Many body perturbation theory within the $GW$ approximation4, 5, 6 is a method increasingly applied to calculate charged excitation energies, that is, IPs and electron affinities (EAs) of molecules.7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17 It can be applied using either full7, 18, 19, 20, 17 or partial5, 8, 15, 21, 22, 17 self consistency, or as a single step correction on top of a DFT or Hartree-Fock (HF) calculation, hereafter referred to as $G_{0}W_{0}$. While fully self-consistent $GW$ calculations ($scGW$) are still scarce, non-self-consistent $G_{0}W_{0}$ is by far the most commonly used variant, as it has been shown to lead to accurate IPs at a fraction of the cost of partially or fully self-consistent schemes.18, 9, 23 However, while the $G_{0}W_{0}$ correction can be readily applied to any DFT calculation, a strong starting point dependence has been found,24, 15, 20, 25, 9, 26, 8, 13, 10, 11, 12, 27, 28, 29, 16, 30, 17 hence making the choice of the exchange-correlation functional in the underlying DFT calculation crucial for the accuracy of the calculated IPs and EAs.

The starting point dependence of $G_{0}W_{0}$ has been connected to mainly two effects. First, $G_{0}W_{0}$ is based on a perturbative correction, so it can only be expected to yield reliable results if the eigenvalues of the DFT starting point are already good approximations to their corresponding quasiparticle energies. Second, the validity of the quasiparticle approach hinges on a correct description of the screened interaction $W$, which is mainly determined by the systems polarizability and, hence, its band gap. Consequently, it is essential that the eigenvalue gap between the highest occupied molecular orbital (HOMO) and the lowest unoccupied molecular orbital (LUMO) of the DFT starting point approximates the true quasiparticle gap, i.e., the difference between the IP and EA.16

Semilocal exchange-correlation functionals fulfill neither of these two requirements, which is why using hybrid functional starting points that include a fraction of HF exchange has been found to be essential for obtaining accurate IPs and EAs from $G_{0}W_{0}$ calculations.25, 13, 29, 9, 10, 11, 12, 27, 28 In recent years, a number of different non-empirical tuning schemes have been proposed, aiming to find an optimal amount of HF exchange for the $G_{0}W_{0}$ starting point.

One of the most promising non-empirical tuning approaches is the so-called IP-tuning, 31, 32, 33, 34 which will be explained in more detail below. Due to their high accuracy for the prediction of IPs, EAs, and fundamental gaps, IP-tuned long-range corrected (LRC) hybrid functionals had long been suggested as a potentially interesting starting point for $G_{0}W_{0}$.35 Recently, it was demonstrated that, indeed, the combination of these functionals with a single step $G_{0}W_{0}$ correction produces excellent results for the IPs and EAs of a organic molecules, clearly outperforming typically used semilocal or global hybrid starting points.27, 36, 29, 28 In addition, IP-tuned hybrid functionals have been shown to be appropriate starting points for excited state $GW$ Bethe-Salpeter equation calculations.37

However, the combination of the IP-tuning procedure with the oligomer approach has been shown to lead to several problematic shortcomings,38, 39, 40 which can be traced back to the fact that IP-tuned functionals are size dependent.41, 42, 43, 44 Especially for highly conjugated molecular chains, such as the ones frequently used in organic electronics and photovoltaics, the amount of exact exchange obtained by the IP-tuning procedure (accounted for by the range-separation parameter $\omega$ in LRC functionals or the amount of HF exchange $\alpha$ in global hybrid functionals) decreases rapidly as the oligomer chain length increases, before gradually leveling off.41, 42, 40, 38, 43, 45, 44 For highly conjugated systems, convergence can only be found for oligomer chains of several nanometers.43 This strong size-dependence has been demonstrated to create unphysical trends for properties that depend heavily on the amount of exact exchange, such as bond-length alternation38 and optical band-gaps.39 The aim of this work is to investigate how this size dependence affects the reliability of non-empirically tuned hybrid functionals as a starting point for $G_{0}W_{0}$ when calculating the IPs of $\pi$-conjugated polymers of increasing chain length.

Polyacetylene (H$-$(HC$=$CH)${}_{N}-$H, see Fig. 1) serves as an ideal model for $\pi$-conjugated polymers. Its simple structure and small repeat units allow calculations for oligomers with many repeat units using highly accurate wavefunction methods such as coupled cluster singles, doubles, perturbative triples (CCSD(T)), thus, providing a reliable theoretical reference for computationally more efficient methods such as DFT and/or $G_{0}W_{0}$. At the same time, it shows the $\pi$-conjugated backbone found in all polymers of this class, exhibiting the typical pattern of single and double bonds. 46, 47, 48 The difference in bond length between single and double bonds, known as bond length alternation (BLA), is a defining property of $\pi$-conjugated polymers. Predicting the BLA and its evolution with chain-length is a notorious problem many standard computational methods struggle with. 49, 50, 51, 38 This problem is particularly severe since the BLA decisively influences the electronic and optical properties of organic $\pi$-conjugated molecular chains.

In the following methodology section we will further elaborate on the aforementioned difficulties of common density functionals to act as a reliable starting point for $G_{0}W_{0}$ and provide the theoretical background for the three non-empirical tuning procedures studied in this work. In the results section, we will then compare the $\Delta$SCF and HOMO IPs obtained from each tuning method to CCSD(T) reference calculations for the oligomers of polyacetylene up to a chain length of N=22. Finally a critical evaluation of the tuning methods’ applicability as a reliable starting point for $G_{0}W_{0}$ as well as their respective behavior with increasing oligomer size will follow.

2 Methodology

The aim of this section is to introduce the methods used in this article, especially highlighting the differences in the non-empirical tuning procedures.

The $\boldsymbol{G_{0}W_{0}}$ method

$GW$ is an approximation made in many-body perturbation theory, expressing the self-energy as a first order expansion in terms of the single particle Greens function, $G$, and the screened Coulomb interaction, $W$.4, 5, 6 Fully self-consistent $GW$ is computationally unviable for medium/large systems and results in comparatively large errors for IPs of organic molecules.29, 9, 12 In contrast, carrying out a single step perturbative correction on top of a DFT calculation, as done in $G_{0}W_{0}$ calculations, is significantly cheaper and has been shown to give excellent results for IPs when combined with a well-considered starting point.18, 9, 36, 29, 10, 12, 27, 28, 13, 52, 53, 23 The $G_{0}W_{0}$ quasiparticle energies are given by

[TABLE]

where ${\epsilon}_{i}^{\mathrm{DFT}}$ are the underlaying DFT eigenvalues, $\phi_{i}^{\mathrm{DFT}}$ are the corresponding orbitals, $\hat{\sum}(E_{i}^{\mathrm{QP}})$ is the $GW$ self-energy calculated using the DFT orbitals and eigenvalues, and $\hat{v}_{xc}$ is the exchange-correlation operator.54 As the eigenvalues and orbitals of the DFT calculation to which the $G_{0}W_{0}$ correction is applied are directly used as input for the correction, it is vital that these are good approximations to actual quasiparticle energies.

Deviation from straight line

The exact total energy of any many-electron system exhibits straight line behavior when going from M electrons to M$\pm$1 electrons via fractional electron numbers ($\delta$), while it shows a derivative discontinuity at integer electron numbers.55 The difference in slope when going from M-1 to M and from M to M+1 electrons equals the difference between the systems’ IP and EA. Standard DFT functionals, however, have been shown to not reproduce this behavior. 56, 57, 58 Fig. 2 shows the deviation from straight line behavior for fractional electron numbers for Octatetraene (N=4) obtained from a semilocal functional and the HF approach. Semilocal functionals, represented here by the functional of Perdew-Burke-Ernzerhof (PBE), show no derivative discontinuity and exhibit a convex slope between occupation numbers, thereby artificially lowering the energy for fractional electron numbers. This is known as the delocalization error of semilocal DFT. In contrast, the HF approach shows an exaggerated derivative discontinuity and a concave slope between the M and M$\pm$1 electron systems. This strongly stabilizes the energy of integer electron occupations, causing what is known as the localization error.

The integral of the curves in Fig. 2 is often referred to as the many-electron self interaction error (MESIE).57, 56

As calculations with fractional electron numbers are both time consuming and, more importantly, not possible for $GW$, it can be advantageous to follow an alternative method to quantify the curvature of the energy.59, 60, 52 This deviation from straight line error (DSLE) is defined as the energy difference between the EA of the cationic system, $\mathrm{EA_{c}}(M-1)$, and the IP of the neutral system, $\mathrm{IP}(M)$.

[TABLE]

Within the framework of DFT, $\mathrm{EA_{c}}(M-1)$ and $\mathrm{IP}(M)$ are represented by the LUMO eigenvalue of the cation, $\epsilon_{\mathrm{LUMO,cat}}$, and HOMO eigenvalue of the neutral molecule, $\epsilon_{\mathrm{HOMO,neut}}$, respectively. This definition of the DSLE relies on the fact that, according to Janak’s theorem,61 the energy slope of the fractionally occupied HOMO is given by the HOMO eigenvalue, i.e.,

[TABLE]

Consequently, the initial and final slopes of the total energy curve are given by the cations’ LUMO ($\delta=-1$) and the neutral molecules’ HOMO ($\delta=0$) eigenvalue, respectively. In Fig. 2, these eigenvalues are represented by short lines with the corresponding slopes, showing excellent agreement with the calculated full total energy difference curves. A table showing the general agreement between the DSLE found for initial/final slopes of fractional electron numbers and the cation LUMO and the neutral HOMO eigenvalues can be found in the supporting information. We therefore conclude that calculating the DSLE based on the neutral molecules’ HOMO and the cations’ LUMO eigenvalues provides for a reliable measure of the deviation from linearity without having to go through the numerically expensive fractional electron number calculations.

IP tuning

IP-tuning is a well established non-empirical tuning procedure based on enforcing the functional to obey the IP-theorem,31, 32, 33, 34 i.e.,

[TABLE]

This fixes the slope of the fractional electron number curve in Fig. 2 to be zero at $\delta=0$. Consequently, IP-tuning ensures not only the IP-Theorem but also (approximately) the straight-line behavior of the total energy curve when going from the neutral molecule to the cation. Global hybrid functionals have been previously discussed in the literature in the context of IP-theorem based tuning procedures.62, 63 Here, we have used the global hybrid functional PBEh for all tuning, which combines PBE with a fraction of exact exchange:

[TABLE]

The IP-tuned fraction of exact exchange, $\alpha$, is found by minimizing

[TABLE]

with $E^{\alpha}_{\mathrm{neut/cat}}$ representing the energy of the neutral and cationic systems calculated on the neutral geometry and $\epsilon^{\alpha}_{\mathrm{HOMO,neut}}$ is the HOMO eigenvalue of the neutral molecule.

IP-tuning is perhaps the most well known non-empirical tuning procedure since it has been extensively discussed in the literature and used for a wide range of applications.64, 65, 66, 67, 68, 40, 45, 33, 69, 70, 71, 72, 73, 74, 75, 63 Importantly, IP-tuning has originally been introduced and is, indeed, most frequently used as a computationally efficient alternative to $GW$ calculations, since it has been demonstrated to provide for very accurate IPs and EAs. In a recent theoretical benchmark study on a test set of organic acceptor molecules76, 36, 29 as well as in a study that directly compared to highly resolved gas-phase photoelectron spectroscopy measurements,27 IP-tuned LRC hybrid functionals outperformed $scGW$ and $G_{0}W_{0}$ calculations based on commonly used semilocal DFT or HF starting points. However, in the same studies it was demonstrated that the results of the IP-tuned functionals could be further improved by applying a $G_{0}W_{0}$ correction to the DFT eigenvalues. Here, we will study the performance of IP-tuned global hybrid functionals both with and without a $G_{0}W_{0}$ correction, where a special emphasis lies on the size dependence of these functionals and the potential consequences this may have for the application to long, $\pi$-conjugated molecular chains.

The system size dependence of the IP-tuned exact exchange parameter can be clearly observed when IP-tuning PBEh for the oligomers of polyacetylene. As demonstrated in Fig. 3, the IP-tuned fraction of exact exchange drops off as the chain length increases; starting at 0.70 for N=4 it drops to 0.51 for N=22. While a slight sloping off can be observed towards longer chain lengths, no convergence can be observed for the oligomer sizes studied here. Since a different amount of HF exchange and, hence, a different exchange-correlation functional is used for different oligomer lengths, the IP-tuning procedure shows strong size dependence and is therefore likely to fail in predicting accurate values for the polymer limit via extrapolation of the results obtained from the individual oligomers.

DSLE tuning

Along with the fractional particle curves and DFT eigenvalue slopes, Fig. 2 also shows the slopes derived from the corresponding quasiparticle energies $\epsilon^{\mathrm{QP}}_{\mathrm{LUMO,cat}}$ and $\epsilon^{\mathrm{QP}}_{\mathrm{HOMO,neut}}$ obtained from a $G_{0}W_{0}$ correction to the PBE and HF eigenvalues, respectively. This highlights that, depending on the starting point, a significant DSLE may still remain after the $G_{0}W_{0}$ correction. In other words, $G_{0}W_{0}$ based on a PBE or HF starting point is not free from DSLE.

This finding has recently lead Rinke et al. to propose an alternative non-empirical tuning procedure for the DFT starting point of $G_{0}W_{0}$.52 While the IP-tuning procedure significantly reduces the DSLE in the DFT starting point, the central idea of the DSLE tuning is to enforce the straight line behavior after the $G_{0}W_{0}$ correction. This is achieved by tuning the amount of HF exchange in the PBEh starting point in order to minimize

[TABLE]

where EA${}_{c}(M-1)$ is the EA of the cationic system and IP$(M)$ is the IP of the neutral system. It is important to remember that, when following the DSLE-based tuning to minimize Eq. (7), the eigenvalues considered are not the DFT eigenvalues but the quasiparticle energies obtained from the $G_{0}W_{0}$ correction.

When tuning $G_{0}W_{0}$ on PBEh for a series of polyacetylene oligomers of increasing length, a decrease in $\alpha$ with increasing chainlength similar to IP-tuning can be observed for the DSLE tuned functional (see Fig. 3). The fractions of exact exchange are, however, much lower for all chain length, with an $\alpha$ of 0.38 for N=4 which falls off to 0.20 for N=22.

$\boldsymbol{G_{0}W_{0}}$ tuning

While both tuning procedures introduced above are constructed around the idea of enforcing straight line behavior, Scheffler et al. came up with an entirely different approach to non-empirically tune the fraction of exact exchange in a hybrid functional.53, 23 Their tuning procedure is based on the idea of minimizing the single-shot $G_{0}W_{0}$ correction to the HOMO energy of a DFT calculation. Following Eq. (4), the HOMO energy of a system is strictly equal to the IP for the exact functional, meaning the DFT calculation describes the systems’ IP accurately and the $G_{0}W_{0}$ correction should vanish. This is, however, not true for common functionals, where the correction found from $G_{0}W_{0}$ can be significant.18, 8, 23 By adapting the fraction of exact exchange in the DFT starting point, $G_{0}W_{0}$ tuning aims to minimize the self-energy $G_{0}W_{0}$ correction to the HOMO IP, i.e.,

[TABLE]

where $\epsilon^{\alpha,\mathrm{DFT}}_{\mathrm{HOMO}}$ is the DFT eigenvalue and $\epsilon^{\alpha,\mathrm{QP}}_{\mathrm{HOMO}}$ is the $G_{0}W_{0}$ quasiparticle energy of the HOMO. When compared to the other tuning methods, one benefit of $G_{0}W_{0}$ tuning is that it only requires calculations on the neutral molecule, not the cation. Similar to the DSLE tuning, it requires several $G_{0}W_{0}$ calculations to find the tuned $\alpha$ value. In contrast, IP-tuning only requires computationally more efficient DFT calculations for determining $\alpha$.

Different from the other tuning methods, the fraction of exact exchange found for polyacetylene oligomers using $G_{0}W_{0}$ tuning shows no clear trend with oligomer length, see Fig. 3. Additionally, the tuned $\alpha$ values are significantly lager than for the other tuning methods. While a very slight fluctuation between $\alpha$=0.82 and 0.83 can be observed from N=4 to N=10, $\alpha$ stabilizes at 0.81 for N$\geq$12. This points towards the $G_{0}W_{0}$ tuned functionals being approximately size-consistent, at least for the particular case of the polyene oligomers studied in this work.

Oligomer approach and size dependence

It has been shown that IP tuning faces difficulties when used to calculate properties of polyacetylene oligomers of increasing length. As the spatial extension of the $\pi$-conjugated system increases, the amount of exact exchange required to minimize the error in Eq. (6) falls off considerably.41, 43 This has important consequences for the calculation of structural, electronic, or optical properties. This is exemplified by the example of polyacetylene, where the amount of BLA strongly depends on the amount of exact exchange present in a functional used for the structure optimization.38, 39 While the extrapolation of CCSD(T) results predicts convergence of the BLA to $\sim$0.08 Å after about 8-10 repeat units, IP-tuned functionals show a continued dropping off for increasing oligomer lengths, resulting in grossly underestimated BLA for long chains.38, 40 A similarly erroneous behavior is encountered for optical excitation energies, which were also shown to depend heavily on the amount of HF exchange, independently of any geometry effects.39

To understand the origin of this problem, it is helpful to consider the change of the systems energy upon addition or removal of an electron in the limit of an infinite chain. Semilocal functionals exhibit a strong delocalization error, which means it is energetically favorable for an electron (or hole) to delocalize over the entire $\pi$-conjugated system. For an infinitely long molecular chain this means the electron would delocalize over an infinite number of atoms, such that only an infinitesimally small fraction of an electron is added to each repeat unit. As a consequence, the energy no longer curves between $M$ and $M\pm 1$ but instead behaves like the initial slope of the curve for all fractions of electron added/removed, resulting in the linearity of the total energy between the $M$ and $M\pm 1$ electron systems.77, 44 It should be noted, however, that although the MESIE vanishes, a delocalization error is still very much present but simply a lot harder to quantify. Another consequence is that the IP-theorem no longer bears significant meaning, as every approximate exchange-correlation functional will give a straight line between $M$ and $M-1$ and, consequently, obey the IP-Theorem. Still, each functional might yield a considerably different IP. Consequently, the IP-tuning procedure can not straightforwardly be applied to infinite systems.

While this argument only strictly holds for infinitely large molecular chains, it has been shown that it has important consequences also for finite molecular chains of increasing size. As the system size increases, the delocalization error found for semilocal DFT steadily decreases as the energy curve for fractional occupation numbers approaches linearity.77, 44 Fig. 4 shows that this effect becomes already apparent for relatively short polyacetylene oligomers with up to 16 double bonds. The inset showing the deviation from linearity for each chain length highlights a substantial reduction in MESIE as the chain length increases.

In contrast to the case of (semi-)local and hybrid functionals with no or only a small fraction of HF exchange, this behavior is not observed when using full HF or hybrid functionals that include a very large amount of HF exchange. As soon as the curvature of the energy becomes concave, the localization error becomes effective (see Fig. 4), meaning it is favorable for the added/removed electron to stay localized. In this situation, increasing the system size has no effect on the shape of the energy curve.77

The appearance of piecewise linearity for large finite systems and the associated loss of the physical relevance of the IP tuning scheme are responsible for the lack of size consistency and the related size dependence and consequential failure of IP tuned functionals to accurately describe the structural, electronic, and optical properties polyacetylene oligomers of increasing length. Our aim is to investigate to which extent these problems are carried forward when using IP tuned functionals as a starting point for $G_{0}W_{0}$. As DSLE tuning is also based on enforcing straight line behavior, it can be expected that it would face similar difficulties. A strong indication for this is given by the size dependence of the fraction of exact exchange found using DSLE tuning. In contrast to the two tuning methods enforcing straight line behavior, $G_{0}W_{0}$ shows no size dependence. How this affects its ability to accurately describe the evolution of IP with increasing polyacetylene chain length remains to be seen.

3 Computational Details

Ground state geometry optimizations of polyacetylene oligomers with N=4 to N=22 (see Fig. 1) were performed in Gaussian0978 using a 6-311G basis set, a (75,302) grid, and PBEh* (45.73% exact exchange), a functional that was specifically parameterized to give excellent results for the BLA of polyacetylene and its evolution with oligomer length.38 Butadiene (N=2) was excluded from the calculations, as it has been shown to posses very different properties from the other oligomers due to its shortness.79 The resulting geometries all show $C_{2h}$ symmetry. The subsequent calculations were carried out using these geometries with no further structure optimization to ensure accurate BLA. The total energy curves for fractional electron numbers used to show the appearance of piecewise linearity were carried out in PSI480 using PBE and HF with a cc-pVTZ basis set. All tuning was done by adapting $\alpha$ in PBEh, considering two decimal places. Details on all of the tuning procedures can be found in the background section. IP, $G_{0}W_{0}$ and DSLE tunings, as well as $G_{0}W_{0}$ corrections and most $\Delta$SCF calculations were carried out using FHI-AIMS81, with a tier 4 basis set for all calculations including $G_{0}W_{0}$ and a tier 2 basis set for all others. The $G_{0}W_{0}$ integration was done following the default two pole model, a comparison with the converged Pade model for the shortest chainlengths can be found in the supporting information. For the cation calculations needed for $\Delta$SCF energies for functionals including a large amount of exact exchange ($G_{0}W_{0}$ tuned and HF), an increasing amount of spin contamination with chain length was found using unrestricted Hartree-Fock (UHF), leading to complete failure to converge for longer chains. In these cases we employed restricted open-shell Hartree-Fock (ROHF) calculations in Gaussian for both the neutral and the cation systems for all chain lengths, using a 6-311G(d) basis set.

Basis-set extrapolated CCSD(T) is generally considered to be a reliable reference for vertical IPs.76, 10, 28 Due to the very high computational cost and memory requirements, however, CCSD(T) IPs could only be calculated for Butadiene (N=2) and Octatetraene (N=4). An alternative method called domain based local pair-natural orbital coupled-cluster (DLPNO-CCSD(T) has recently been developed and improved to obtain linear scaling with system size, while reproducing the CCSD(T) total energy with very high accuracy.82, 83, 84 For the case of the polyacetylene oligomers studied here, we found the deviation of the DLPNO-CCSD(T) IPs from the full, basis-set extrapolated CCSD(T) IPs to be 0.02% for N=2 and 0.05% for N=4. Based on this excellent agreement, we have therefore chosen to use DLPNO-CCSD(T) IPs as reference data for all chain lengths up to N=22. All CCSD(T) and DLPNO-CCSD(T) calculations were performed using ORCA.85 An unrestricted reference wavefunction was used for both the neutral and cationic system to obtain $\Delta$SCF IPs. All calculations were performed within the resolution of identity combined with a chain of spheres exchange (RIJCOSX)86 approximation, using ORCAs automatic basis set limit extrapolation. For CCSD(T) N=2 and DLPNO-CCSD(T) N=2-16 the extrapolation was based on cc-pVTZ/QZ87, for CCSD(T) N=4 and DLPNO-CCSD(T) N=22 the extrapolation was based on cc-pVDZ/TZ due to memory constraints. To refer to this (minor) difference in how the numbers were obtained, the DLPNO-CCSD(T) data point for N=22 has been marked with a star in all figures. The DLPNO-CCSD(T) IPs for all chainlengths can be found in the supporting information.

4 Results

A summary of the (negative) DFT HOMO eigenvalues and $G_{0}W_{0}$ HOMO quasiparticle energies found using each of the tuning methods is provided in Fig. 5, along with the reference data obtained from DLPNO-CCSD(T). For the $G_{0}W_{0}$ tuned functional only one set of results is given as, by construction, the DFT and $G_{0}W_{0}$ eigenvalues are equal. For short oligomers, all $G_{0}W_{0}$ calculated IPs as well as the HOMO eigenvalue from the IP-tuned functional show good agreement with the DLPNO-CCSD(T) reference. Only the HOMO eigenvalue obtained from the DFT calculation underlaying the DSLE tuning shows poor agreement with the reference data, severely underestimating IPs by more than 1 eV. Given the small fractions of HF exchange included in these functionals (see Fig. 3), this result is, of course, not surprising. Importantly, the DSLE tuning procedure focuses on reducing the DSLE after the $G_{0}W_{0}$ correction; the DFT calculation merely provides a starting point and, therefore, is not meant to yield accurate HOMO eigenvalues.

At medium and longer chain lengths, a clear difference between the performance of the different tuning procedures arises. The $G_{0}W_{0}$ tuned IPs fall off too slowly, resulting in overestimated IPs with an error increasing with chain length. While the DSLE tuned $G_{0}W_{0}$ results show good agreement with the reference for very short chains, the IPs fall off too fast and an increasing underestimation of the IP with increasing chain length can be observed. IP tuning by itself underestimates the IPs, even for the shorter chain lengths, and slopes off too fast, resulting in an increasing deviation from the reference data for longer chains. This result is a direct consequence of the size dependence of the IP-tuning procedure and resembles the results previously obtained for the BLA38 and optical band gaps43 of polyacetylene. This underlines once again that, while IP-tuning has been demonstrated to yield excellent results for the IPs and EAs of small to medium-sized organic $\pi$-conjugated molecules, great care has to be taken when applying it to large, conjugated molecules, in particular in the context of the oligomer approach. Surprisingly, however, the IPs obtained from applying a $G_{0}W_{0}$ correction to the IP-tuned eigenvalues reproduce the reference data extremely accurately, showing both good agreement for absolute values as well as the slope with increasing chain length.

IPs can also be calculated as the difference in total energy between the neutral and cationic system, known as the $\Delta$SCF approach. While this should be equal to the HOMO eigenvalue for the exact functional according to Eq. (4), there can be a significant difference between the HOMO eigenvalue and the $\Delta$SCF IP for approximate functionals. For the infinite chain limit, however, the HOMO eigenvalue equals its corresponding $\Delta$SCF IP for all functionals that show a delocalization error, as discussed above. Therefore, it is interesting for us to also study the evolution of the $\Delta$SCF IP with increasing chain length. These results are provided in Fig. 6.88

IP-tuned and DSLE tuned functionals lead to very similar $\Delta$SCF IPs, both underestimating the reference data and falling off too rapidly with increasing chain length. This error can be traced back to their size dependence: The strong decrease in exact exchange with increasing chain length leads to increasingly inaccurate $\Delta$SCF IPs. Since the the IP-tuning procedure forces the HOMO eigenvalue to equal the $\Delta$SCF IPs, this result further illustrates why the HOMO IP obtained from IP-tuned functionals becomes more and more inaccurate for $\pi$-conjugated molecular chains of increasing length. In contrast to IP- and DSLE tuning, $G_{0}W_{0}$ tuning results in a very good agreement with the reference, with only minor deviations from both the absolute IP values and their slope. We attribute this to the approximate size-consistency of the $G_{0}W_{0}$ tuning procedure, resulting in an almost constant amount of exact exchange and, thus, allowing the $\Delta$SCF IPs to evolve correctly with increasing chain length.

5 Further Discussion

To interpret the different tuning methods’ behavior with increasing chain length and to highlight their respective strengths and shortcomings, it is helpful to study both their DSLE and the difference between the $\Delta$SCF IP and the negative HOMO eigenvalues, as shown in Fig. 7 on the left and right side, respectively. These two values represent the errors minimized by DSLE and IP tuning, respectively. A clearer comparison between the two quite similar errors can be drawn following the approach by Bruneval59 and plotting half of the DSLE as described by equation 2. This follows from the fact that the energy curve for fractional particle numbers can be expanded as a second order polynomial and hence the error can be expressed as the mean value of HOMO and LUMO eigenvalues.

PBE and HF have been included to provide the limiting cases of 0% HF and 100% HF, as all our tuned functionals are a combination of these two. PBE is a semi-local functional, with a strong delocalization error, represented by a positive DSLE and a $\Delta$SCF IP significantly larger than the negative HOMO eigenvalue. Both errors, however, show a strong decrease with chain length. This behavior of PBE can be rationalized by the appearance of piecewise linearity with increasing system size demonstrated in Fig. 4. HF on the other hand shows a large negative DSLE and $\Delta$SCF IP smaller than the negative HOMO eigenvalue, with almost no change in error with increasing chain length. This shows the localization error, which prevents the appearance of piecewise linearity for increasing system size, hence leading to a consistent performance of the functional for all oligomer lengths. The difference between the 0.5 * DSLE and the $(\Delta$SCF - $\epsilon_{\mathrm{HOMO}})$ seen in figure 7 can be explained by the behavior of the cation LUMO in the ROHF approach, as can be seen in the supporting information.

The largest fraction of exact exchange is found for the $G_{0}W_{0}$ tuned functional, where $\alpha\approx 0.81$. Fig. 7 shows a similiar behavior of the $G_{0}W_{0}$ tuned and HF methods, with both approaches displaying a negative DSLE, a $\Delta$SCF IP smaller than the negative HOMO eigenvalue and almost constant errors with chain length. It follows that $G_{0}W_{0}$ tuned functionals exhibit a localization error. This is also reflected in the density difference of the neutral and cationic systems obtained from the $G_{0}W_{0}$ tuned functional when compared to the DLPNO-CCSD(T) result (see supporting information). The tendency of $G_{0}W_{0}$ to exhibit a localization error for DFT starting points with a large amount of exact exchange has been previously discussed in literature.52

The small discrepancies found for the DSLE of the $G_{0}W_{0}$ tuned functional before and after the $G_{0}W_{0}$ correction can be traced back to a change in the cation LUMO, as only the HOMO is tuned to give identical DFT and $G_{0}W_{0}$ HOMO energies. A figure showing how the $G_{0}W_{0}$ correction affects the 5 highest occupied and 5 lowest unoccupied orbitals of the neutral system can be found in the supporting information.

DSLE tuning leads to the smallest fractions of exact exchange (0.2 $<\alpha<$ 0.4), and the corresponding DFT starting point shows a behavior similar to PBE. Again, a strong delocalization error is apparent through a positive DSLE and a $\Delta$SCF IP much larger than the negative HOMO eigenvalue. Again, this delocalization error can be directly observed in the density difference between the neutral and cationic systems (see supporting information). Even though a small amount of exact exchange is present, the appearance of piecewise linearity for long chains can be observed, the errors seemingly decreasing with chain length. By construction, the tuning minimizes the DSLE after the $G_{0}W_{0}$ correction, which is clearly seen in Fig. 7 (left). As the DFT starting point appears to have decreasing DSLE with increasing chain length, however, the correction needed from $G_{0}W_{0}$ to enforce linearity decreases (see Fig. 8).

For IP tuned functionals linearity is already enforced approximately in the DFT starting point, and the $\Delta$SCF and eigenvalue IPs are equal by construction. The DSLE is almost zero, showing a slight deviation with increasing chain length caused by the behavior of the cation LUMO (see supporting information for details). Applying the $G_{0}W_{0}$ correction on top of IP tuned functionals leads to a slightly negative DSLE indicating a localization error. However, the $G_{0}W_{0}$ correction increases with chain length, shown in Fig. 8, thus working against the error in the DFT starting point introduced by the appearance of piecewise linearity.

6 Conclusions

In summary, we conclude that the ability of $G_{0}W_{0}$ calculations to accurately describe the evolution of the HOMO IP with chain length for $\pi$-conjugated molecular chains such as polyacetylene relies on a number of factors:

Firstly, the amount of exact exchange present in a hybrid functional determines the slope of the IP with increasing chain length, with low fractions of exact exchange causing the IPs to fall off too rapidly and high fractions of exact exchange causing the IPs to not fall off fast enough. This corresponds to similar trends observable for other properties such as BLA and excitation energies.41, 39 Eliminating the localization/delocalization error in either the DFT or the $G_{0}W_{0}$ step (as resulting from IP and DSLE tuning, respectively) does not guarantee accurate IPs or a correct evolution with chain length.

Secondly, $G_{0}W_{0}$ tuning leads to a slight localization error, observable in the eigenvalue IPs. In contrast to the other methods it is, however, not strongly size dependent and hence leads to $\Delta$SCF IPs in good agreement with the reference data.

Finally, a $G_{0}W_{0}$ correction improves the DFT HOMO IPs for all tuning methods. The best agreement with the CCSD(T) reference data, however, is found when using an IP tuned starting point, yielding highly accurate IPs for all oligomer sizes. Given the strong size dependence of the IP-tuned starting point and, related to that, the decreasing accuracy of the IP-tuned HOMO eigenvalues for extended $\pi$-conjugated systems, this result could not be expected. However, it seems that the $G_{0}W_{0}$ correction can counteract the error introduced by the size dependence, at least for the example of polyacetylene studied in this work. This is likely to be caused by a favourable cancelation of errors between the IP-tuned starting point and the $GW$ approximation. We conclude that, although there are known shortcomings of the IP-tuning approach for large $\pi$-conjugated molecules due to its size dependence, this apparently does not diminish the suitability of such functionals as a good starting point for $G_{0}W_{0}$ calculations on large, $\pi$-conjugated molecules.

{suppinfo}

The supporting information includes the IPs obtained from the DLPNO-CCSD(T) calculations, as well as additional data showing the difference in DSLE obtained from fractional occupation curves and HOMO/LUMO eigenvalues, and data comparing the $G_{0}W_{0}$ eigenvalues of the HOMO and LUMO for N=4 and 6 found following the two pole model for $GW$ integration compared to the converged Pade results. Also included are figures showing the $G_{0}W_{0}$ correction to the five highest occupied as well as the five lowest unoccupied molecular orbitals following the $G_{0}W_{0}$ tuning procedure, and the straight line error of the cation LUMO and neutral system HOMO for ROHF and after IP tuning. Finally, density plots for the difference between the neutral and the cationic system as well as for the neutral HOMO are included for all tuning methods and the DLPNO-CCSD(T) reference.

{acknowledgement}

The authors thank L. Gallandi and P. Rinke for helpful discussions.