**On the Structure and Geometry of Biomolecular Binding**

On the Structure and Geometry of Biomolecular Binding

(Parte **1** de 4)

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 nonempiricalhybridfunctional,withthe dispersioncoefficientsand van der Waalsradiidetermined 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 interactionspossess a significant component of electrostaticattractionbetweenpermanentor 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 approximation, have been proposed. Many of those can be divided into three broad categories: (1) nonempirical methods, typically relying on the adiabatic connection theorem,3 wherethelong-rangecorrelationiseithercomputedexplicitly4–1 or integrated with traditionalxc functionals;12,13 (2) semiempirically parametrized xc functionals, calibrated for data sets that include noncovalently interacting systems;14–18 (3) pairwise addition of C6/R6 corrections to the internuclear energyexpression,dampedin the short-rangewhileproviding 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 computations of dispersivelybound 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. X, x, 0 A

10.1021/ct900410j X 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-basedorganicsemiconductors.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 functionals can be highly inaccurate, as discussed, e.g., in ref 29. At the same time, an accuratedescriptionof 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-interactionerrors (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 exchangecorrelation 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,3 Hybrid xc functionals were found to mitigate the effect of the SIE significantly via the inclusion of a fractionof Fock exchange.31–33Therefore,a desirablescheme 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-rangecorrelation,whereasdispersiveinteractionsmainly 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 nonempirical 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 representativeof the second strategy (semiempiricalmethods based on hybrid functionals) is the M06 family of functionals,17 a family of meta-GGA functionals (i.e., functionals whose “semi-local”component includes kinetic energy spindensities, in addition to the spin-densities and their gradients31) 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 coefficients 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-hybridfunctionalsand to pertinentexperiments.

We havechosentwomembersof themetal-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 selfinteraction 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 transitionmetal Pc’s, such as NiPc, π-d interactionsaffect the intermolecular distance in the stack.4 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. x, No. x, X Marom et al.

that PBE+vdW, PBE-hybrid+vdW, and M06 all yield similar geometries, but the electronic structure is well described only with the PBE-hybrid+vdW, the B3LYP+ vdW,47 and the M06 approaches. Moreover, we find significant differences in the binding energy between PBE- hybrid+vdW and M06. We attribute these differences to the long-rangebehaviorof these two methodsand show that they can be reduced by applying the TS-vdW C6/R6 correction to M06.

2. Methodology a. TS-vdW Correction Scheme. In the TS-vdW28 C6/R6 correction scheme used here, the pairwise vdW interaction,

Edisp, added to the internuclear energy term, is given by

where C6ij is the dispersioncoefficientfor the ij pair of atoms, Rij is the interatomic distance, Rij0 is the sum of equilibrium vdW radii for the pair, and fdamp is a damping function discussed below. The novel feature of the TS-vdW scheme is that the parameters C6ij and Rij0 are determined from first principles. The method yields significantly lower errors for the S22 database of molecular binding energies than empiri- cal C6/R6 methodsand has been recentlyshownto outperform the latter for water clusters.48

Briefly, the TS-vdW scheme is based on accurate ab initio computed reference values for free atom static dipole polarizabilities and C6 coefficients,49 a combination rule for deriving heteronuclear C6 coefficients from homonuclear static dipole polarizabilities, and Hirshfeld partitioning50–52 of the DFT electron density to calculate the relative polarizability of an atom inside a molecule. In this way, different atomic hybridization states are inherently taken into account for different molecular geometries. Complete details are given in ref 28.

The damping function in eq 1 is needed to avoid the divergence of the R-6 term at short distances and reduce the effect of the correction on covalent bonds. A Fermi-type function was used here, in the form where d determines the “steepness” of the damping function and sR reflects the range of interaction covered by the chosen DFT exchange-correlationfunctional.25 By fitting to the S22 database of Jurecka et al.,53 which contains binding energies of 2 differentweaklyboundsystemscloseto CCSD(T)basis set limit, the value of d was set to 20 and sR was set to 0.94 for PBE and 0.96 for the PBE-hybrid.28,54 A similar procedure for B3LYP yielded an sR of 0.84, indicating that a smaller range of dispersion interaction is covered, likely due to a somewhat more repulsiveexchange componentthan that of PBE or the PBE-hybrid. Finally, as discussed by Karton et al.,5 the M06 functional yields attraction at intermediaterange but still does not possess the correct longrange behavior. Therefore, the same fitting procedure was performed for M06 as well, yielding the sR value of 1.16. b. Computational Details. The routines for evaluation of energies and forces using the TS-vdW method have been implementedin the FHI-aims code56 for consistent geometry optimizations.FHI-aimsis an all-electronelectronicstructure code developed at the Fritz Haber Institute. It uses efficient numerical atomic-centered orbitals (NAO) as a basis set and allows one to achievehighly convergedresultswith optimum efficiencyin computerresources.In this work, the tier2 NAO basis set, which yields results that are similar in accuracy to those of the aug-c-pVQZ Gaussian basis set for the S22 database, has been employed for geometry relaxation.

Additional calculations of single molecule geometry and dimer geometry of MgPc and NiPc were performed using the Gaussian57 code with the PBE, PBE-hybrid,B3LYP, and M06 functionals. Calculations of the electronic structure of a single NiPc molecule were performed using the revPBE functional, the functionals M06L and M062X of the M06 family, and BLYP58-based functionals with similar fractions of exactexchange.The Def2-TZVPWeigend-Ahlrichsbasis set59 was used for all calculations, except for the MgPc M06 calculations, for which a larger all-electron c-pVTZ basis60 was used. Throughout this work, the single molecule geometry was optimized independently for each functional and basis set.

Binding energy curves were constructed using the PBE, PBE-hybrid,B3LYP, and M06 functionals,with and without

C6/R6 corrections. The counterpoise (CP) method61,62 was used to correct for basis set superposition errors (BSSE). In order to obtain dimer binding energy curves as a function of a single parameter, the intermolecular distance was varied under the assumption that the metal atom of one molecule lies directly above the azamethine nitrogen of the other molecule,4–46 and that the metal atom is in the molecular plane. The latter assumption is consistent with experimental

Figure 1. Schematic top-view and side-view of the MPc dimer. The metal atom of one molecule lies above the azamethinenitrogenof the othermolecule.In the MgPcdimer, the Mg atom is shifted from the molecular plane, as shown in the side-view.

Electronic Structure Using DFT J. Chem. Theory Comput., Vol. x, No. x, X C

observations for the molecular stacks in crystalline NiPc.4 Therefore, for NiPc the equilibrium geometry was deduced from the minimum of the binding energy curve. A full relaxation with the PBE+vdW functional has indeed confirmed that the monomers remain planar. As noted above, for the MgPc dimer the monomers do not remain planar.45,46 Therefore, for this system a full geometry relaxation has additionally been carried out with all functionals used, so as to obtain realistic geometries.

The basis set convergence of our calculations was verified by direct comparisonof the eigenvaluesand binding energies obtainedfrom PBE calculationsof a MgPc dimer comprising two planar MgPc molecules at an interplanar distance of 4 Å, using both FHI-aims and Gaussian with the basis sets specified above. The two spectra were in good agreement with a maximal difference of 0.0065 eV and a mean error of 0.002 eV (the latter is equivalent to a relative mean error of 0.06%), for all eigenvalues larger than -15 eV. The binding energy obtained using the Gaussian code with the CP-correctedcc-pVTZ basis set was smaller by 20 meV than the value obtained using the FHI-aims code with the tier2 basis set, and smaller by 13 meV than that obtained with the tier3 basis set, which essentially recovers the complete basis set limit.

For additional insights into basis set convergence issues,

Figure 2 shows binding energy curves of the MgPc dimer, calculated with PBE using a smaller, double- (DZ) level basis set, consisting of an all electron c-pVTZ basis set for the Mg atoms and the 6-31G(d,p) basis set for the H,C, and N atoms. Using a DZ basis set leads to an overbinding of 0.25 eV relative to the FHI-aims tier2 basis set and yields an intermolecular distance of 3.7 Å, as compared to 4.0 Å with the tier2 basis set. The CP procedure overcorrects this overbinding and results in an underbinding of 0.06 eV and an intermolecular distance of 4.3 Å. The larger, triple- (TZ) basis set yields a distance of 4.0 Å, in agreement with the FHI-aims tier2 result. As expected, its CP correction is much smaller and reduces the binding energy by 0.05 eV, overcorrecting by only 0.02 eV compared to the tier2 basis. This close agreement between results obtained using different types of basis sets within different codes shows that our results are well-converged. Basis set convergence tests for the M06 functional are discussed in the Supporting Informa- tion (SI). Importantly, we note that reliance on DZ basis sets mayleadto spuriousagreementbetweenM06andPBE+vdW binding energy curves, whereas the CP-corrected TZ basis set calculations reveal pronounced differences between the two curves, which are elaborated below. We note that while our TZ and tier 2 NAO calculations are sufficiently converged, BSSE errors can also be reduced substantially using diffuse functions. We have not utilized this route here because for the large systems studied in this work we have found that this introduces severe convergence difficulties.

3. Results and Discussion

Because we aim at a treatment of both the equilibrium electronicstructureand the long-rangedispersiveinteractions, we start our analysis by choosing which functionals are the most promising candidates for providing such a comprehen- sive treatment. For examining the TS-vdW C6/R6 correction scheme, we focus primarily on PBE and the PBE-hybrid as a prototypical semilocal and hybrid functional, respectively. We are well aware that B3LYP is likely the most popular choice for a hybrid functional but even so prefer the PBE- hybrid for several reasons. First, the PBE-hybridand B3LYP yield essentiallyindistinguishablespectra for metal-phthalocyanines (ref 32 and cf. Figures 5 and 7 below), so we prefer to introduceas little empiricismas possible.This is especially so given that the PBE-hybrid has other advantages over B3LYP, e.g., yielding the correct result for the uniform electron gas limit and doing significantly better at predicting solid state atomization energies.31,63 Furthermore, as noted abovethe rangeof dispersioninteractioncoveredby the PBE- hybridis somewhathigherthan that of B3LYP.Nevertheless, because of its prominence in applications we do provide B3LYP results as well.

Next, we assess, using the NiPc electronicstructure,which of the M06 family of functionals we should pursue. This family, constructed by Zhao and Truhlar,17 consists of four different functionals, denoted as M06 (fractional Fock exchange), M06L (fully semilocal treatment of exchange, i.e., no fractional Fock exchange), M06-2X (with twice as much Fock exchange as in M06), and M06-HF (with 100% Fock exchange). Of those, M06 was recommended by Zhao and Truhlar for systems involving both transition metal chemistry and noncovalent interactions,17 but it is instructive to consider the accuracy of the electronic structure obtained with other functionals of the M06 family. We additionally examine the electronic structure obtained from the revPBE GGA, primarily because the above-discussed “vdW-DFT” functional is based on it.

Figure 3 shows calculated eigenvalue spectra of the NiPc monomer, as well as the same spectra broadened by convolutionwith a 0.35 eV Gaussianto simulatethe effective experimental resolution of ultraviolet photoemission spectroscopy(UPS)experimentsperformedby Elliset al. on NiPc thin films,64 also shown in the figure. We note that, strictly speaking, Kohn-Sham eigenvalues are not equivalent to quasiparticle excitation energies even if the exact xc functional is used.31 Nevertheless, if a suitable approximate xc functional is used, they are often good approximations to electron removal energies.31,33c,65 The figure compares the

Figure 2. Binding energy curves for the MgPc dimer, obtained with the PBE functional using various basis sets. CP denotes use of the counter-poise correction method.

D J. Chem. Theory Comput., Vol. x, No. x, X Marom et al.

results of three M06 variants with corresponding oneparameterhybrids66 basedon BLYP,58 a semiempiricalGGA 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 structureis 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 fromthe hybridfunctionals(BLYP+27% Fockexchangeand 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 intermolecular 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, 2, 26). For both dimers, M06 significantly improves upon the semilocal and hybrid functionals, yielding binding energies that are higher by about 1.0 eV. However, the binding energies obtained with PBE+vdW,PBE-hybrid+vdW,andB3LYP+vdWarehigher yet, by ∼0.7 eV as compared to M06. This difference betweenthe TS-vdWcorrectedresultsand M06 is largerthan the level of accuracy found in recent M06 studies of smaller dispersively bound systems,5,68 likely due to the sheer size of the MPc molecules and the contribution of the π-d

Figure3. ComputedNiPcsinglemoleculespectra,calculated with revPBE, different M06-variants, and BLYP-based singleparameterhybridexchange-correlationfunctionals.Raw eigenvalue 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 1.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. x, No. x, X 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 applications, the latter difference may be practically unimportant if the near-equilibrium region is well-described. However, 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.

(Parte **1** de 4)