Density functional theory approach to gold-ligand interactions: Separating true effects from artifacts

Interactions of gold clusters and surfaces with ligands are of great interest to chemistry for several reasons. These interactions are the key to the self-assembly of ligand-protected gold nanoparticles.^1,2 Ligand-protected gold nanoparticles are intensely investigated for their potential uses in sensing and plasmonics.³ These interactions are also crucial to the unique catalytic properties of gold: e.g., in heterogeneous catalysis by nano-sized gold⁴ and in homogeneous catalysis by gold compounds.⁵ Finally, there is a need for the fundamental understanding of the rapidly developing field of the coordination chemistry of gold.^6,7 Study of gold-ligand interactions are the key to this understanding.

Gold-ligand interactions are examples of donor-acceptor (DA) interactions, sometimes referred to as charge transfer (CT) interactions, where gold in many oxidation states acts as Lewis acid and ligands mainly act as Lewis bases. Of course, it is possible for these roles to reverse.⁸ 

From the outset these interactions are surrounded by controversy regarding the necessity for and the meaning of the charge-transfer term in the partitioning of the interaction energy into physically meaningful components. For example, the symmetry-adapted perturbation theory (SAPT) has “no room” for the CT term⁹ although there have been recent attempts to make such an accommodation.^10,11 Other partitioning schemes^12,13 distinguish between the polarization and charge-transfer contributions (the latter are sometimes referred to as orbital contributions). In the Reed and Weinhold scheme the CT terms are larger than the total interaction energies.¹⁴ 

There is considerable literature evidence that conventional density functional theory (DFT) methods have difficulties in dealing with DA interactions be it in main-group or in transition metal chemistry. A variety of causes for these problems have been presented in the literature. To name just a few, the Yang group pointed to localization error in common functionals as the reason behind large errors in thermochemistry of Diels-Alder cyclization and dimerization of Al complexes.¹⁵ In a similar vein, Goerigk and Grimme constructed a separate database of test systems where the self-interaction error (SIE) is expected to be the problem.¹⁶ Typical DA complexes ethylene-F₂ and NH₃-ClF figure prominently in this database. In a paradigm photovoltaic system, TTF-TCNQ, Sinni et al. brought to the fore the importance of the energy gaps between HOMO of the electron donor and LUMO of the electron acceptor to interactions and charge transfer in DA complexes.¹⁷ They also identified another source of uncertainty in interaction results – the choice of an exchange functional. HOMO-LUMO gaps narrow considerably in complexes of base-donor and radical-acceptor, thereby exacerbating the problem of erroneous charge transfer in the ground state, as the work of Johnson et al. noted.¹⁸ They attributed this problem to the lack of piecewise linearity of DFT energy of both donor and acceptor with respect to an addition and removal of fractional electrons.¹⁹ The implications of SIE on the charge flow across the metal-ligand DA interface have been extensively discussed^20,21 (and references therein).

Steinmann et al. took a different view by arguing that the reason why DA interactions are challenging for the DFT is that it does not properly account for dispersion.²² While the importance of dispersion is unquestionable in π-π stacking DA interactions,^23,24 in DA complexes involving the σ*-acceptor halogen molecules Kozuch and Martin found that the addition of dispersion was not critical.²⁵ Our own results show that in addition to dispersion it is the many-electron SIE in conventional DFT methods artificially altering the donating and accepting abilities that contribute to the problem.²⁶ 

