Baixe Describing Both Dispersion Interactions and Electronic e outras Notas de estudo em PDF para Engenharia Elétrica, somente na Docsity! Describing Both Dispersion Interactions and Electronic Structure Using Density Functional Theory: The Case of Metal-Phthalocyanine Dimers Noa Marom,† Alexandre Tkatchenko,‡ Matthias Scheffler,‡ and Leeor Kronik*,† Department of Materials and Interfaces, Weizmann Institute of Science, RehoVoth 76100, Israel, and Fritz-Haber-Institut der Max-Planck-Gesellschaft, Faradayweg 4-6, 14195 Berlin, Germany Received August 7, 2009 Abstract: Noncovalent interactions, of which London dispersion is an important special case, are essential to many fields of chemistry. However, treatment of London dispersion is inherently outside the reach of (semi)local approximations to the exchange-correlation functional as well as of conventional hybrid density functionals based on semilocal correlation. Here, we offer an approach that provides a treatment of both dispersive interactions and the electronic structure within a computationally tractable scheme. The approach is based on adding the leading interatomic London dispersion term via pairwise ion-ion interactions to a suitably chosen nonempirical hybrid functional, with the dispersion coefficients and van der Waals radii determined from first-principles using the recently proposed “TS-vdW” scheme (Tkatchenko, A.; Scheffler, M. Phys. Rev. Lett. 2009, 102, 073005). This is demonstrated via the important special case of weakly bound metal-phthalocyanine dimers. The performance of our approach is additionally compared to that of the semiempirical M06 functional. We find that both the PBE-hybrid+vdW functional and the M06 functional predict the electronic structure and the equilibrium geometry well, but with significant differences in the binding energy and in their asymptotic behavior. 1. Introduction Noncovalent interactions, of which London dispersion is an important special case, are essential to many fields of chemistry. Such interactions possess a significant component of electrostatic attraction between permanent or instantaneous dipoles and higher order multipoles and dominate in regions where there is little overlap of charge densities, i.e., at medium-range to long-range, as compared to the short-range chemical bond. In principle, exact density functional theory (DFT) should include accurate treatment of the long-range correlation, which is essential for describing noncovalent interactions.1 However, van der Waals (vdW) interactions (a term that we use here interchangeably with London dispersion) are inherently outside the reach of (semi)local approximations to the exchange-correlation (xc) functional as well as of conventional hybrid functionals, based on semilocal correlation.1,2 Many strategies toward inclusion of van der Waals interactions in DFT calculations, at various levels of ap- proximation, have been proposed. Many of those can be divided into three broad categories: (1) nonempirical meth- ods, typically relying on the adiabatic connection theorem,3 wherethelong-rangecorrelationiseithercomputedexplicitly4–11 or integrated with traditional xc functionals;12,13 (2) semiem- pirically parametrized xc functionals, calibrated for data sets that include noncovalently interacting systems;14–18 (3) pairwise addition of C6/R6 corrections to the internuclear energy expression, damped in the short-range while providing the desired long-range asymptotic behavior.19–28 Such C6/ R6 corrections are usually semiempirical but can be derived from first-principles considerations.28 Understandably, most of the literature on DFT computa- tions of dispersively bound systems has focused on obtaining * Corresponding author phone: +972-8-934-4993; e-mail: leeor.kronik@weizmann.ac.il. † Weizmann Institute of Science. ‡ Fritz-Haber-Institut der Max-Planck-Gesellschaft. J. Chem. Theory Comput. XXXX, xxx, 000 A 10.1021/ct900410j  XXXX American Chemical Society correct geometries and binding energies. There are very important classes of systems, however, for which it is crucial to obtain a correct prediction of the electronic structure as well. An important example, on which we elaborate here, is that of small-molecule-based organic semiconductors. In such materials, intermolecular interaction in the molecular crystal is typically dispersive (or at least has a significant dispersion component), and geometry predicted using standard func- tionals can be highly inaccurate, as discussed, e.g., in ref 29. At the same time, an accurate description of the electronic structure is essential to understanding the relations between the chemical nature of the constituent molecules and their function in organic electronic devices. A key question, then, is whether one can systematically obtain a sufficiently accurate theoretical treatment of both noncovalent interactions and the electronic structure, within a computationally tractable scheme that is preferably widely applicable and involves as little empiricism as possible. This is challenging because the electronic structure can be very sensitive to the type of functional used. A recurring reason for inadequate treatment of the electronic structure is the presence of self-interaction errors (SIE),30,31 i.e., the spurious Coulomb interaction of an electron with itself in the Hartree term of the Kohn-Sham equation, which is not fully canceled out by approximate expressions for the exchange- correlation term. Local and semilocal functionals, e.g., the local-density approximation (LDA) and various flavors of the generalized gradient approximation (GGA), respectively, often exhibit significant SIE that results in a poor description of the electronic structure of organic molecules and crystals.32,33 Hybrid xc functionals were found to mitigate the effect of the SIE significantly via the inclusion of a fraction of Fock exchange.31–33 Therefore, a desirable scheme would combine a successful description of van der Waals interactions with a hybrid functional based description of the electronic structure.34 This should be possible because the electronic structure is mostly sensitive to exchange and short-range correlation, whereas dispersive interactions mainly affect the total energies and geometries. In principle, such a successful combination may be achieved within each of the three above-discussed strategies for treating van der Waals interactions. The most practical and successful representative of the first strategy (a nonem- pirical method relying on the adiabatic connection theorem) is the “vdW-DFT” functional of Dion et al.13 (see ref 35 for some recent applications). It is based on a GGA (specifically revPBE36) exchange functional, combined with LDA for the local part of the correlation, on top of which the nonlocal correlation component is added. Although this nonlocal correlation can be combined with other functionals, results for, e.g., the binding energy may depend significantly on the underlying “parent” functional.37 Therefore, we will not be discussing this approach here. Currently the most popular representative of the second strategy (semiempirical methods based on hybrid functionals) is the M06 family of function- als,17 a family of meta-GGA functionals (i.e., functionals whose “semi-local” component includes kinetic energy spin- densities, in addition to the spin-densities and their gradi- ents31) with varying fractions of exact exchange. This approach provides some flexibility in the choice of an appropriate functional, an issue elaborated below. However, the correct long-range R-6 behavior is still absent from such functionals even if medium-range noncovalent binding is well-achieved. The third strategy, addition of pairwise C6/ R6 terms to the internuclear energy term, allows for the highest degree of flexibility in choosing independently the appropriate description of the electronic structure, on top of which a suitable dispersion correction is performed. Obviously using C6/R6 corrections is not free from limitations either. First, the approach assumes that nonco- valent interactions have little direct effect on the electron density and affect the system mainly by influencing the equilibrium geometry. Second, screening by the conduction electrons has to be addressed for metallic systems. Third, the short-range damping function may be problematic for the accurate description of short bond lengths. Fourth, Dobson et al.38 have shown that summation over pairwise interactions may result in incorrect asymptotic behavior in certain special cases, e.g., low-dimensional (semi)metallic systems. Here, we examine the degree to which a quantitative treatment of both the electronic structure and the dispersion interactions is achieved in practice. We show that this is indeed possible using the recently presented “TS-vdW” correction scheme,28 in which the leading-order C6 coef- ficients and vdW radii are determined in a first principles manner from the DFT ground-state electron density. These corrections are combined with the GGA of Perdew, Burke, and Ernzerhof (PBE)39 with the one-parameter nonempirical PBE-hybrid (also known as PBEh or PBE0),40 or with the three-parameter semi-empirical hybrid functional B3LYP.41 We compare our results to those obtained from the M06 functional,17 as well as to those obtained from the standard PBE and PBE-hybrid functionals and to pertinent experiments. We have chosen two members of the metal-phthalocyanine (MPc) family as case studies for the above comparison, NiPc and MgPc. MPc’s are highly stable organic semiconductors with a broad range of applications in, e.g., light emitting diodes, solar cells, gas sensors, thin film transistors, and even single molecule devices.42 Specifically, their electronic structure has been shown to be highly sensitive to self- interaction errors.32 Furthermore, it is known that π-π and π-d interactions, which possess a dispersive component and are attributed to nonlocal electron correlations that occur in systems with spatially close-lying π orbitals,43 play an important role in the stacking of molecules in MPc crystals. In transition metal Pc’s, such as NiPc, π-d interactions affect the intermolecular distance in the stack.44 In crystalline MgPc, π-π interactions not only affect the intermolecular distance but also lead to a structural change in the molecular subunit as the Mg atom deviates from the molecular plane and shifts toward the azamethine N of the adjacent molecule (see also Figure 1), so that the basic unit of the MgPc crystal is, in fact, a dimer.45,46 Thus, both NiPc and MgPc provide stringent test cases for a treatment of both geometrical and electronic structure. Here, we calculate the binding energy curves, geometry, and electronic structure of NiPc and MgPc dimers. We find B J. Chem. Theory Comput., Vol. xxx, No. xx, XXXX Marom et al. results of three M06 variants with corresponding one- parameter hybrids66 based on BLYP,58 a semiempirical GGA functional. To examine the role of exchange, in each case the M06-variant result is compared to a BLYP-based hybrid that has the same the fraction of Fock exchange. Two trends are immediately obvious: First, whereas the M06 spectrum agrees well with experiment, M06L and M06-2X yield spectra that do not. This agrees with the recommendation of Zhao and Truhlar. Second, the M06-variant spectra are remarkably similar (though, of course, not identical) to the corresponding BLYP-hybrid ones. This shows that, despite the many additional fitting parameters used in any M06 variant, the dominant factor in determining the electronic structure is the fraction of Fock exchange. In turn, the spectra obtained with BLYP and with BLYP+27% Fock exchange are remarkably similar to previously published spectra (ref 32 and cf. Figure 5 below), obtained with the nonempirical PBE and PBE-hybrid (i.e., 25% Fock exchange) functionals, respectively, further underscoring the dominant role of Fock exchange. Therefore, of the entire M06 family, only M06 is considered hereafter. Interestingly, the leading (HOMO) and second peak of the revPBE spectrum are much closer to the spectra obtained from the hybrid functionals (BLYP+27% Fock exchange and M06) than to those obtained from the semilocal functionals (BLYP and M06L). We have observed a similar behavior for other MPc’s (not shown for brevity). Likely, this is at least partly because the exchange enhancement factor of revPBE was constructed by fitting to exact exchange-only calculations of total atomic energies.67 This compensates to some extent for self-interaction errors and thus improves the fit to experiment in the higher-lying part of the spectrum. However, this comes at the price of distorting the shape of the third and fourth peaks. Because revPBE, while better than other GGAs in this respect, still fails to yield a satisfactory electronic spectrum, we do not discuss it further here. We now turn to the binding energy curves of NiPc and MgPc dimers, shown in Figure 4, obtained using the PBE, PBE-hybrid, B3LYP, and M06 functionals, with and without the C6/R6 correction. Clearly, the uncorrected PBE and PBE- hybrid calculations underestimate considerably the strength of the noncovalent interaction and overestimate the inter- molecular distance in both dimers. The B3LYP calculations reveal no net attraction at all. This is a known tendency of semilocal and conventional hybrid functionals that has been demonstrated repeatedly for various systems (see, e.g., refs 1, 2, 13, 14, 19, 20, 22, 26). For both dimers, M06 significantly improves upon the semilocal and hybrid func- tionals, yielding binding energies that are higher by about 1.0 eV. However, the binding energies obtained with PBE+vdW, PBE-hybrid+vdW, and B3LYP+vdW are higher yet, by ∼0.7 eV as compared to M06. This difference between the TS-vdW corrected results and M06 is larger than the level of accuracy found in recent M06 studies of smaller dispersively bound systems,55,68 likely due to the sheer size of the MPc molecules and the contribution of the π-d Figure 3. Computed NiPc single molecule spectra, calculated with revPBE, different M06-variants, and BLYP-based single- parameter hybrid exchange-correlation functionals. Raw eigen- value data, as well as the same data broadened by a 0.35 eV Gaussian, are shown. This facilitates comparison with the UPS data of Ellis et al.,64 obtained for a 11.8 Å NiPc film at θ ) 70°, also shown in the figure. Figure 4. Binding energy curves, obtained with different exchange-correlation density functionals, for (a) NiPc and (b) MgPc dimers, composed of planar molecules. Electronic Structure Using DFT J. Chem. Theory Comput., Vol. xxx, No. xx, XXXX E interaction. Furthermore, this difference can be traced back to the long-range behavior of the M06 functional, which is essentially the same as that of PBE or PBE-hybrid at intermolecular distances larger than 5.5 Å. For some ap- plications, the latter difference may be practically unimpor- tant if the near-equilibrium region is well-described. How- ever, it is fundamentally important to realize that the hybrid meta-GGA approach does not possess the correct asymptotic behavior. This limitation may manifest itself practically as well, e.g., in systems where a cumulative effect of many long-range dispersive interactions is important. For both dimers, the addition of the C6/R6 correction to the M06 functional recovers the correct long-range behavior and does not affect the equilibrium intermolecular distance. However the binding energy increases by ∼1.0 eV, becoming ∼0.3 eV larger than with PBE+vdW or PBE-hybrid+vdW. We note that the remaining difference may be attributed to the employed damping function. Since M06 already provides considerable attraction at the intermediate range, it may require a different model for the damping function. Without experimental or high-level quantum-chemical data for the binding energy, it is hard to say which functional yields a more accurate binding energy. However, the difference between PBE+vdW or PBE-hybrid+vdW and M06+vdW is significantly reduced (0.3 eV) as compared to the differ- ence between uncorrected M06 and PBE(-hybrid)+vdW (0.7 eV). To understand how well the approaches studied here do at geometry prediction, we have computed the equilibrium intermolecular distances obtained for the NiPc and MgPc dimers, as well as the shift of the Mg atom from the molecular plane for the latter. The computed values, com- pared to experimental data, are given in Table 1 (additional data on single molecule bond lengths and angles are given in the SI). As discussed above, PBE and PBE-hybrid significantly overestimate the equilibrium intermolecular distance of both dimers, whereas M06 and PBE+vdW yield values that are in good agreement with the experimental ones. The distance of the Mg atom from the molecular plane is underestimated by PBE and PBE-hybrid by ∼0.3 Å, whereas it is within ∼0.1 Å from experiment with PBE+vdW and within ∼0.15 Å with M06. Full relaxation was not performed for M06+vdW, since the latter functional is not implemented in FHI-aims. However, on the basis of the binding energy curve of the planar dimer, we expect only minor changes from the uncorrected M06 dimer geometry. Having accounted for dispersive interactions such that the correct equilibrium geometry was obtained, we now return to the electronic structure. Figure 5 shows calculated eigenvalue spectra of the NiPc dimer, as well as the same spectra broadened by convolution with a 0.35 eV Gaussian to simulate the effective experimental resolution of ultraviolet photoemission spectroscopy (UPS). The calculated spectra are compared to the single molecule spectrum calculated with PBE-hybrid, as well as to the thin film UPS data of Ellis et al.,64 also shown in the figure. As expected, the dimer PBE-hybrid spectrum is similar to that of the single molecule, with some level splitting due to the interaction between the two molecules. In previous work, it was shown that for the NiPc single molecule, as well as for other transition metal Pc’s, the PBE functional fails qualitatively, primarily because of underbinding of localized orbitals due to self-interaction errors.32 A similar picture is revealed for the NiPc dimer, where the spectra calculated with the hybrid functionals, PBE-hybrid and M06, Table 1. Intermolecular Distance in the NiPc and MgPc Dimers and Mg Atom Shift out of the Molecular Plane for the Lattera NiPc intermolecular distance [Å] MgPc intermolecular distance [Å] Mg atom shift [Å] expt 3.24 [44] 3.172 (120 K), 3.185 (260 K) [45] 0.613 (120 K), 0.454 (260 K) [45] PBE 4.2 3.79 0.88 PBE-hybrid 4.1 3.73 0.86 M06 3.30 3.30 0.61 PBE+vdW 3.22 3.29 0.56 a Calculated with the PBE, PBE-hybrid, M06, and vdW-corrected PBE functionals, compared to experimental values. Figure 5. NiPc single molecule (orange online) and dimer (blue online) spectra, calculated with different exchange- correlation functionals and broadened by a 0.35 eV Gaussian, compared to the UPS data of Ellis et al.,64 obtained for a 11.8 Å NiPc film at θ ) 70°. The dimer eigenvalues shown are those obtained for the equilibrium geometry specified in Table 1, except for the PBE-hybrid+vdW and B3LYP+vdW eigen- values that were calculated for the geometry obtained with PBE+vdW. F J. Chem. Theory Comput., Vol. xxx, No. xx, XXXX Marom et al. agree with experiment even at the overestimated intermo- lecular distance of 4.1 Å (obtained with the uncorrected PBE- hybrid). Contrary to the hybrid spectra, the PBE spectra are quite different from experiment. An obvious difference from experiment is that the PBE spectrum is “compressed”; i.e., there is a general narrowing of the gaps between peaks and more energy levels are “squeezed” into a given energy window. “Compression” of experimental spectra is a well- known tendency of semilocal functionals, which can be attributed to the comparison of Kohn-Sham eigenvalues with quasiparticle excitation energies.31,69,70 (Note that in a hybrid calculation, unlike in a “true” Kohn-Sham one, one makes use of a nonlocal potential that can mimic the nonlocal self-energy. This may avoid the “compression” problem.31,32) Moreover, in the PBE spectrum there is a spurious peak between the experimentally observed first and second peaks and the subfeatures for the second peak are missing. This PBE distortion of the line shape remains with PBE+vdW geometry, but the PBE-hybrid+vdW and B3LYP+vdW retain the correct electronic structure.47 Figure 6 shows the eigenvalues of the MgPc molecule, obtained with different functionals, together with selected molecular orbitals and their energy positions. In agreement with trends observed for other transition metal Pc’s,32 the PBE spectrum of MgPc also appears to be affected by SIE. The b2g, eu, b1g, and a1g orbitals, localized over the central region of the molecule, are shifted to higher energies compared to the hybrid spectra, leading to a distortion of the PBE spectrum. Figure 7 shows the calculated eigenvalue spectra of the MgPc dimer, as well as the same spectra, broadened by convolution with a 0.35 eV Gaussian to simulate the effective experimental resolution of UPS. Single molecule spectra obtained with PBE and PBE-hybrid are also shown. As expected, the dimer spectra are similar to those of the single molecule, obtained with the same functional, with some level splitting due to the interaction between the two molecules. Similarly to NiPc, the PBE spectra obtained for the PBE and PBE+vdW geometries appear compressed compared to the hybrid spectra. However, the differences in the line shape between the PBE and the hybrid calculations are not as visually obvious for MgPc as they are for NiPc, at least at the broadening level used.71 Still, on the basis of the qualitative differences in molecular orbital ordering shown Figure 6. Energy and ordering of selected MgPc molecular orbitals, calculated with different exchange-correlation functionals. All spectra were shifted so as to align the highest occupied molecular orbital (HOMO). For clarity, only one example of each doubly degenerate eg and eu orbital is shown. Figure 7. MgPc single molecule (orange online) and dimer (blue online) spectra, calculated with different exchange- correlation functionals and broadened by a 0.35 eV Gaussian. The dimer eigenvalues shown are those obtained for the equilibrium geometry specified in Table 1, except for the PBE- hybrid+vdW and B3LYP+vdW eigenvalues that were calcu- lated for the geometry obtained with PBE+vdW. Electronic Structure Using DFT J. Chem. Theory Comput., Vol. xxx, No. xx, XXXX G