In such circumstances a typical literature response is to examine dozens of different functionals for a particular type of interaction and recommend those that agree best with a benchmark of some sort. For metal-ligand complexes this strategy often results in some eye-popping spreads of binding energies. For example, Seth et al.²⁷ show wide variations in the binding energy of phosphine to the Grubbs II catalyst (Ru metal center): from 1/10 (BLYP²⁸) to 1/2 (PBE0²⁹) to a bizarre four times (LC-M06L) the experimental estimate. A several-functional compilation for AuCO³⁰ provides binding energies from 1/3 (M06-2X³¹) to 3/2 (PBE0) to over twofold overbinding by the favorites of the gold community: BP86,^32,33 PBE,³⁴ and PW91.³⁵ A similar, albeit narrower, comparison for Au-methylthiol³⁶ shows twofold disparity between PW91 (PBE) and B3LYP.³⁷ Arguably, with such uncertainty on every rung of DFT, the structural and energetical predictions of ligand protected gold nanoparticles could vary from functional to functional. Indeed, Häkkinen et al. observed that PBE and BLYP lead to different structures of Au₃₈(SCH₃)₂₄.³⁸ The predictions for ligand exchange reactions may also suffer. Therefore, benchmarking of these interactions with trustworthy methods is needed.

There are several reasons for the problems with computing gold-ligand DA interactions. The first one concerns a proper description of electron affinity (EA) and ionization potential (IP) of the acceptor and donor molecules.³⁹ Inaccuracies in the computational treatment of these quantities, such as, neglecting relativistic effects or insufficient account of correlation, particularly for gold, could lead to order of magnitude differences in binding energy. The second reason originates from the DFT treatment of these quantities (and is not specific to gold-ligand complexes). In an approximate Kohn-Sham (KS) scheme IP and EA are not related to KS frontier eigenvalues and often the magnitude of HOMO energy of the donor is smaller than its IP and the magnitude of LUMO energy of the acceptor is larger than its EA.⁴⁰ This mismatch alters a flow of charge between a donor and acceptor and poses a serious limitation on the application of DFT to metal-molecule interfaces.²⁰ An equivalent view of this problem involves a delocalization error, that is, the nonlinearity of DFT energy between the integer numbers of electrons.⁴¹ 

The Generalized Kohn-Sham (GKS) formalism involving range-separated hybrid functionals may help ameliorate some of these problems. The presence of nonlocal, orbital dependent exchange via the range-separation helps the HOMO and LUMO eigenvalues to acquire the meaning of IP and EA, respectively, and aids in absorbing the derivative discontinuity.⁴² An optimal choice of the distance at which the long range is to commence may include setting the value of range separation parameter to satisfy the condition that a negative of the HOMO energy (ɛ_HOMO) equals IP⁴³ or evaluating the average distance between an outer electron and its exchange hole.⁴⁴ Such tuning leads to two desirable properties of a functional: HOMO and LUMO of the GKS scheme become the IP and EA, respectively, and the functional gains the piecewise linear dependence with respect to fractional electrons. As the added benefit, the exchange-correlation (xc) potential acquires the proper asymptote. The tuning strategy dramatically improves the quality of CT excitations.⁴³ In this work we intend to show that it could be equally effective for the ground-state interactions involving partial electron transfer from a donor to an acceptor.

Our second aim is to apply the ideas of optimized range-separated functionals to the computations of gold-ligand interactions which are uncontaminated by what Johnson et al. call an artificial charge transfer.¹⁸ Our models are sufficiently small to be validated by coupled-cluster singles, doubles, and perturbative triples (CCSD(T)) computations and include the interaction of Au₄ with common ligands such as pyridine, trimethylphosphine, SCN⁻, and benzenethiolate and benzenethiol. They are representative of species used in the functionalization of gold clusters and surfaces. The aurophilic interactions,⁶ important ingredients of the nanoparticle protection, will also be revisited in this context.

There exists substantial DFT (see Refs. 38 and 45–48 and references therein) and ab initio^49,50 literature on the subject of these and related interactions. The comparison of our results with these previous data will be discussed below.

The functionals used in the present work employ the range separation of the general form:

where erf and erfc stand for the error function and its complement, and ω is the range separation parameter. This expression partitions the total exchange energy into short-range and long-range contributions. α defines the constant fraction of HF exchange in the short range. The functional employed in this work is composed of PBE correlation and the short-range exchange ωPBE based on Henderson, Janesko, and Scuseria (HJS) model of the exchange hole.⁵¹ The choice of the DFT exchange functional is important in the context of DA interactions.¹⁷ Too repulsive an exchange, such as B88,³² makes interaction potentials too shallow, whereas an insufficiently repulsive one, such as PBE, makes interaction potentials artificially deep.^52,53 The exchange functional ωPBE based on HJS exchange hole is one of the most reliable in this regard.⁵⁴ 

A two-parameter tuning⁵⁵ of this functional involves the optimization of ω to minimize the condition that −ɛ_HOMO equals the (vertical) ionization potential⁵⁶ 

for six different values of α = 0.0, 0.1, 0.2, 0.3, 0.4, 0.5. A similar family of optimized functionals was previously proposed in Refs. 57 and 58 recently employed in Ref. 59 Condition (2) imposes Janak's theorem that the exact functional satisfies.⁶⁰ The functionals will be denoted LRC-ωPBEα. The parameter ω, and indeed α as well, are system specific. This may create a problem for two different subsystems with different donor-acceptor abilities. The family of optimized functionals increases the likelihood that there will exist a parameter set common to both subsystems.

The optimized LRC-ωPBEα functionals can also be characterized (and referred to) by the total percent contribution of HF exchange, i.e., originating from both the short and the long-range, in the exchange energy. Such a contribution can be computed for a given geometry by performing three computations: first with the full functional and saving the vectors, then switching off the DFT and the HF exchange, respectively, and performing no-SCF calculations with the saved vectors.

A DFT interaction energy lacks the long-range dispersion component. It was computed through the C₈/R⁸ term from the DFT-D3 approach of Grimme and coworkers⁶¹ and denoted E_disp(8). We used the same D3 parameters as for LC-ωPBE functional⁶² of Vydrov and Scuseria. This choice has been justified in our recent work.⁶³ 

All calculations employed aug-cc-pVTZ basis set⁶⁴ or its equivalent aug-cc-pVTZ-PP, which in the case of Au, combined the 19-electron effective core potential of Figgen et al.⁶⁵ with valence basis set of Peterson and Puzzarini⁶⁶– all denoted avtz. Our CCSD(T) calculations were performed in the same avtz basis set as DFT except for one instance: In the complete basis set extrapolation (CBS) of (AuH-PH₃)₂ where both avtz and aug-cc-pVQZ (avqz) calculations were needed for the expression of Halkier et al.,⁶⁷ 

where E_corr refers to CCSD(T) correlation energy and X = 4. The HF component of the CBS interaction energy was obtained without extrapolation in the avqz basis set.

The interaction energies were computed via the supermolecular method with the counterpoise correction (CP)⁶⁸ while holding the monomer geometries rigid. The monomer geometries of ligands were optimized using the B3LYP functional. The geometry of Au₄ was optimized within LC-ωPBE functional of Vydrov and Scuseria.⁶² The calculations were performed using NWChem⁶⁹ and MOLPRO⁷⁰ program suites.

To test the idea of tuned functionals for their suitability to provide interaction energies we first consider the water dimer since H₂O can be both an electron donor and acceptor. Other motivation stems from the recent observation of Ziegler, Grimme, and coworkers that the employment of range separation in the water dimer can lead to erratic results in the computed interaction energies.²⁷ 

Table I shows the variations of the interaction energy with the optimized LRC-ωPBEα functionals. E_int(DFT+D3) values are very reasonable given the consensus value is −5.0 kcal/mol.⁷² The variations in E_int(DFT) are very small in a wide range of the parameter space, e.g., its value varies only by 0.05 kcal/mol between α = 0, ω = 0.503 and α = 0.3, ω = 0.351. Furthermore, E_int(DFT) remains virtually independent of the total percentage of HF exchange (%HF) up to about 50%. Above that, the interaction energy begins to drop. This drop at larger %HF and small values of ω may have two causes: (i) 1/ω, i.e., the range of interaction described by the short-range exchange, becomes comparable with the length of the hydrogen bridge H—O; (ii) too high percentage of HF exchange may upset the balance between exchange and correlation within the functional. The main conclusion of this comparison is that with the help of optimized LRC-ωPBEα functionals we can obtain interaction energies which are nearly constant in a broad range of α and ω parameter space. Such interaction energies are free from delocalization error and remain physically sensible as long as the percentage of HF exchange (or 1/ω) is not too large.

To determine how these ideas apply to gold, as our second test case, we examine an aurophilic dimer of two HAu-PH₃ subunits arranged in a configuration that precludes electrostatic dipole-dipole interactions (Fig. 1(a)). There has been great deal of discussion in the literature about the origin of the aurophilic attraction attraction,^6,73–77 and to what degree DFT can be trusted in reproducing it. Other motivation follows from the work of Werner and collaborators who demonstrated the importance of ionic charge-transfer-type configurations in aurophilic bonding.⁷⁸ 

The interaction energies for the aurophilic dimer are shown in Table II for the six optimized functionals. The E_int values evaluated at the repulsive wall (R = 2.5 Å), near the minimum (R = 3.0 Å) and in the long range (R = 5.0 Å) remain nearly constant, to within 0.1 kcal/mol in the minimum, and to 3% on the repulsive wall. Let us recall that at the HF level the interaction is purely repulsive.⁷³ The present results show that the interaction energy, free from delocalization and BSSE errors, is attractive thus indicating the presence of a short-range intermolecular correlation contribution. The interaction energy still lacks the long-range dispersion contribution. This is estimated from D3 method of Grimme⁶¹ to be −4.80 kcal/mol bringing the total binding energy to a reasonably good agreement with the CCSD(T) benchmark value of −6.13 kcal/mol. By comparison, Johnson's value of DFT+disp equals −5.59 kcal/mol.⁷⁷ 

The frontier orbital energies of the HAu-PH₃ monomer and the HOMO-LUMO gap (E_LH), shown as a function of %HF in the optimized functionals, remain constant (see Fig. S1 in supplementary material). Upon the interaction, E_LH narrows somewhat (by ∼1.3 eV at the minimum) and it too remains fairly constant as %HF increases (see Table II).

The parameters for the six optimized LRC-ωPBEα functionals for Au₄ and five ligands are presented in Table III. Although the optimized parameters for Au₄ and the ligands are not the same, as required by the supermolecular method, there exists in each case an α and ω combination which is common or nearly so for Au₄ and each of the ligands.

The dependence of these functionals on fractional numbers of electrons in Au₄ is shown in Fig. 2. Functionals PBE0 and BP86 are included for comparison. It is seen that the optimized functionals show practically indistinguishable performance.

Figure 3 displays the eigenvalue variations of Au₄ (top) and SCN⁻ (bottom) with optimized functionals. The two sets present an interesting contrast between molecules consisting of transitions metals and these of main-group atoms. In Au₄ only the three upper orbitals appear distinct while the lower, chiefly d-composed eigenvalues, form a band-like structure. By contrast, in SCN⁻ the occupied orbitals are well separated (apart from degeneracy). In SCN⁻ several occupied orbitals remain constant, whereas in Au₄ only the frontier orbitals remain so and the inner ones slope down as the functionals acquire more HF exchange (see. Refs. 57 and 58 for similar observations). The frontier orbitals, HOMO of the donor and the LUMO of the acceptor, are separated by a very small energy gap of ∼1.5 eV, so one expects a strong interaction.

SCN⁻ ligand is a soft base which binds to a soft acid via its S end. The orientation of the monomers in the complex is determined by mutual alignment of high-lying HOMO of SCN⁻ (π* maximally extended at S) and low-lying LUMO of Au₄ (see Figs. 1(b) and 1(c)). Consequently, the SCN⁻ fragment and the other sulfur ligands are nearly perpendicular to the plane of Au₄. The interaction energies are shown in Table IV. The monomers were held fixed at their geometries and only the intersystem parameters R and θ were roughly optimized,⁸⁰ where R is the Au-S distance and θ denotes the angle between the shorter axis of Au₄ and S-C.

For the Au₄-SCN⁻ interaction, Table IV shows the results for three LRC-ωPBEα functionals where the parameters for the donor and acceptor roughly overlap. The purpose of this comparison is to demonstrate that the results are fairly independent of the choice of optimized functional. Indeed, the binding energies remain stable with respect to the functional, similarly to the case of the water dimer (Table I). The interaction of this strength induces changes in electron density which manifests themselves in the opening of the HOMO-LUMO gap and a partial charge transfer from SCN⁻ to Au₄ of 0.44e, according to Mulliken population analysis.

Let us assume for a moment that one electron is fully transferred to Au₄. In Au₄⁻ one can isolate the interaction energy between one Au (²S) atom on the shorter diagonal with the Au₃⁻, all considered in the geometry of Au₄. Such E_int computed with TPSS/avtz amounts to −47.3 kcal/mol, i.e., the value closely resembling the interaction energy between Au₄ and SCN⁻ in Table IV. The binding of S-C₆H₅⁻ is some 10 kcal/mol stronger still. One can expect that such an interaction is capable of pulling out an Au atom from the Au₄ cluster. Indeed, as the geometry of Au₄-SC₆H₅⁻ is allowed to fully relax, the complex begins to split into Au₃⁻ and Au-SC₆H₅ subunits (see Fig. 4). The ability of thiolates to etch Au atoms³⁸ or even a monoatomic nanowire⁸¹ from Au_n clusters has been shown in the literature.

The neutral benzenethiol (PhSH) provides interesting context for the interactions of thiolates. The bonding between PhSH and Au₄ in two optimized functionals shows to be of intermediate strength (23.7–23.4 kcal/mol) with 4.2 kcal/mol coming from the D3 dispersion. This is again a complex in which thiol acts as an electron donor and gold as an acceptor. Although PhSH has been observed adsorbed on gold surface,⁸² there is also a prevailing view that S-containing ligands are electron acceptors from gold (see, e.g., Ref. 38 and 83). Let us examine if the full electron transfer is possible from Au₄ to PhSH. The energy gap to overcome would involve LUMO(thiol)-HOMO(Au₄) of 0.4 eV + 7.6 eV = 8 eV. This energy gap is certainly larger than LUMO(Au₄)-HOMO(thiol) of −1.9 eV + 8.4 eV = 6.5 eV. The energy balance in the first scenario can be partially offset by 5.8 eV originating from the interaction of ion pair separated by 2.48 Å, but this still leaves a 2.2 eV deficit. This scenario also calls for the formation of metastable PhSH⁻ in a dissociative state PhS⁻ + H. We conclude that in the Au₄ + PhSH complex the thiol molecule can only act as an electron donor. To make it an electron acceptor would require a gold cluster with much lower IP (e.g., Au_n, n-odd, have lower IPs) or the dissociation of S-H. The latter issue has been extensively discussed in the literature.^{36,45,47,46,82}

A less convoluted explanation takes into account that electrons tend to flow from a region of high chemical potential (low electronegativity) to a region of low chemical potential (high electronegativity). Both of these quantities are well defined within DFT.⁸⁴ The electronegativity χ = ½ (IP+EA) = −μ can be described in the context of GKS using the HOMO and LUMO eigenvalues. The values listed above yield χ(Au₄) = 4.7 eV and χ(PhSH) = 4.4 eV, thus, indicating the direction of charge flow from PhSH to Au₄.

The remaining results in Table IV refer to the complexes of Au₄ with pyridine and trimethylphosphine. θ represents the angle between the short axis of Au₄ and the C₂(C₃) axes of the ligands. Pyridine binds to Au₄ in the coplanar fashion and P(CH₃)₃ in a co-axial one. The optimized LRC-ωBPEα functionals for both complexes feature α = 0.2. At this DFT level the interaction energy of pyridine is −26.7 kcal/mol, whereas that of P(CH₃)₃ is 21.6 kcal/mol stronger (Table IV). This difference cannot be explained by the differences in the IP of the donor alone as the IP of pyridine is 9.26 eV⁸⁵ (vertical IP from the present monomer calculations yields 9.72 eV), while the IP of P(CH₃)₃ is 8.12 eV⁸⁵ (vertical IP from our calculations yields 8.33 eV). The energy difference is due to the different character of the two upper-most occupied orbitals of the two complexes.

In Fig. 5 we compare the HOMO-1, HOMO, and LUMO orbitals of both ligands bound to Au₄. The HOMO-1 orbital of Au₄-pyridine shows a node between Au d-orbital and the N lone pair. In Au₄-P(CH₃)₃, the more extended lone pair on P envelops this node to overlap with the same-phase lobes of the d-orbital on the nearest Au. The HOMO orbital of Au₄-P(CH₃)₃ extends further from Au₄ toward P revealing hints of back-bonding.

In both complexes there is a strong enhancement of the dipole moment upon complex formation (4.6 D with pyridine and 4.2 D with P(CH₃)₃). A difference emerges from the Mulliken charge rearrangement. In the complex with pyridine there is a net transfer of 0.4e to Au₄. In the complex of P(CH₃)₃ a full electron is transferred to Au₄, but in a non-uniform way: the Au atom bound to P loses 0.86e, while the three remaining Au atoms gain 1.94e among them.

To examine the performance of conventional DFT approaches in the treatment of the gold cluster-ligand interactions we compare in Table V the interaction energies of Au₄+SCN⁻ computed with functionals from increasing rungs: LDA, pure GGA, hybrid GGA,⁸⁶ meta GGA, and a standard range-separated ωB97X. The result of ωB97X closely matches the value from our optimized LRC-ωPBEα functionals in Table IV. This is the result of the gap between the donor HOMO and acceptor LUMO (1.7 eV) closely matching that of our optimized functionals. M06HF yields a weaker interaction which can be related to the increased LUMO(A)-HOMO(D) gap to 2.34 eV. By contrast, both PBE0 and BP86 which yield stronger binding, also provide negative LUMO(A)-HOMO(D) gaps, i.e., an incorrect balance of donor-acceptor properties of the constituents. The same is true, to the extreme, of the LDA result. Our conclusion is that even though BP86 and PBE0 results appear to agree with CCSD(T) (BP86: −51.8 kcal/mol; PBE0: −49.8 kcal/mol; CCSD(T): −52.9 kcal/mol), this apparent agreement is illusive as the underlying chemistry is wrong. Moreover, as seen in Fig. 2, both PBE0 and BP86 functionals strongly depart from a piecewise linearity of DFT energy with fractional electrons.

The recognition of the importance of frontier orbital gaps of donor and acceptor to the ability of DFT to describe correctly the ground state donor-acceptor binding and charge transfer is due to the seminal works on the TTF-TCNQ dimer.^17,87 The authors argued for the removal of many-electron SIE by increasing the percent HF in an xc functional. They observed that these frontier orbital gaps linearly correlate with %HF, thus allowing them to recommend the set percentage that prevents excessive ground-state charge transfer. By contrast, in our optimized functionals the HOMO(D)-LUMO(A) gaps are almost independent of the %HF. The restoration of the straight-line behavior of functionals with respect to fractional electrons due to such an optimization (see Fig. 2) prevents the artificial charge transfer and affords the DFT binding energies in good agreement with the wave function methods (see below) for the right reasons.

The energetics from the optimized LRC-ωPBEα functionals can be compared to the CCSD(T) results obtained in the same basis set (see Table IV). Before performing this comparison, some accuracy issues concerning the CCSD(T) results for gold-ligand interactions should be clarified. The wave function methods, such as CCSD(T), are much more strongly dependent on the basis set than DFT methods are. For example, our tests for Au₄-pyridine show that the difference in the binding energy obtained in vtz and avtz amounts to 4 kcal/mol in CCSD(T) but only 1 kcal/mol in the DFT method employed here. Therefore, obtaining a basis set saturated result in CCSD(T) would require a much larger basis set than DFT. A second issue concerns whether or not to correlate the outer-core 5s electrons of gold.⁸⁸ We tested the effect of correlating 5s electrons on the interaction energy on the Au₄-pyridine complex (since it contains the shortest intersystem bond). The effect again proved to be strongly basis set dependent amounting to ∼1 kcal/mol stabilization in vtz and 0.1 destabilization in avtz.

It should be stressed that MP2 cannot be reliably applied for these interactions (for reasons see Ref. 76) while SCS-MP2 can. For example, in Au₄-P(CH₃)₃ the MP2 approach yields 14 kcal/mol stronger binding than CCSD(T) of Table IV. On the other hand, SCS-MP2 overbinds only by 1.5–4 kcal/mol (see Table S1 in supplementary material⁷⁹).

The CCSD(T) binding energies are between 3 and 5 kcal/mol stronger than our DFT values. We consider this a very reasonable result because the DFT interaction energy, while accounting for short-range intermolecular correlation, which provides the bulk of the bonding (see Ref. 50), neglects the non-local correlation contribution. The D3 dispersion results are shown in the penultimate column of Table IV. The dispersion term contributes a moderate 5%–15% to the interaction energy. The sum of DFT and dispersion contributions falls within a close range of CCSD(T).

Our result for benzenethiolate anion interacting with Au₄ indicates a strong bonding (DFT+disp in excellent agreement with CCSD(T)) of covalent nature because of vanishingly small HOMO(D)-LUMO(A) gap. Another way to achieve an analogous covalency would be by the interaction of benzenethiyl radical with odd-number Au_n clusters. Indeed, such a possibility was investigated by Pasteka et al.⁵⁰ between CH₃S• and a single Au atom on the singlet and triplet surfaces. Their all-electron scalar-relativistic Douglas-Kroll (DK) result for the singlet state bonding is −56.8 kcal/mol. A qualitative agreement with our result (anion with even-numbered cluster) is noteworthy. Among other conventional DFT studies^49,38,2 there is a consensus that the Au-S bond is stronger than an Au-Au bond. Kruger et al.⁸¹ demonstrated that CH₃S• can pull a monoatomic gold wire from a gold surface. Häkkinen group showed that SR groups tend to etch gold atoms from Au_n clusters to form RS-Au-SR “staple” motives that protect the gold core. These results were obtained with PBE functional and 11-electron pseudopotential.

Our result for benzenethiol interacting with Au₄ indicates moderately strong bonding of about 28 kcal/mol (23.4–23.7 kcal/mol from DFT and 4.2 kcal/mol from dispersion). This is the interaction with the most significant dispersion contribution of 15%. By comparison, in a related system, recent work by Dufour et al.⁴⁶ using CAM-B3LYP⁸⁹ finds the binding of methanethiol with Au₁₁ to be weaker (6 kcal/mol) but of the same nature, i.e., a DA complex with the thiol acting as a donor. These authors also remark on the very strong sensitivity of the results to the choice of functional prior to choosing CAM-B3LYP.

Trimethylphosphine binds strongly to Au₄ at our DFT level by 48.3 kcal/mol plus 5.3 kcal/mol for dispersion. The literature value for the related systems includes Au-P(CH₃)₃ of Pasteka et al.⁵⁰ with binding of 22.6 kcal/mol (all electron DK). Goel et al. studied Au₄ bound to two and four P(CH₃)₃ ligands.⁴⁸ Their reported TPSS binding energies per bond are ∼43 kcal/mol in the former and 29 kcal/mol in the latter. The spread of binding energies in different functionals (as read from their graph) appears to be about 8 kcal/mol with PBE being the strongest and B3LYP being the weakest. They also find LDA to be an outlier overbinding by 15 kcal/mol compared to the rest.

The results presented in this paper argue that GKS treatment is essential to the description of donor-acceptor bonds, in general, and gold-ligand interactions, in particular. Tuning the long-range corrected functionals within the GKS scheme ensures that the IP of the donor molecule and the EA of the acceptor molecule are correctly represented within the functional by the gap between the HOMO of the donor and the LUMO of the acceptor. However, one acknowledged problem is that different monomers may require different tuning parameters.⁹⁰ The two-parameter tuning for the monomer properties of the family of LRC-ωPBEα functionals help ensure the balanced treatment of the donating and accepting properties in these complexes. The optimized functionals also show nearly linear behavior of DFT energy with fractional electrons.

Hybrid functionals also belong to the GKS class.⁹¹ However, these functionals do not guarantee that the gap between the HOMO(D) and LUMO(A) will be correctly represented. In fact, the PBE0 hybrid in our calculations for Au₄-SCN⁻ yields a negative gap. Previous work advocated the variation in the constant percent of exact exchange in hybrid functionals to match this gap. We believe that the (two-parameter) optimization of range-separated functionals is a better route to achieving this goal. Specifically, we showed the examples of homo-dimers where the HOMO-LUMO gaps are largely independent of the total percent of the HF exchange arising from both the short and long ranges in the total exchange.

The functionals that align frontier eigenvalues with IP and EA for donor-acceptor pairs are ideally suited for the predictions of the electron flow between a gold cluster or surface and ligands. A property which controls the flow, Mulliken electronegativity χ = 1/2(IP+EA), is an average of IP and EA. In the case of a metal, χ corresponds to the position of the Fermi level whereas in the case of ligands, such as these studied in the present work, it is the middle of the E_HL gap.

A number of insights into the gold-ligand interactions were obtained. In a model aurophilic dimer bound only by the correlation effects, our optimized functionals yield the interaction energy in the range of −1.7 to −1.8 kcal/mol. This dimer (in this particular configuration) is purely repulsive at the HF level. Our attractive interaction, free from artifacts masquerading as interaction energy, indicates the dimer is indeed bound by intermolecular correlation effects that are partially accounted for by the DFT. Together with the remaining D3 dispersion of ∼−4.8 kcal/mol our estimate is reasonably close to the CCSD(T)/CBS well depth of 6.13 kcal/mol.

In Au₄ interacting with a number of ligands (SCN⁻, PhS⁻, PhSH, pyridine, and P(CH₃)₃) monomer-optimized LRC-ωPBEα functionals have led to very sound binding energies. In combination with moderate contribution from the dispersion term, our values are in very good agreement with CCSD(T), i.e., 1-2 kcal/mol (pyridine, PhSH, PhS⁻, and P(CH₃)₃) to up to 3 kcal/mol (SCN⁻). It should be emphasized that the dispersion effect, as estimated by D3 method, contributes between 5% and 15% to the DA bonding in these systems. Benzenethiol forms a moderately strong DA bond with Au₄ of ∼28 kcal/mol with thiol acting the electron donor. The binding of the benzenethiolate anion is much stronger (over 60 kcal/mol) exceeding the Au-Au bonding strength. Trimethylphosphine interacts much more strongly with Au₄ than pyridine and benzenethiol.

The affinity of organo-sulfur ligands for gold stems from the larger variety of bonding situations from medium strength DA bonds of neutral thiols to strong coordinate bonds exhibited by negatively charged thiolates to covalent bonds of radical thiyls. We showed that the strength of thiolate bonding to gold cluster is capable of inducing an incipient transfer of gold atom from the cluster to ligand.

Our results for these ligands are consistent with other wave function predictions⁵⁰ and may serve as a platform for testing wide range of functionals. Such tests should take the dispersion energy into consideration.

This work was supported by the National Science Foundation (NSF) (Grant No. CHE-1152474) and by the Polish Ministry of Science and Higher Education (Grant No. N204 248440). M.H. was supported by the project “Towards Advanced Functional Materials and Novel Devices: Joint UW and WUT International PhD Programme” operated within the Foundation for Polish Science MPD Programme, implemented as a part of the Innovative Economy Operational Programme (EU European Regional Development Fund). We wish to thank Professor Janusz Zachara and Professor Wojciech Grochala for reading and commenting on the paper.