Density-functional tight-binding for phosphine-stabilized nanoscale gold clusters

Abstract

We report a parameterization of the second-order density-functional tight-binding (DFTB2) method for the quantum chemical simulation of phosphine-ligated nanoscale gold clusters, metalloids, and gold surfaces. Our parameterization extends the previously released DFTB2 “auorg” parameter set by connecting it to the electronic parameter of phosphorus in the “mio” parameter set. Although this connection could technically simply be accomplished by creating only the required additional Au–P repulsive potential, we found that the Au 6p and P 3d virtual atomic orbital energy levels exert a strong influence on the overall performance of the combined parameter set. Our optimized parameters are validated against density functional theory (DFT) geometries, ligand binding and cluster isomerization energies, ligand dissociation potential energy curves, and molecular orbital energies for relevant phosphine-ligated Au_n clusters (n = 2–70), as well as selected experimental X-ray structures from the Cambridge Structural Database. In addition, we validate DFTB simulated far-IR spectra for several phosphine- and thiolate-ligated gold clusters against experimental and DFT spectra. The transferability of the parameter set is evaluated using DFT and DFTB potential energy surfaces resulting from the chemisorption of a PH₃ molecule on the gold (111) surface. To demonstrate the potential of the DFTB method for quantum chemical simulations of metalloid gold clusters that are challenging for traditional DFT calculations, we report the predicted molecular geometry, electronic structure, ligand binding energy, and IR spectrum of Au₁₀₈S₂₄(PPh₃)₁₆.

---

1 Introduction

Atomically precise ligated gold clusters of nanometer dimension receive continued attention due to their unique catalytic properties^1–3 and well-defined discrete electronic energy levels⁴ that potentially offer greater flexibility and control over the more metal-like states of the corresponding bare gold nanoparticles and complexes.^3,5 Nanoscale gold clusters have shown high catalytic activity and selectivity for certain reactions at low temperature, such as the oxidation of carbon monoxide, propene and alcohols, or the hydrogenation of acetylene.^3,6–8 Commonly chosen “capping ligands” employed to stabilize atomically precise nanoscale gold clusters are typically thiolates and phosphines, which prevent aggregation, coalescence and unlimited growth during synthesis.⁹ Post-treatment is usually required to remove some (or all) of the ligands to allow interaction with the substrate or reactants.^10–14 To achieve control over such complex catalytic systems it is vitally important to understand the relationship between molecular and electronic structure, often studied by a combination of experimental and theoretical approaches.^15,16 In addition, a better knowledge of the energetics associated with ligand removal is required to identify how post-treatment can be done without inducing concomitant side-effects such as agglomeration. Density functional theory (DFT) methods are most often employed in theoretical investigations, as they are capable to accurately describe electronic, geometrical, and vibrational structure of gold clusters and nanoparticles.^16–29 In particular, the DFT-based simulation of IR, Raman, and UV-vis spectra has been achieved for a range of gold clusters, thus providing useful theoretical fingerprints to distinguish between bonding arrangements and orientations between gold atoms and ligands, and ligand–ligand interactions within clusters.^23,30–34

Unfortunately, DFT calculations can become prohibitively expensive with system size,³⁵ and routine theoretical investigations are limited to moderate system sizes. One way of reducing the computational cost is to entirely neglect the ligands and only perform DFT calculations on the gold cores. An alternative way is to employ simplified ligand models where e.g. triphenylphosphine (PPh₃) is replaced by phosphine (PH₃) or trimethylphosphine (PMe₃).^{17,18,22,25,26,28,29,36} This, as with the complete ligand removal approach, has the additional benefit that conformational searching is simplified, as the torsions of the three phenyl groups per ligand give rise to a large number of local minima with similar energies. Nevertheless, it is well known that the electronic effects of the larger ligands are different from those of the smaller ones, for instance inductive effects,^36–39 and a computational truncation of the ligands will influence the chemistry and therefore description of the catalytic properties in calculations. Integrated schemes such as ONIOM⁴⁰ may be used to capture such electronic effects in calculations on the untruncated “real” systems; however, the choice for high and low levels of theory and the definition of the interface between them is not straightforward. A rigorous ONIOM study requires benchmarks of the selected methods against a high-level calculation for the “real” system,⁴¹ which is often computationally unfeasible. In practice therefore, integrated methods are often difficult to employ in the context of ligated nanoscale metal clusters. Besides the computational effort related to the proper modeling of the ligands, the size of the metal cluster itself can become problematic for conventional DFT studies, severely impacting the size range of gold clusters to be investigated for property control and fine-tuning. Exacerbating this problem is the recent emergence of larger, so called metalloid, clusters such as Au₁₀₈S₂₄(PPh₃)₁₆.²⁵ To make matters even worse, interactions of deposited gold clusters with substrate surfaces such as SiO₂ and TiO₂ may play an important role in the catalytic reaction,^10–14,42 which further increases the computational expense of DFT studies. A recent review on connections between theory and experiment for gold nanoclusters has thus posed the question as to how theoretical calculations can be expanded to treat larger sizes and length scales.¹⁵

Semi-empirical electronic structure methods offer the capability to simulate large systems with explicit inclusion of electronic structure by introducing empirical parameters and methodological approximations to rigorous ab initio or first principles methods.^43–48 Among them, density-functional tight-binding (DFTB),^45,49–53 an approximation to DFT formulated in the framework of non-orthogonal tight binding, has emerged as one of the most accurate, and potentially versatile choices. DFTB is capable of simulating systems containing many thousands of atoms with an accuracy comparable to traditional DFT methods.^46,54,55 As the DFTB method takes advantage of the two-center approximation, tabulated Hamiltonian and overlap integrals within the Slater–Koster scheme,^51,56 it is two to three orders of magnitudes faster than DFT. In order to apply DFTB into the theoretical study of ligated gold clusters, it is necessary to provide accurate parameters for all binary chemical element interactions, most notably Au–S and Au–P. The Au–S interaction in combination with the “mio” parameter set^51,52,57 was included in the “auorg” parameter set for gold–thiolates clusters,^58,59 but the parameters for the Au–P interaction have not been developed.

In this work, we report a parameterization of the Au–P interactions for the second-order DFTB (DFTB2) method for the sake of compatibility with the previously developed parameters. The accuracy of these new DFTB parameters is probed by comparing the corresponding properties against DFT-calculated values and available experimental results in terms of: (1) root mean square deviations of optimized geometries, (2) energetic properties, i.e. ligand binding energies, relative isomer energies, and ligand dissociation energy profiles, (3) electronic structure, (4) vibrational normal modes of Au–P containing clusters, and (5) the adsorption of a PH₃ molecule on the gold (111) surface. Finally, we present DFTB-based predictions for structural, energetic, and vibrational properties of the recent experimentally reported metalloid gold cluster Au₁₀₈S₂₄(PPh₃)₁₆.²⁵

2 Methodology and computational details

2.1 Brief overview of DFTB2

A comprehensive review of DFTB methods can be found elsewhere⁶⁰ and will not be repeated here. In this work, we only focus on generating parameters for the DFTB2 method, which is also referred to as self-consistent-charge (SCC)-DFTB.⁵¹ The DFTB2 total energy can be viewed as a 2nd order Taylor expansion of the Kohn–Sham energy with respect to a reference initial electron density ρ₀ and electron density fluctuations Δρwhere is the initial Hamiltonian constructed from the superposition of neutral atomic densities in a two-center approximation,⁴⁹ and the |Ψ_i〉 are occupied valence molecular orbitals (MOs), expanded as a linear combination of optimized pseudo-atomic orbitals |ϕ_μ〉. The optimization of these pseudo-atomic valence orbitals and orbital densities for a given chemical element constitutes the determination of the electronic parameters. Δq_A is a point charge⁶¹ on atom A, and γ_ABΔq_AΔq_B represents the Coulomb interaction energy between the two point charges;⁵¹ when A = B, γ_AA is the chemical hardness or second derivative of the total energy with respect to the charge on atom A. γ_AB is defined aswhere S is an exponentially decaying short-range function that depends on the distance between the two atoms r_AB = |r_A − r_B| and their chemical hardness, given in form of the so-called Hubbard parameter U. The latter is calculated prior to molecular DFTB calculations for each chemical element using the DFT method, typically employing the Perdew–Burke–Ernzerhof (PBE) functional.⁶²

In the framework of the two-center approximation, the Hamiltonian integrals and the overlap integrals 〈ϕ_μ|ϕ_ν〉 are pre-tabulated for each chemical element pair. As in all tight binding approaches, the DFTB2 total energy consists of an electronic energy, which is the sum of the first two terms in eqn (1), and a summation over all unique repulsive potentials between two atoms E^rep_AB. The latter are formulated as a two-center term that depends only on the chemical element type of atoms A and B and their interatomic distance r_AB.⁶⁰ In principle, these pairwise repulsive potentials can be pre-calculated analytically from DFT for diatomic molecules. However, it was found that the performance of such DFT-based repulsive potentials is usually not sufficiently accurate for general purposes.⁵⁹ Therefore, in practice, one computes DFT- or wave function theory (WFT)-based reference relative energies for model systems that contain various bond lengths of the chemical element pair in question, and fits the total repulsive energy such as to minimize the difference between reference and resulting DFTB relative energies. For this purpose, the repulsive potentials are approximated by a combination of exponential and spline functions

where “r⁰_AB,n” is a spline knot at the n^th interval, and the “a_AB,n” are the polynomial spline coefficients. These variables are considered the free empirical parameters and optimized either by hand or automatically^63–65 to minimize the difference between reference and DFTB relative energy for a series of training systems. In addition to relative energies, repulsive potentials can also be optimized to fit molecular geometries by minimizing the difference between the reference DFT or WFT energy gradient and the DFTB energy gradient for a set of equilibrium and non-equilibrium geometries in the training set. The possible constraints on the repulsive potentials include a cutoff radius, a limit on the number of allowable extrema, and a continuity requirement up to a given derivative order.

2.2 Gold–phosphorus parameterization

As indicated above, there are two groups of DFTB2 parameters that need to be determined: (1) the electronic parameters, and (2) the repulsive potentials for pairs of chemical elements. The electronic parameters are comprised of the radii r^wf used in the definition of atomic confinement potentials for generating pseudo-atomic orbitals |ϕ_μ〉 as well as the confinement radii r^dens for the atomic density; additional electronic parameters are the atomic orbital energies and the atomic Hubbard parameters U for each chemical elements.⁵¹ While the AO energies and the Hubbard parameters U are normally taken from DFT calculations of the free atom, the other electronic parameters and pairwise repulsive potentials are subject to optimization with the goal to reproduce certain desired properties; for instance, electronic band structure, atomization energies, reaction energies, and geometries (energy gradients or atomic forces). In order to parameterize the Au and P interactions for DFTB2, we adopted the Au electronic parameters from the “auorg” set published by Fihey et al. (referred to as auorg^α),⁵⁸ and a modified version of “auorg” by Oliveira et al. (referred to as auorg^χ).⁵⁹ The difference between these two parameter sets lies in the Au 6p-orbital energy; in the auorg^α set it was taken as the true PBE orbital energy, while in the auorg^χ set it was empirically shifted upward by ≈+0.0279 hartree. The main purpose for this orbital energy shift was to obtain improved values for cohesive energies of pure gold nanoclusters with respect to PBE.⁵⁹ Following the work of “auorg”, only 5d and 6 s valence electrons are considered, in total 11 valence electrons per Au atom. The parameters of the other elements in the auorg^α and auorg^χ sets were taken from the “mio” parameter set.^51,52 Consequently, we adopted the electronic parameters for P taken from the “mio” parameter set as well. It is important to mention that the 3d orbital energy of the P atom had been shifted by +0.5 hartree from its PBE-calculated value of 0.02044 hartree, according to Gaus et al. to reduce the overbinding in phosphate compounds, as the shift had also been adopted in the “3ob” parameter set for the DFTB3 method.^52,66 However, in the case of Au–P interactions, such a drastic shift of the P 3d virtual orbital energy introduces significant underbinding. We therefore decided to investigate the effect of the P 3d orbital energy level in detail. We systematically increased the value of the P 3d orbital energy by increments of 0.1 hartree from its PBE computed and the “mio” shifted value, and used our genetic algorithm (GA) optimization tool⁶⁵ to automatically generate repulsive potentials for the Au–P interaction with these orbital energy shifts. In this way we tested the performance of the resulting Au–P parameterization with special consideration of geometries and binding energies for the adsorption of a PH₃ on Au (111). An energy value of 0.12044 hartree was determined as the optimal compromise for the P 3d orbital. We refer to Table S1 in the ESI† for the performance of DFTB2/mio with different P 3d orbital energies for a selected test set of chemical reactions involving H, C, N, O, P, and S containing compounds.

We decided to introduce a nomenclature for denoting the different choices of Au 6p and P 3d orbital energies in our parameters. Following Oliveira et al.⁵⁹ we denote the original PBE Au 6p energy value by ‘α’ and its modified value by ‘χ’. The shifted “mio” P 3d energy will leave these notations unchanged, while the use of our new optimized P 3d orbital energy value of 0.12044 hartree will be denoted with a prime, ‘′’. Hence, auorg^α coupled with the original “mio” P is unchanged, while the auorg^α and auorg^χ coupled with the new P 3d-orbital energy are referred to as auorg^α′ and auorg^χ′, respectively. Because auorg^χ exhibits the largest underestimation of the Au–P electronic binding energy (both Au 6p and P 3d virtual orbital energies are shifted upwards from PBE values), we did not generate the Au–P repulsive potential for this orbital energy combination. Table 1 summarizes the nomenclature of the three different parameter sets and the relationship with their employed Au and P virtual atomic orbital energies.

Table 1 Energies in hartree for Au 6p and P 3d virtual atomic orbitals in different parameter sets

	auorg^α	auorg^α′	auorg^χ′
ε ^6pAu	−0.02786	−0.02786	−0.00001
ε ^3dP	0.52044	0.12044	0.12044

---

The Au–P repulsive potentials were optimized on the basis of shell-resolved self-consistent charge electronic energies to fit ligand binding energies and forces for the training set listed in Table 2 using our in-house genetic algorithm (GA)-based parameterization tool.⁶⁵ After several preliminary tests, we employed a cutoff radius of 4 Å and 5 spline knots (the maximum for n in eqn (3) was equal to 5), one allowable extremum, and a continuity requirement up to the third derivative. For the GA optimization, population sizes of 3000 and 5000 generations were employed with two-point crossover and random mutation rates of 0.9 and 0.2, respectively. The GA was used to minimize a scoring function F^score defined as the fitness of the parameter sets with respect to DFTB values and reference data according to the formula

where ΔE^bind = E^bind_eDFTB − E^bind_ref is the deviation in ligand binding energies, ΔF^force = F^force_eDFTB − F^force_ref is the deviation in forces, W^bind_i and W^force_i are weight factors of the i^th binding energy and the force of the i^th structure, respectively, N_eq is the number of fitting data points, N_i is the number of atoms in the i^th compound, and eDFTB stands for the DFTB energy without the repulsive potential term. In the training, we empirically shifted the reference ligand binding energies of PH₃ by −24 kcal mol⁻¹ for auorg^α′ and −22 kcal mol⁻¹ for the auorg^χ′ to better reproduce ligand binding energies of larger ligands and the adsorption of PH₃ on the Au (111) surface. Surprisingly, auorg^α′ and auorg^χ′ optimized Au–P repulsive potentials are almost identical, see Fig. S1 in the ESI.† In principle, the Au–P repulsive potential for auorg^α can also be optimized in a similar way to maximize its performance, however, for the sake of simplicity and transferability we decided to use the same Au–P repulsive potential for auorg^α′ and auorg^α.

Table 2 Training set for Au–P repulsive potential fitting and weights for each complex used in the fitting with eqn (4)

Forces (equilibrium geometries)
Complexes	Structures	Weights
Au1–PH2		1.0
[Au1–PH3]^+		1.0
[Au3–PH3]^+		1.0
[Au4–PH2]^+		1.0
[Au4–PH3]^2+		1.0
[Au6–PH3]^2+		1.0
[Au6(planar)–PH3]^2+		1.0

---

Forces (distorted geometries)
Complexes	ΔRAu–P (Å)	Weights
[Au1–PH3]^+	−0.3	0.5
[Au1–PH3]^+	−0.2	0.5
[Au1–PH3]^+	−0.1	0.5
[Au3–PH3]^+	+0.1	0.5
[Au3–PH3]^+	+0.2	0.3
[Au3–PH3]^+	+0.3	0.1

---

Ligand binding energies
Reactions	ΔE (kcal mol^−1)	Weights
[Au3]^+ + PH3 ⇒ [Au3–PH3]^+	−55.53	1.0
[Au4]^2+ + PH3 ⇒ [Au4–PH3]^2+	−76.27	1.0
[Au6]^2+ + PH3 ⇒ [Au6–PH3]^2+	−65.01	1.0
[Au6(planar)]^2+ + PH3 ⇒ [Au6(planar)–PH3]^2+	−64.58	1.0

---

2.3 Computational details

In the original parameterization of the auorg set,⁵⁸ the generalized gradient approximation (GGA) PBE density functional was selected to generate reference geometries and cluster binding energies. In this work, we opted for the TPSS density functional⁶⁷ because it was noted by Kepp⁶⁸ and Goel et al.⁶⁹ that this meta-GGA functional reproduces experimental or high-level theory bond energies and ligand–gold distances better than the standard GGA PBE functional. For the training set of small clusters in Table 2, the reference data were computed by the TPSS in combination with Ahlrich's triple-zeta valence polarized basis set (def2-TZVP).⁷⁰ No dispersion correction was employed in the calculations generating the training set in order to avoid complications originating from a possible convolution of DFTB repulsive energy terms and the long-distance dispersion term.

In order to benchmark the accuracy of the new parameters, ligand binding energies and optimized geometries were compared to their TPSS counterparts for various complexes of PH₃, PMe₃, PPh₃ and small- to moderate-sized gold clusters. For the larger complexes, which are listed in Table 3, the experimental structures from the Cambridge Structural Database (CSD) were used as the reference geometries. The different phosphorus-containing ligands considered in this test set are triphenylphosphine (PPh₃), 1,1-bis(diphenylphosphino) propane (dppm), 1,3-bis(diphenylphosphino) propane (dppp), methyldiphenylphosphine (PMePh₂), tris(2-(diphenylphosphino) ethyl) phosphine (PP₃) and 1,8-bis(diphenylphosphino) octane (dppo). Additionally, there are thiol-containing ligands in some clusters comprised of reduced S₂⁻ and meta-methylbenzenethiol (m-MBT). For these test sets, the reference calculations were performed at the TPSS/def2-SVP⁷⁰ level of theory. Here we included the empirical D3 dispersion contribution^71,72 in both DFT and DFTB2 calculations.⁷³ The D3 dispersion correction was used to improve description of ligand–ligand interactions and ligand effects which are deemed to be important factors on the structural, electronic, and vibrational properties of ligated gold clusters.²³ In the DFT calculations we employed an effective core potential (ECP)⁷⁴ for the Au atoms. In order to accelerate the DFT calculations, we employed the resolution-of-the-identity (RI) approximation with the corresponding auxiliary basis sets.⁷⁵ All non-periodic DFT calculations were carried out using the implementation of the ORCA code.⁷⁶ DFTB2 single point energy and geometry optimization calculations were performed with the DFTB+ code.⁷⁷ All DFTB calculations were carried out with the shell-resolved SCC option “OrbitallyResolvedSCC = Yes”. The auorg parameter was designed such that this option mostly affects the charge distribution on the gold atoms themselves and their interactions with the other elements.

Table 3 Experimental crystal structures taken from the CSD database for the test set

Complexes	CSD codes
[Au6(dppp)4]^2+	BOTSOS^78
[Au6(PPh3)6]^2+	CATPAO10 (ref. 79)
[Au7(PPh3)7]^+	BIXZAK^80
[Au8(PPh3)7]^2+	BASWUN10 (ref. 80)
[Au8(PPh3)8]^2+	OPAUPF^81
[Au8S2(dppm)4]^2+	LEVKIJ^27
[Au9(PPh3)8]^3+ (D2h)	MIVPOX^82
[Au11(PMePh2)10]^3+ (C3v)	ZUCMAL^83
[Au11(PMePh2)10]^3+ (D4d)	ZUCMEP^83
[Au13(dppm)6]^5+	LEVKAB^27
[Au20(PP3)4]^4+	POFPUX^84
Au22(dppo)6	TOCFIC^85
[Au38(m-MBT)20(PPh3)4]^2+	CEMZIG^86
Au70S20(PPh3)12	TELMUV^26

---

To test the accuracy of the DFTB parameters for the prediction of vibrational spectroscopic data, we compared DFTB- and PBE-calculated far-IR spectra to the experimental spectra for [Au₆(dppp)₄]²⁺, [Au₈(PPh₃)₈]²⁺ and [Au₉(PPh₃)₈]³⁺ clusters. The DFTB2 IR vibrational spectroscopy calculations were computed using the GAMESS-US code.^87,88 The simulated IR spectra were obtained by convoluting the calculated stick spectra using a Lorentzian line shape function with 3 cm⁻¹ full width at half maximum.

The transferability of the DFTB parameters was evaluated by comparing DFT and DFTB energy landscapes of PH₃ adsorption on the Au(111) surface. A 4 × 4 supercell consisting of 4 layers of Au atoms was cut from bulk and a vacuum layer of 20 Å was added to effectively suppress through-space slab–slab interactions. For the adsorption energy scans, all Au atoms were fixed. Selective geometry relaxation was used by constraining the P atom of the PH₃ molecule in all but the z-direction, while the H atoms the were fully optimized. For the DFT calculations, the PBE functional was used with the projector augmented wave (PAW) approach,⁸⁹ the kinetic energy cutoff was 450 eV, and k-point grids were generated dynamically using a 3 × 3 × 1 Monkhorst–Pack scheme. The DFT calculations were carried out with the Vienna Ab initio simulation package (VASP) program in conjunction with the provided “PAW_PBE” pseudopotentials.^90,91 The convergence criteria were set to 10⁻⁶ eV for achieving self-consistent field energies, and 0.005 eV Å⁻¹ for the maximum force in case of geometry optimizations.

For the demonstration application, a case study on the large Au₁₀₈S₂₄(PPh₃)₁₆ metalloid was performed. The initial positions of Au, S, and P atoms were taken from the experimental crystal structure (CSD code: DAFLOO). The PH₃ and PPh₃ ligands were added based on the position of the P atoms in the CIF file. Then, these complexes went through a two-step geometry optimization: first, only C and H atoms were optimized while Au, S and P atoms were fixed, and second, all atoms were optimized. The fully relaxed complexes were used to calculate binding energy, electronic properties and the IR spectrum.

3 Results and discussion

3.1 Performance for small- and moderate-sized clusters

Small-sized clusters. The differences between DFTB and TPSS/def2-SVP levels of theory for optimized geometries and ligand binding energies were evaluated for nine model complexes, which were constructed from three small-sized gold clusters Au_n (n = 2, 3, and 4) with three phosphorus-based donor ligands, including PH₃, PMe₃ and PPh₃. Overall, Fig. 1 shows that all three parameter sets are performing very well, as measured by the root mean square deviation (RMSD) over Au and P atom positions, with all values below 0.13 Å. auorg^α is superior among the three sets with a maximum RMSD of only 0.09 Å, while auorg^α′ produces typically the largest error. The averaged and normalized DFTB ligand binding energy generally shows overbinding for these small-sized complexes relative to the reference TPSS/def2-TZVP level of theory, strongest for auorg^α′ and least for auorg^α (see Fig. 1 and Table S2 in the ESI†). We note that the deviation in ligand binding energies decreases with the size of gold clusters as well as the size of the ligands: errors are largest for the PH₃ ligand and smallest for the PPh₃ ligand. The variation of the binding energy deviation can be attributed to the electronic interaction rather than the repulsive potential: while all Au–P bond lengths in these complexes are almost the same, ≈2.35 Å, the averaged and normalized ligand binding energies increase with the Au cluster size as well as ligand size. For instance, the ligand binding energies of PH₃, PMe₃, PPh₃ with Au₄ cluster are −33.91, −46.44, and −48.69 kcal mol⁻¹ respectively. Since the overbinding can be caused by low virtual orbital energies, it makes sense that the auorg^α′ is mostly affected by it. As for the ranking in performance for averaged and normalized ligand binding energies, we find the same order as for the geometries: auorg^α is the best parameter in terms of binding energy, while auorg^α′ performs the worst.

Fig. 1 RMSD over atomic positions (upper panel), and deviation in averaged and normalized ligand binding energies (lower panel) for the small-sized gold clusters with different phosphine ligands (L = PH₃, PMe₃, and PPh₃). The RMSD over atomic positions only considers Au and P atoms.

Moderate-sized clusters. Fig. 2 displays the same performance data for moderate-sized clusters Au_n (n = 6–22) with phosphine or trimethylphosphine ligands. The DFTB optimized geometries agree acceptably well with those from TPSS/def2-SVP, with RMSDs less than 0.5 Å in most cases. Exceptionally, three cases of [Au₇(PH₃)₇]⁺, Au₂₂(PH₃)₁₂ and [Au₂₀(PMe₃)₁₆]⁴⁺ complexes exhibit RMSD values larger than 0.7 Å. While the RMSD of auorg^α′ and auorg^χ′ for these three clusters, and auorg^α for Au₂₂(PH₃)₁₂ [Au₂₀(PMe₃)₁₆]⁴⁺ remain in an acceptable range around ≈0.7 Å, the auorg^α RMSD value of 1.18 Å appears problematic for [Au₇(PH₃)₇]⁺.

Fig. 2 RMSD over atomic positions (upper panel), and deviation in averaged and normalized ligand binding energies (lower panel) for the moderate-sized phosphine-stabilized gold clusters with phosphine ligands (L = PH₃ and PMe₃). The RMSD over atom positions only considers Au and P atoms.

After closer inspection, we conclude that this large value results from a strongly distorted Au₇ core shape that can only be stabilized when larger ligands such as PMe₃ or PPh₃ are used. It is worth mentioning that these model systems are experimentally not stable in general, because the gold core prefers a planar structure rather than a 3D geometry, unless ligands are present.^92,93 The model geometries here were constructed by replacing the experimentally used ligands (listed in Table 3) with PH₃ or PMe₃ ligands. For these hypothetical models, large RMSDs are expected because the presence of many local minima on the potential energy surfaces complicates the geometry optimization, and can cause the DFT and DFTB geometry optimization to converge to different local minima. For the [Au₇(PH₃)₇]⁺ cluster, lower RMSDs were observed for auorg^α′ and auorg^χ′ because these two parameters increase the contribution of the P 3d-orbital in stabilizing the complex due to its lower orbital energy.

Compared to the situation in small-sized clusters, ligand overbinding is less prominent in the medium-sized systems for both auorg^α′ and auorg^χ′ parameters, with energy deviations as high as −18 kcal mol⁻¹ for the smallest PH₃-ligated clusters (see Fig. 2 and Table S3 in the ESI†). Again we find generally that the averaged and normalized ligand binding energy deviations become smaller when the larger PMe₃ ligands are used. Again, auorg^α has smaller overbinding than the other two parameters for PH₃ ligands, and this tendency turns to underbinding when the larger PMe₃ ligands are used. The performance ranking for the three parameter sets is less clear-cut as in the case of the small-sized clusters, but especially in terms of ligand binding energies we find the tendency confirmed that auorg^α tends towards least overbinding and auorg^α′ towards greatest overbinding, with auorg^χ′ somewhere in between but closer to auorg^α′. Since small ligands themselves consistently tend to increase overbinding, we find that auorg^α performs best for small ligands and small Au clusters, while auorg^α′ or auorg^χ′ parameters are better for larger ligands and larger gold clusters.

3.2 Performance for large-sized clusters

The accuracy of the new DFTB2 parameters for geometries and ligand binding energies was further assessed for a series of gold clusters (Au_n, n = 6–70) in complexes with their experimentally used larger ligands including PPh3, PMePh₂, PP₃, dppp, dppm, dppo, and m-MBT, see Table 3. In this test, the experimental crystal structures were used as the reference geometry for RMSD evaluation, rather than DFT geometries. In both DFT and DFTB calculations, all counterions were removed and no symmetry constrains were applied. In this section, only the performance of auorg^α and auorg^α′ parameters is presented and discussed in the main text. The performance of the auorg^χ′ parameter is presented in the ESI.†

Geometry. The RMSDs between experimental and computed cluster geometries shown in Fig. 3 demonstrate that the DFTB2 methods are able to reproduce X-ray structures, with maximum RMSDs equal to 0.46 Å and 0.37 Å for auorg^α and auorg^α′, respectively. The same RMSDs of auorg^χ′ are shown in Fig. S2 in the ESI.† In fact, the DFTB2 geometries even outperform TPSS/def2-SVP geometries for many clusters in the test set, as for instance in the case of [Au₆(dppp)₄]²⁺ (BOTSOS), [Au₆(PPh₃)₆]²⁺ (CATPAO10), and [Au₈S₂(dppm)₄]²⁺ (LEVKIJ). In the previous section, [Au₇(PH₃)₇]⁺, Au₂₂(PH₃)₁₂ and [Au₂₀(PMe₃)₁₆]⁴⁺ model clusters were found to be the most problematic cases. Here, with the experimentally used ligands, DFTB-optimized structures of [Au₇(PPh₃)₇]⁺, [Au₂₀(PP₃)₄]⁴⁺, and Au₂₂(dppo)₆ have lower RMSD values than their previously described truncated model clusters with PH₃ and PMe₃. auorg^α′ slightly outperforms auorg^α in most of the cases (by only up to 0.1 Å for the [Au₂₂(dppo)₆] neutral complex). It is important to keep in mind that the RMSDs over atom position convolutes deviations in angles and torsions with the bond lengths. If one only considers bond length comparison, Au–Au and Au–P deviations are within 0.06 Å from experimental values.

Fig. 3 RMSD over atomic positions for the large-sized phosphine-stabilized gold clusters. The RMSD of atomic positions considers Au, and P atoms for all large-sized phosphine-based gold clusters, [Au₁₁(PMePh₂)₁₀]³⁺_# denotes [Au₁₁(PMePh₂)₁₀]³⁺ (C_3v), [Au₁₁(PMePh₂)₁₀]³⁺_* denotes [Au₁₁(PMePh₂)₁₀]³⁺ (D_4d), [Au₃₈(L)₂₀(PPh₃)₄]²⁺ denotes [Au₃₈(m-MBT)₂₀(PPh₃)₄]²⁺.

In order to graphically illustrate the differences between experimental and computed DFT and auorg^χ′ geometries, Fig. 4 shows X-ray and optimized structures for the following ligated Au clusters: [Au₆(dppp)₄]²⁺ (BOTSOS), [Au₇(PPh₃)₇]⁺ (BIXZAK), [Au₈(PPh₃)₈]²⁺ (OPAUPF), and [Au₉(PPh₃)₈]³⁺ (MIVPOX-D_2h). The analogous comparison for auorg^α and auorg^χ′ is shown in Fig. S3 in the ESI.† The overlapped structures were determined by a minimization procedure, which includes recentering and rotation to minimize the RMSD using the quaternion algorithm.⁹⁴ The figure highlights the good performance of the theoretical methods in the description of the cluster core geometries relative to experiment. We note that ligand orientations can be strongly impacted by crystal field effects that are not present in our gas phase theoretical calculations, and thus we will not discuss differences in ligand geometries.

Fig. 4 Overlap of experimental crystal structure (Au in gold, P in orange and C in grey) and optimized auorg^α′ and DFT structures. auorg^α′ and DFT structures are represented by light red and sky blue, respectively. The gold nanoclusters considered in this figure are (A) [Au₆(dppp)₄]²⁺ (BOTSOS), (B) [Au₇(PPh₃)₇]⁺ (BIXZAK), (C) [Au₈(PPh₃)₈]²⁺ (OPAUPF), and (D) [Au₉(PPh₃)₈]³⁺ (MIVPOX-D_2h).

Binding energy. The evaluation of predicted averaged and normalized ligand binding energies follows the schemes used above for the smaller-sized and moderate-sized clusters. The ligand binding energies are listed in Table S4 in the ESI,† the relative deviations of the DFTB methods from TPSS/def2-SVP data are shown in Fig. 5 for auorg^α and auorg^α′, and Fig. S4 in the ESI† for auorg^χ′.

Fig. 5 Deviation in averaged and normalized ligand binding energies for the large-sized phosphine-stabilized gold clusters in reference to the TPSS/def2-SVP binding energies, [Au₁₁(PMePh₂)₁₀]³⁺_# denotes [Au₁₁(PMePh₂)₁₀]³⁺ (C_3v), [Au₁₁(PMePh₂)₁₀]³⁺_* denotes [Au₁₁(PMePh₂)₁₀]³⁺ (D_4d), [Au₃₈(L)₂₀(PPh₃)₄]²⁺ denotes [Au₃₈(m-MBT)₂₀(PPh₃)₄]²⁺.

All predictions by DFTB tend towards underbinding, with only three minor exceptions. From the smallest [Au₆(dppp)₄]²⁺ up to [Au₁₃(dppm)₆]⁵⁺, the auorg^α′ binding energy deviations of ≤5 kcal mol⁻¹ are in very good agreement with the DFT reference. auorg^α trends towards strong underbinding in all cases. In the case of clusters with more core Au atoms (n ≥ 20), auorg^α and auorg^α′ have both noticeable underbinding as high as 20 and 15 kcal mol⁻¹, respectively. The results are consistent with the overall trend that was already seen in the small-sized and moderate-sized ligand complexes above. To further investigate the effect of geometry on the ligand binding energies, TPSS single point energy calculation using the DFTB optimized geometries were carried out, and the corresponding deviations in ligand binding energies are shown in Fig. 5 and S4 in the ESI,† indicated as usual by the at-the-geometry-of symbol “//” symbol. It becomes immediately obvious that the deviation of ligand binding energies with TPSS single point energy refinement is reduced in all cases, with a maximum absolute deviation of ≤9 kcal mol⁻¹. It follows that, if highly accurate ligand binding energies are required, DFTB geometry optimizations followed by TPSS single point energy calculations can provide a reasonable “shortcut” over straightforward DFT calculations. In addition to the binding energies, DFTB isomerization energies are compared to the DFT values as well for [Au₉(PPh₃)₈]³⁺ and [Au₁₁(PMePh₂)₁₀]³⁺ clusters, see Fig. S5 in the ESI† for more details.

3.3 Gold cluster–ligand bond dissociation energy curves for the [Au₈(PPh₃)₈]²⁺ complex

Ligand removal of nanoscale gold clusters is a key step in making the clusters catalytically more active by increasing the gold core interaction with the substrate or reactants. Previous studies have reported the partial removal of ligands of atomically-precise gold clusters Au_n(PPh₃)_m (n = 8, 9, 11, and 101) on titania after undergoing different treatments such as calcination and acid-washing.^10,11 One of the primary applications of the DFTB methodology developed here is the theoretical study of cluster fragmentation and the catalytic reaction mechanisms of clusters with dissociated ligands. For validation against DFT, we therefore compare rigid energy scans for Au–P bond dissociation between TPSS/def-SVP and DFTB2 methods for the aforementioned [Au₈(PPh₃)₈]²⁺ cluster as a representative for the experimentally relevant complexes with larger ligands. Starting from the DFT optimized geometries of the cluster, the Au–P bond of the ligand in question was simply stretched up to a distance of 20 Å. The relative energy ΔE of this practically dissociated geometry was defined in all methods as 0 kcal mol⁻¹. In reality, complex structural relaxation of the cluster and ligand would occur obscuring methodological differences, which is the reason for presenting rigid scans. Fig. 6 shows four different bond dissociation energy curves that correspond to four ligand detachment scans.

Fig. 6 Gold cluster–ligand rigid bond dissociation energy curves of [Au₈(PPh₃)₈]²⁺.

Overall, the DFTB curves mimic the TPSS curves closely in both energy of the binding region and the shape of the energy curves. The rigid scan does not include a barrier and converges within about 12 Å to the dissociation limit. As for the DFTB curves, we observe underbinding in the region from 2.5 to 6 Å in all of the plots (cf.Fig. 5). DFTB is strongly underbinding in the case of one ligand by more than 30 kcal mol⁻¹ (upper right curve in Fig. 6). Here, the Au atom involved in the Au–P bond has a more “surface”-like binding with the other Au atoms as opposed to a pyramidal shape as in the other ligand cases. For more “cluster”-like detached Au atoms, the deviations are in the range of ≈5–10 kcal mol⁻¹. In all cases, DFTB2/auorg^α′ and DFTB2/auorg^χ′ are very similar to each other and outperform DFTB2/auorg^α. As is typical for the DFTB method in general, the entrance region of the binding well, here around 6 Å, is underbinding since the atomic orbitals are compressed with fixed electronic confinement radii optimized to describe geometries and binding energies.⁹⁵ The underestimation of phosphine ligands binding to gold “surface” will be discussed further below regarding the adsorption of a PH₃ molecule on the gold (111) surface.

3.4 Electronic structures of phosphine-stabilized gold clusters

To evaluate the ability of DFTB to accurately describe electronic structures of phosphine-stabilized gold clusters, we investigated the frontier orbital energies (HOMO-1, HOMO, LUMO, LUMO+1) and HOMO–LUMO gaps (HLGs) of a number of medium-sized gold clusters using DFT and DFTB methods. In this discussion, we concentrate on the auorg^α′ parameter, but show the corresponding performance of the other parameters in the ESI.†Fig. 7 shows the TPSS/def2-SVP and DFTB/auorg^α′ calculated energy levels of the frontier molecular orbitals for selected gold clusters. The plots showing frontier orbitals computed with auorg^α′ and auorg^χ′ are shown in Fig. S5 in the ESI.† The total charge in the selected clusters changes from +1 to +3, with the orbitals of the least charged clusters being highest in absolute orbital energies, and the highest charged clusters having lowest absolute orbital energies. Within the energy range spanned by these four orbitals among these clusters, the DFTB orbital energy levels reproduce TPSS orbital energy levels very well, as indicated by the dashed lines to mark the changes of electronic structure with molecular structure and total charge. Individual orbital energy level shifts are appreciable and may reach ≈0.9 eV (see Table S6†), with a general trend towards lower energies. Regarding HLGs, DFTB/auorg^α′ tends to underestimate those for the singly charged [Au₇(PPh₃)₇]⁺ (BIXZAK), the doubly charged [Au₈(PPh₃)₈]²⁺ (OPAUPF), and the triply charged [Au₁₁(PMePh₂)₁₀]³⁺ (ZUCMEP) clusters. In these particular three clusters, the DFTB HOMO energy shifts are −0.28, −0.32, and −0.40 eV for Au₇, Au₈ and Au₁₁ respectively, while their LUMO shifts are −0.63, −0.89, −0.92 eV, obviously quantitatively larger. This imbalance results in their smaller HOMO–LUMO gap values. Nonetheless, the HLG values of DFTB/auorg^α′ are quantitatively close to the TPSS calculated values with deviations no greater than 0.6 eV.

Fig. 7 Energy level diagram for the frontier orbitals of various clusters as calculated by (A) TPSS/def2-SVP and (B) DFTB2/auorg^α′, [Au₁₁(PMePh₂)₁₀]³⁺_# denotes [Au₁₁(PMePh₂)₁₀]³⁺ (C_3v), [Au₁₁(PMePh₂)₁₀]³⁺_* denotes [Au₁₁(PMePh₂)₁₀]³⁺ (D_4d). Dashed lines are included to guide the eye.

Fig. S7, in the ESI,† displays the HOMO and LUMO orbital shapes for [Au₆(dppp)₄]²⁺, [Au₇(PPh₃)₇]⁺ (BIXZAK), [Au₈(PPh₃)₈]²⁺ (OPAUPF) and [Au₉(PPh₃)₈]³⁺ (MIVPOX-D_2h) clusters for the TPSS and DFTB2/auorg^α′ methods. Despite some inevitable quantitative differences in the orbital topology, the amplitudes of the HOMO and LUMO orbitals obtained with the DFTB2/auorg^α′ method are very close to the ones obtained by TPSS/def2-SVP. Noticeable qualitative differences are the LUMOs of Au₇ and Au₈ where DFTB both overestimated the contribution from the gold core and less contributions from the phosphine ligands, possibly providing a rationale for the previously discussed large LUMO energy shifts of Au₇ and Au₈. The HOMO and LUMO plots of DFTB2/auorg^α and DFTB2/auorg^χ′ are shown in Fig. S8 and S9 in the ESI.† Both of these two sets of molecular orbitals have similar shapes with DFTB2/auorg^α′ orbitals; with DFTB2/auorg^α also overestimating contribution from the gold core while DFTB2/auorg^χ′ underestimates the gold and overestimates the contribution from phosphine ligands in Au₇. We note that these differences within the different DFTB parameter sets originate mostly from the different choices of the virtual orbital energies (see Table 1) and a minor degree from the different geometries.

3.5 IR spectra of [Au₆(dppp)₄]²⁺, [Au₈(PPh₃)₈]²⁺ and [Au₉(PPh₃)₈]³⁺ clusters

The IR spectra of [Au₆(dppp)₄]²⁺ (BOTSOS), [Au₈(PPh₃)₈]²⁺ (OPAUPF), and [Au₉(PPh₃)₈]³⁺ (MIVPOX-D_2h) clusters have been studied previously both experimentally and computationally using DFT¹⁹ and therefore represent good choices for the performance of our DFTB parameters for predicting vibrational spectra. The calculated DFTB IR spectra were evaluated by comparing them with experimental data and DFT-calculated IR, see Fig. 8 and S10 in the ESI.† The experimental far-IR spectra of the three clusters consists of two main parts: gold core distortion between the range 90–250 cm⁻¹ and Au–P modes above 400 cm⁻¹. Both PBE/def2-SVP and DFTB2/auorg^α′-calculated IR spectra match the previously reported experimental and M06/LANL2DZ predicted spectra, having the same mode description for each vibration. Also, the DFTB and PBE normalized spectral shapes are in good agreement, especially concerning the main peaks involving Au–P interactions. A noticeable difference between the experimental and calculated spectra of Au₈ and Au₉ is the presence of an experimentally distinct peak at 398 cm⁻¹. While all DFT and DFTB methods are able to predict peaks that are comparable with the position of this feature, the calculated intensities in DFTB are considerably weaker. Nevertheless, both methods have shown that this peak can be attributed to phenyl group twisting vibrations. Another difference is the presence of a peak at 317 cm⁻¹ in the Au₈ and Au₉ DFTB calculated spectra that is missing in the DFT spectra. The DFTB peak matches a seemingly broad experimental peak in the same region. In the previous study,¹⁹ this peak was not assigned to any vibrational mode, as DFT was not able to predict this feature. DFTB revealed that these peaks near 317 cm⁻¹ of Au₈ are related to vibrations with weaker intensity assigned to a phenyl rocking vibration that is attached to protruding Au₁–P₈. In the case of Au₉, this peak has several contributing vibrations that are also assigned to a phenyl rocking of the complex's several phosphine ligands. The summary of contributing transitions for all clusters are given in Tables S7–S9 in the ESI.† We conclude that DFTB is very useful in the prediction of vibrational spectra for nanoscale gold clusters with stabilizing ligands; this had not been discussed in the previous auorg parameterization work.^58,59

Fig. 8 Experimental (in black),¹⁹ computed DFTB2/auorg^α′ (in red), PBE/def2-SVP (in blue) and M06/LANL2DZ (in cyan) far-IR spectra for (A) [Au₆(dppp)₄]²⁺, (B) [Au₈(PPh₃)₈]²⁺, and (C) [Au₉(PPh₃)₈]³⁺ clusters. The additional red and blue dashed lines are for the scaled up plots of the region 100–400 cm⁻¹ for DFTB2/auorg^α′ and PBE/def2-SVP.

3.6 Chemisorption of PH₃ on the Au (111) surface

To evaluate the transferability of the new DFTB parameters, we explored the energy landscape of PH₃ chemisorption on the Au (111) surface and compared it to DFT with the PBE functional, which is popular in the solid state community to simulate surface processes. The computed DFTB2/auorg^α, DFTB2/auorg^α′, DFTB2/auorg^χ′ and PBE energy landscapes are shown in Fig. 9. In this Au (111) surface model, there are three important adsorption locations: top, hollow fcc, and hollow hcp (see Fig. S11 in the ESI†). The corresponding top and hollow adsorption sites are marked in Fig. 9 for reference. This test is particularly useful for evaluating the transferability of the new parameters, as the number of Au atoms coordinating with the P atom of PH₃ molecule varies from one at the top site, to two at the bridge sites (see Fig. S11 in the ESI†), to three at the hollow sites. Fig. 9 shows that the adsorption energies vary in a range of approximately 5 kcal mol⁻¹, and both auorg^α′ and auorg^χ′ parameters are able to reproduce qualitatively the relative PBE energy landscape. The top site is the most preferred adsorption location of PH₃ adsorption on the Au (111) surface according to PBE, DFTB2/auorg^α′ and DFTB2/auorg^χ′. However, DFTB2/auorg^α predicts the most preferred adsorption location is somewhere in the middle of top, hollow fcc, and hollow hcp positions.

Fig. 9 Energy landscape in kcal mol⁻¹ of PH₃ adsorption on the Au (111) surface obtained at the PBE and DFTB2 methods. The energy profiles show the three important adsorption locations of the Au (111) surface for clear comparison.

Considering the absolute binding energy, all three DFTB parameter sets underestimate the Au–P interaction energy compared to the PBE method. Among them, the underestimation is smallest with the deviation ΔE ≈ 4.6–7.4 kcal mol⁻¹ for auorg^α′, the deviation ΔE ≈ 7.6–9.8 kcal mol⁻¹ for auorg^χ′, and largest with ΔE ≈ 8.9–9.8 for auorg^α. The trends in the deviations can again be rationalized by correlating them with the orbital energies of the virtual Au 6p and P 3d orbitals: the larger the upward shift, the greater the underbinding. This test shows that DFTB2/auorg^α′ is the best option to study phosphine ligands chemisorbed on Au surfaces, and that our parameter has difficulties for the interaction of phosphines with flat Au surfaces, as mentioned above in Section 3.3.

3.7 Case study for Au₁₀₈S₂₄(PPh₃)₁₆: geometry, ligand binding energy, frontier orbitals, and IR spectra

In order to demonstrate the capability of the new DFTB parameters for the study of large-scale phosphine-stabilized gold clusters, we performed a case study for the recently synthesized metalloid complex Au₁₀₈S₂₄(PPh₃)₁₆.²⁵ The core of the metalloid consists of a Russian doll motif with an Au₆ inner octahedron enclosed by a second-shell Au₃₈ octahedron. This core is then capped at its six tips by novel Au₄S₄ planar rings, which are then connected by 4 more gold atoms, resulting in a (Au₂S)₂₄ outer shell motif. Finally, sixteen Au(PPh₃) groups complete this metalloid and are bound to each side of the Au₄₄ octahedral Russian-doll core. The available experimental crystal structure (CSD code: DAFLOO) provides a clear picture of Au, S and P atomic positions in 3D. However, it does not depict the positions of C atoms well due to overlapping rings of the triphenylphosphine groups, with unreasonable C–C distances in the CIF file, sometimes being shorter than 0.2 Å. To prepare the initial geometrical inputs for our DFTB geometry optimization, we constructed the initial structure based on the experimental positions of Au, S and P atoms, and then added manually the sixteen PPh₃ ligands. The geometry was then optimized using the DFTB2/auorg^α′ method. Fig. 10 compares the final DFTB optimized geometry to the initial, experimental crystal structure. DFTB2/auorg^α′ predicts the averaged and normalized ligand binding energies of Au₁₀₈S₂₄(PPh₃)₁₆ to be −63.2 kcal mol⁻¹, in good agreement with TPSS/def2-SVP single point energies performed using our DFTB optimized geometries, which predicts the averaged ligand binding energy to be −72.9 kcal mol⁻¹. We note that a single point energy + gradient calculation on 16 “Intel Xeon E5-2697 2.30 GHz CPU cores” took 39 935 s (11.09 h) with DFT and only 27 s with DFTB, corresponding to a speed-up factor of nearly 1500×. The DFTB method is therefore the only practical option if one considers routinely carrying out geometry optimizations of large nanoscale gold clusters, or a more challenging molecular dynamics (MD) simulation with a few million steps to reach nanosecond time scales.

Fig. 10 X-ray experimental structure of Au₁₀₈S₂₄(PPh₃)₁₆ taken from Cambridge Crystallographic Database (CSD code DAFLOO) with overlapping phenyl rings and ultra short ≈0.2 Å C–C bond lengths between triphenylphosphine groups (left panel), and DFTB2/auorg^α′ optimized structure without overlapping phenyl rings (right panel).

Computationally, the electronic structure of this cluster has been approximately investigated previously at the BP86/TZVPP level of theory, using a simplified model in which the triphenylphosphine ligands were replaced by PH₃ groups, and the entire structure assumed to take T_d point group symmetry.

The HOMO–LUMO gap of this simplified Au₁₀₈S₂₄(PH₃)₁₆ cluster was found to be 0.68 eV. Fig. 11 shows the DFTB2/auorg^α′ HOMO and LUMO plots of Au₁₀₈S₂₄, Au₁₀₈S₂₄(PH₃)₁₆, and Au₁₀₈S₂₄(PPh₃)₁₆ clusters. The calculated HOMO–LUMO gaps for these three clusters are 0.200 eV, 0.681 eV, and 0.634 eV, respectively. The increase in the HLG from the Au₁₀₈S₂₄ core to fully ligated cluster denotes that the phosphine groups stabilize the core, much like in Au₇₀S₂₀(PPh₃)₁₂.²⁶ The quantitative difference in the energy gaps between Au₁₀₈S₂₄(PH₃)₁₆ and Au₁₀₈S₂₄(PPh₃)₁₆ clusters are smaller than 0.05 eV. However, it is worth noting that partial density of state (PDOS) plots, see Fig. 12, portray that there is an undeniable difference of sulfur-orbital contribution in the bonding of the Au₁₀₈S₂₄(PH₃)₁₆ and Au₁₀₈S₂₄(PPh₃)₁₆ clusters. In addition, both HOMO and LUMO are shifted by ≈−0.80 eV when PPh₃ ligands are replaced by PH₃ ligands. In the case of the HLG, employing the simplified model Au₁₀₈S₂₄(PH₃)₁₆ yields a similar value to the HLG of the Au₁₀₈S₂₄(PPh₃)₁₆ cluster, however, the Fermi level has shifted dramatically by −0.8 eV! Our calculations therefore indicate that the use of a simplified model is inadequate and can indeed significantly affect the electronic properties of nanoclusters.

Fig. 11 HOMO and LUMO plots of Au₁₀₈S₂₄, Au₁₀₈S₂₄ (PH₃)₁₆, and Au₁₀₈S₂₄(PPh₃)₁₆ clusters as calculated by DFTB2/auorg^α′; isosurface value = 0.015 a.u.

Fig. 12 Density of states (DOS) and partial density of states (PDOS) of Au₁₀₈S₂₄(PH₃)₁₆ (top) and Au₁₀₈S₂₄(PPh₃)₁₆ (bottom) clusters as calculated by DFTB2/auorg^α′.

The calculation of IR spectra of large-scale systems within the normal mode approximation, requiring the calculation of the Hessian matrix, is beyond the capability of most contemporary DFT codes on current computer systems within a reasonable amount of time. In this case, DFTB can be particularly useful as it has been shown in the benchmark sections that DFTB can reproduce the experimental IR spectra very well. In this case study, we carried out DFTB2/auorg^α′ Hessian calculations to predict the IR spectrum for the Au₁₀₈S₂₄(PPh₃)₁₆ cluster. Since the performance of DFTB in simulating infrared spectra of thiolate-ligated gold nanoclusters has not been tested before, in this work, we validated the accuracy of DFTB in predicting IR spectra for various thiolated gold clusters. We compare the DFTB simulated IR spectra of Au₄(SCH₃)₄, [Au₂₅(SCH₃)₁₈]⁻ clusters, and six Au₁₈ clusters protected with various types of ligands to the corresponding DFT and experimental spectra in the case of Au₁₈(S-c-C₆H₁₁)₁₄.^23,31,32 The comparison is presented in the ESI, Fig. S13 and S14.† In summary, similar to phosphine-stabilized nanoclusters in Section 3.5, the DFTB-calculated IR spectra for thiolate-protected gold clusters agree well with previous DFT simulations and experimental IR spectra. These new results, in combination with the results presented in Section 3.5, indicate that DFTB is a reliable method in predicting IR spectra of thiolate- and phosphine-stabilized gold nanoclusters. Fig. 13 shows the DFTB calculated far-IR spectrum of Au₁₀₈S₂₄(PPh₃)₁₆. We focus our attention on the far-IR region because it is of particular interest in understanding the core vibrations of Au, S, and P atoms. The calculated IR spectrum has several large peaks in the 500–600 cm⁻¹ region which are attributed to PPh₃ distortions. The multiple peaks in the 200–350 cm⁻¹ region are caused by various normal modes of both Au–P and Au–S vibrations (see Tables S10–S12 in the ESI†), the peaks labelled in the figure are the normal modes related to the novel Au₄S₄ planar rings: (1) Au–S stretching at 241.1 cm⁻¹, (2) S–Au–S–Au–S symmetric stretching at 277.2 cm⁻¹ and (3) S–Au–S symmetric stretch at 517.05 cm⁻¹. The inset of Fig. 13 shows several peaks <150 cm⁻¹ which are mostly attributed to Au₄₄ core distortion.

Fig. 13 Predicted far-IR spectra for Au₁₀₈S₂₄(PPh₃)₁₆ clusters calculated using DFTB2/auorg^α′. The inset of the figure shows the scaled up plots of the region 0–550 cm⁻¹ with labelled peaks pertaining to normal modes of Au₄S₄ planar rings. The additional black dashed line in the inset is for the scaled-up plot of the region 0–150 cm⁻¹.

Since the novel Au₄S₄ planar ring motif of Au₁₀₈S₂₄(PPh₃)₁₆ has not yet been studied using IR spectroscopy, there are no experimental validations so far about its normal modes. The calculated Au–S stretches found in this planar ring motif are compared with that of the Au₄(SCH₃)₄, [Au₂₅(SCH₃)₄]⁻ and [Au₃₇(PPh₃)₁₀(SR)₁₀Cl₂]⁺ clusters.^31–33 The S–Au–S–Au–S asymmetric and symmetric stretching values are found to be 274 cm⁻¹ and 277 cm⁻¹ (see Peak 2 of the inset in Fig. 13). The S–Au–S–Au–S symmetric stretching value is similar to the DFT-calculated breathing mode (293 cm⁻¹) of Au₄(SCH₃)₄ cluster. Other IR studies of thiolate-containing Au nanoclusters such as [Au₂₅(SCH₃)₁₈]⁻ and [Au₃₇(PPh₃)₁₀(SR)₁₀Cl₂]⁺ have reported their Au–S symmetric stretches at 293 cm⁻¹ and 239 cm⁻¹, respectively, which are relatively close to our DFTB-calculated Au–S stretch values; 241 cm⁻¹ as shown as Peak 1 in the inset of Fig. 13, and 287 cm⁻¹ as presented in Tables S11 and S12,† for the Au₁₀₈S₂₄(PPh₃)₁₆ cluster. Peak 3 in the inset of Fig. 13 shows another peak at 514 cm⁻¹, which is another characteristic S–Au–S symmetric stretch. One can depict the similarity of the Au–S stretches in the 200–300 cm⁻¹ regions with that of other sulfur-capped Au clusters.^{20,23,31–33,96,97} Based on the good agreement between our DFTB calculated IR spectra with DFT calculated as well as the experimental spectra for gold clusters in the benchmark, the predicted IR spectrum of Au₁₀₈S₂₄(PPh₃)₁₆ is reliable and provides a useful fingerprint to identify this cluster in future studies.

4 Conclusions

Parameters for the density-functional tight-binding (DFTB) method were generated to describe gold–phosphorus interactions for simulations of phosphine-stabilized nanoscale gold clusters. We build on preceding works reporting second-order DFTB (DFTB2) parameters for hybrid gold–thiolates compounds,^58,59 which were themselves based on the mio parameter set.^51,52,57 The present effort thus expands the applicability of the DFTB2 method to organometallic gold complexes with ligands containing the full set of chemical elements contained in the mio parameter set: H, C, N, O, P, and S. In the construction of the repulsive potentials we considered three different combinations of the gold 6p and phosphorus 3d virtual atomic orbital energies. The performance of our parameters was evaluated using density functional theory (DFT) geometries, ligand binding and cluster isomerization energies, ligand dissociation potential energy curves, molecular orbital energies, and simulated far-IR spectra. We further compare predicted geometries with X-ray crystallographic structures and experimental far-IR spectra, and evaluate parameter transferability for phosphine chemisorption on a gold surface.

In general it is found that the absolute ligand binding energies increase with decreased virtual orbital energies, i.e. auorg^α < auorg^χ′ ≈ auorg^α′. The total ligand binding energy increases both with cluster and ligand size. This means that for a given gold cluster, the variation in the ligand binding energy in DFTB needs to be almost completely described by the electronic energy, as the Au–P bond distances only marginally change with different ligands. However, the current DFTB electronic parameters are not flexible enough to accurately describe the variations in the electronic ligand binding energies. Similarly, the same trend holds for increasing gold cluster size. According to the benchmark results, we found that DFTB2/auorg^α is the best option to study small ligands and small Au clusters. DFTB2/auorg^α′ is the best option to study large ligands, large Au clusters, or Au surfaces. For the moderate-sized ligands or gold clusters, both DFTB2/auorg^α and DFTB2/auorg^α′ are good options to be used. The performance of auorg^χ′ for Au–P interactions is similar to that of auorg^α′ because the effects of shifting the Au 6p orbital energy from −0.02786 to −0.00001 hartree is less significant than shifting the P 3d orbital energy from 0.12044 to 0.52044 hartree. However, the effect of shifting the Au 6p orbital energy has not been investigated for Au–(C, H, N, O and S) interactions, thus the use of auorg^α′ is more preferable over auorg^χ′. Besides consideration of cluster and ligand sizes, the geometrical environment of the Au–P bond also influences the performance of our parameters. Surface-like Au atoms present a challenge while more pyramidalized Au atoms are better described. Possible reasons are the minimum basis set or the missing multipolar charge contribution. Such situations are rare, however, and our auorg^α′ parameter performs overall well in theoretical studies of relevant nanoscale gold clusters that experimentally feature large capping ligands. In particular, our DFTB parameterization enables the simulation of ligand dissociation, reliably predicting how post-treatment can be done in experiment without inducing concomitant side-effects such as agglomeration. The auorg^α′ parameter set will be publicly distributed via the http://www.dftb.org website. To switch from DFTB2/auorg^α′ to DFTB2/auorg^α, one can simply modify the P 3d-orbital energy according to Table 1.

We employed the new DFTB parameters to determine the geometric structure of an Au₁₀₈S₂₄(PPh₃)₁₆ nanocluster²⁵ and investigate its molecular and electronic structure. The optimized Cartesian coordinates of Au₁₀₈S₂₄(PPh₃)₁₆ are provided in the ESI.† We found that for large ligand-protected gold clusters, using a simplified model in simulations is inadequate and can substantially affect the bonding and electronic structures of the clusters. Finally, we predicted the IR spectrum of the cluster. Both optimized geometric structure and the simulated IR spectrum of the Au₁₀₈S₂₄(PPh₃)₁₆ nanocluster will serve as useful reference for future studies.

DFTB approaches with appropriately developed parameters for the description of molecular and electronic structure and energetics as well as vibrational spectra are promising to provide insight as to how to tune catalyst reactivity for both product selectivity and reaction specificity. Further development of DFTB parameters for the interaction of gold clusters with various substrate surfaces will advance the development of transformative catalytic systems. This work is vital to herald the promised potential of unprecedented reaction tunability using cluster-based catalysts.

Notice of copyright

This manuscript has been authored by UT-Battelle, LLC under Contract No. DE-AC05-00OR22725 with the U.S. Department of Energy. The United States Government retains and the publisher, by accepting the article for publication, acknowledges that the United States Government retains a non-exclusive, paid-up, irrevocable, worldwide license to publish or reproduce the published form of this manuscript, or allow others to do so, for United States Government purposes. The Department of Energy will provide public access to these results of federally sponsored research in accordance with the DOE Public Access Plan (http://energy.gov/downloads/doe-public-access-plan).

Conflicts of interest

There are no conflicts to declare.

Acknowledgements

VQV acknowledges support by an Energy Science and Engineering Fellowship of the Bredesen Center for Interdisciplinary Research and Graduate Education at the University of Tennessee, Knoxville. JMLM acknowledges financial support from Adelaide Scholarship International (ASI) and Future Fuels Cooperative Research Center (FFCRC). SI acknowledges support by the Laboratory Directed Research and Development (LDRD) Program of Oak Ridge National Laboratory. ORNL is managed by UT-Battelle, LLC, for DOE under Contract No. DE-AC05-00OR22725. This research used resources of the National Energy Research Scientific Computing Center (NERSC), a U.S. Department of Energy Office of Science User Facility operated under Contract No. DE-AC02-05CH11231, as well as the Compute and Data Environment for Science (CADES) at Oak Ridge National Laboratory. Supercomputing resources at the University of Adelaide were provided by the Phoenix HPC service. BGS performed work at the Center for Nanophase Materials Sciences which is a US Department of Energy Office of Science User Facility. We would like to thank Junda Li for insightful assistance and comments during the initial stage prior to this concrete work.

Notes and references

1. W. Jahn, J. Struct. Biol., 1999, 127, 106–112.
2. R. Meyer, C. Lemire, S. K. Shaikhutdinov and H. J. Freund, Gold Bull., 2004, 37, 72–124.
3. Y. Du, H. Sheng, D. Astruc and M. Zhu, Chem. Rev., 2020, 120, 526–622.
4. T. G. Schaaff, M. N. Shafigullin, J. T. Khoury, I. Vezmar, R. L. Whetten, W. G. Cullen, P. N. First, C. Gutiérrez-Wing, J. Ascensio and M. J. Jose-Yacamán, J. Phys. Chem. B, 1997, 101, 7885–7891.
5. A. Heuer-Jungemann, N. Feliu, I. Bakaimi, M. Hamaly, A. Alkilany, I. Chakraborty, A. Masood, M. F. Casula, A. Kostopoulou, E. Oh, K. Susumu, M. H. Stewart, I. L. Medintz, E. Stratakis, W. J. Parak and A. G. Kanaras, Chem. Rev., 2019, 119, 4819–4880.
6. J. Kilmartin, R. Sarip, R. Grau-Crespo, D. Di Tommaso, G. Hogarth, C. Prestipino and G. Sankar, ACS Catal., 2012, 2, 957–963.
7. M. Haruta, Gold Bull., 2004, 37, 27–36.
8. V. R. Choudhary and D. K. Dumbre, Top. Catal., 2009, 52, 1677–1687.
9. G. E. Johnson, T. Priest and J. Laskin, Chem. Sci., 2014, 5, 3275–3286.
10. D. P. Anderson, J. F. Alvino, A. Gentleman, H. A. Qahtani, L. Thomsen, M. I. J. Polson, G. F. Metha, V. B. Golovko and G. G. Andersson, Phys. Chem. Chem. Phys., 2013, 15, 3917–3929.
11. D. P. Anderson, R. H. Adnan, J. F. Alvino, O. Shipper, B. Donoeva, J.-Y. Ruzicka, H. Al Qahtani, H. H. Harris, B. Cowie, J. B. Aitken, V. B. Golovko, G. F. Metha and G. G. Andersson, Phys. Chem. Chem. Phys., 2013, 15, 14806–14813.
12. S. Vajda and M. G. White, ACS Catal., 2015, 5, 7152–7176.
13. S. Tesana, G. F. Metha, G. G. Andersson, C. Ridings and V. Golovko, Int. J. Nanotechnol., 2018, 15, 669–675.
14. G. Krishnan, N. Eom, R. M. Kirk, V. B. Golovko, G. F. Metha and G. G. Andersson, J. Phys. Chem. C, 2019, 123, 6642–6649.
15. K. D. M. Weerawardene, H. Häkkinen and C. M. Aikens, Annu. Rev. Phys. Chem., 2018, 69, 205–229.
16. K. A. Parrish, M. King, M. R. Ligare, G. E. Johnson and H. Hernández, Phys. Chem. Chem. Phys., 2019, 21, 1689–1699.
17. M. Zhu, C. M. Aikens, F. J. Hollander, G. C. Schatz and R. Jin, J. Am. Chem. Soc., 2008, 130, 5883–5885.
18. S. Goel, K. A. Velizhanin, A. Piryatinski, S. Tretiak and S. A. Ivanov, J. Phys. Chem. Lett., 2010, 1, 927–931.
19. J. F. Alvino, T. Bennett, D. Anderson, B. Donoeva, D. Ovoshchnikov, R. H. Adnan, D. Appadoo, V. Golovko, G. Andersson and G. F. Metha, RSC Adv., 2013, 3, 22140–22149.
20. B. Varnholt, P. Oulevey, S. Luber, C. Kumara, A. Dass and T. Bürgi, J. Phys. Chem. C, 2014, 118, 9604–9611.
21. A. Fernando, K. L. D. M. Weerawardene, N. V. Karimova and C. M. Aikens, Chem. Rev., 2015, 115, 6112–6216.
22. F. Muniz-Miranda, D. Presti, M. C. Menziani and A. Pedone, Theor. Chem. Acc., 2016, 135, 5–13.
23. A. Tlahuice-Flores, Prog. Nat. Sci.: Mater. Int., 2016, 26, 510–515.
24. S. Yamazoe, S. Matsuo, S. Muramatsu, S. Takano, K. Nitta and T. Tsukuda, Inorg. Chem., 2017, 56, 8319–8325.
25. S. Kenzler, C. Schrenk and A. Schnepf, Angew. Chem., Int. Ed., 2017, 56, 393–396.
26. S. Kenzler, C. Schrenk, A. R. Frojd, H. Häkkinen, A. Z. Clayborne and A. Schnepf, Chem. Commun., 2018, 54, 248–251.
27. S.-S. Zhang, L. Feng, R. D. Senanayake, C. M. Aikens, X.-P. Wang, Q.-Q. Zhao, C.-H. Tung and D. Sun, Chem. Sci., 2018, 9, 1251–1258.
28. J. C. Burgos, S. M. Mejía and G. F. Metha, ACS Omega, 2019, 4, 9169–9180.
29. Z. Tian, Y. Xu and L. Cheng, Nanomaterials, 2019, 9, 1132–1143.
30. A. Tlahuice-Flores, R. L. Whetten and M. Jose-Yacaman, J. Phys. Chem. C, 2013, 117, 20867–20875.
31. A. Tlahuice-Flores, Mol. Simul., 2013, 39, 428–431.
32. A. Tlahuice-Flores, R. L. Whetten and M. Jose-Yacaman, J. Phys. Chem. C, 2013, 117, 12191–12198.
33. O. Palacios-Álvarez and A. Tlahuice-Flores, J. Raman Spectrosc., 2019, 50, 52–62.
34. A. Wing-Bocanegra and A. Tlahuice-Flores, Phys. Chem. Chem. Phys., 2019, 21, 23855–23864.
35. A. V. Akimov and O. V. Prezhdo, Chem. Rev., 2015, 115, 5797–5890.
36. O. D. Haeberlen and N. Roesch, J. Phys. Chem., 1993, 97, 4970–4973.
37. P. Schwerdtfeger and P. D. W. Boyd, Inorg. Chem., 1992, 31, 327–329.
38. R. L. DeKock, E. J. Baerends, P. M. Boerrigter and R. Hengelmolen, J. Am. Chem. Soc., 1984, 106, 3387–3396.
39. D. M. P. Mingos, T. Slee and L. Zhenyang, Chem. Rev., 1990, 90, 383–402.
40. L. W. Chung, W. M. C. Sameera, R. Ramozzi, A. J. Page, M. Hatanaka, G. P. Petrova, T. V. Harris, X. Li, Z. Ke, F. Liu, H.-B. Li, L. Ding and K. Morokuma, Chem. Rev., 2015, 115, 5678–5796.
41. T. Vreven and K. Morokuma, J. Comput. Chem., 2000, 21, 1419–1432.
42. B. G. Donoeva, D. S. Ovoshchnikov and V. B. Golovko, ACS Catal., 2013, 3, 2986–2991.
43. W. Thiel, Wiley Interdiscip. Rev.: Comput. Mol. Sci., 2014, 4, 145–157.
44. T. Husch, A. C. Vaucher and M. Reiher, Int. J. Quantum Chem., 2018, 118, 25799–25831.
45. Q. Cui and M. Elstner, Phys. Chem. Chem. Phys., 2014, 16, 14368–14377.
46. A. S. Christensen, T. Kubař, Q. Cui and M. Elstner, Chem. Rev., 2016, 116, 5301–5337.
47. S. Grimme, C. Bannwarth and P. Shushkov, J. Chem. Theory Comput., 2017, 13, 1989–2009.
48. C. Bannwarth, S. Ehlert and S. Grimme, J. Chem. Theory Comput., 2019, 15, 1652–1671.
49. D. Porezag, T. Frauenheim, T. Köhler, G. Seifert and R. Kaschner, Phys. Rev. B: Condens. Matter Mater. Phys., 1995, 51, 12947–12957.
50. G. Seifert, D. Porezag and T. Frauenheim, Int. J. Quantum Chem., 1996, 58, 185–192.
51. M. Elstner, D. Porezag, G. Jungnickel, J. Elsner, M. Haugk, T. Frauenheim, S. Suhai and G. Seifert, Phys. Rev. B: Condens. Matter Mater. Phys., 1998, 58, 7260–7268.
52. M. Gaus, Q. Cui and M. Elstner, J. Chem. Theory Comput., 2011, 7, 931–948.
53. T. A. Niehaus and F. Della Sala, Phys. Status Solidi, 2012, 249, 237–244.
54. M. Elstner, J. Phys. Chem. A, 2007, 111, 5614–5621.
55. M. Gaus, Q. Cui and M. Elstner, Wiley Interdiscip. Rev.: Comput. Mol. Sci., 2014, 4, 49–61.
56. J. C. Slater and G. F. Koster, Phys. Rev., 1954, 94, 1498–1524.
57. T. Niehaus, M. Elstner, T. Frauenheim and S. Suhai, J. Mol. Struct.: THEOCHEM, 2001, 541, 185–194.
58. A. Fihey, C. Hettich, J. Touzeau, F. Maurel, A. Perrier, C. Köhler, B. Aradi and T. Frauenheim, J. Comput. Chem., 2015, 36, 2075–2087.
59. L. F. L. Oliveira, N. Tarrat, J. Cuny, J. Morillo, D. Lemoine, F. Spiegelman and M. Rapacioli, J. Phys. Chem. A, 2016, 120, 8469–8483.
60. M. Elstner and G. Seifert, Philos. Trans. R. Soc., A, 2014, 372, 20120483–20120494.
61. R. S. Mulliken, J. Chem. Phys., 1955, 23, 1833–1840.
62. J. P. Perdew, K. Burke and M. Ernzerhof, Phys. Rev. Lett., 1996, 77, 3865–3868.
63. M. Gaus, C.-P. Chou, H. Witek and M. Elstner, J. Phys. Chem. A, 2009, 113, 11866–11881.
64. C.-P. Chou, Y. Nishimura, C.-C. Fan, G. Mazur, S. Irle and H. A. Witek, J. Chem. Theory Comput., 2016, 12, 53–64.
65. V. Q. Vuong, J. Akkarapattiakal Kuriappan, M. Kubillus, J. J. Kranz, T. Mast, T. A. Niehaus, S. Irle and M. Elstner, J. Chem. Theory Comput., 2018, 14, 115–125.
66. M. Gaus, X. Lu, M. Elstner and Q. Cui, J. Chem. Theory Comput., 2014, 10, 1518–1537.
67. J. Tao, J. P. Perdew, V. N. Staroverov and G. E. Scuseria, Phys. Rev. Lett., 2003, 91, 146401–146404.
68. K. P. Kepp, J. Phys. Chem. A, 2017, 121, 2022–2034.
69. S. Goel, K. A. Velizhanin, A. Piryatinski, S. Tretiak and S. A. Ivanov, J. Phys. Chem. Lett., 2010, 1, 927–931.
70. F. Weigend and R. Ahlrichs, Phys. Chem. Chem. Phys., 2005, 7, 3297–3305.
71. S. Grimme, J. Antony, S. Ehrlich and H. Krieg, J. Chem. Phys., 2010, 132, 154104.
72. S. Grimme, S. Ehrlich and L. Goerigk, J. Comput. Chem., 2011, 32, 1456–1465.
73. J. Řezáč and P. Hobza, J. Chem. Theory Comput., 2012, 8, 141–151.
74. D. Andrae, U. Häußermann, M. Dolg, H. Stoll and H. Preuß, Theor. Chim. Acta, 1990, 77, 123–141.
75. F. Weigend, Phys. Chem. Chem. Phys., 2006, 8, 1057–1065.
76. F. Neese, Wiley Interdiscip. Rev.: Comput. Mol. Sci., 2018, 8, 1327–1332.
77. B. Hourahine, B. Aradi, V. Blum, F. Bonafé, A. Buccheri, C. Camacho, C. Cevallos, M. Y. Deshaye, T. Dumitrică, A. Dominguez, S. Ehlert, M. Elstner, T. van der Heide, J. Hermann, S. Irle, J. J. Kranz, C. Köhler, T. Kowalczyk, T. Kubař, I. S. Lee, V. Lutsker, R. J. Maurer, S. K. Min, I. Mitchell, C. Negre, T. A. Niehaus, A. M. N. Niklasson, A. J. Page, A. Pecchia, G. Penazzi, M. P. Persson, J. Řezáč, C. G. Sánchez, M. Sternberg, M. Stöhr, F. Stuckenberg, A. Tkatchenko, V. W.-z. Yu and T. Frauenheim, J. Chem. Phys., 2020, 152, 124101–124119.
78. J. W. A. Van der Velden, J. J. Bour, J. J. Steggerda, P. T. Beurskens, M. Roseboom and J. H. Noordik, Inorg. Chem., 1982, 21, 4321–4324.
79. C. E. Briant, K. P. Hall, D. M. P. Mingos and A. C. Wheeler, J. Chem. Soc., Dalton Trans., 1986, 687–692.
80. M. Schulz-Dobrick and M. Jansen, Angew. Chem., Int. Ed., 2008, 47, 2256–2259.
81. J. W. A. Van der Velden, J. J. Bour, W. P. Bosman and J. H. Noordik, Inorg. Chem., 1983, 22, 1913–1918.
82. F. Wen, U. Englert, B. Gutrath and U. Simon, Eur. J. Inorg. Chem., 2008, 106–111.
83. R. C. B. Copley and D. M. P. Mingos, J. Chem. Soc., Dalton Trans., 1996, 479–489.
84. X.-K. Wan, S.-F. Yuan, Z.-W. Lin and Q.-M. Wang, Angew. Chem., Int. Ed., 2014, 53, 2923–2926.
85. J. Chen, Q.-F. Zhang, T. A. Bonaccorso, P. G. Williard and L.-S. Wang, J. Am. Chem. Soc., 2014, 136, 92–95.
86. X.-K. Wan, J.-Q. Wang, Z.-A. Nan and Q.-M. Wang, Sci. Adv., 2017, 3, 1701823–1701828.
87. Y. Nishimoto and S. Irle, Chem. Phys. Lett., 2017, 667, 317–321.
88. G. M. J. Barca, C. Bertoni, L. Carrington, D. Datta, N. De Silva, J. E. Deustua, D. G. Fedorov, J. R. Gour, A. O. Gunina, E. Guidez, T. Harville, S. Irle, J. Ivanic, K. Kowalski, S. S. Leang, H. Li, W. Li, J. J. Lutz, I. Magoulas, J. Mato, V. Mironov, H. Nakata, B. Q. Pham, P. Piecuch, D. Poole, S. R. Pruitt, A. P. Rendell, L. B. Roskop, K. Ruedenberg, T. Sattasathuchana, M. W. Schmidt, J. Shen, L. Slipchenko, M. Sosonkina, V. Sundriyal, A. Tiwari, J. L. Galvez Vallejo, B. Westheimer, M. Włoch, P. Xu, F. Zahariev and M. S. Gordon, J. Chem. Phys., 2020, 152, 154102–154127.
89. P. E. Blöchl, Phys. Rev. B: Condens. Matter Mater. Phys., 1994, 50, 17953–17979.
90. G. Kresse and J. Furthmüller, Comput. Mater. Sci., 1996, 6, 15–50.
91. G. Kresse and J. Furthmüller, Phys. Rev. B: Condens. Matter Mater. Phys., 1996, 54, 11169–11186.
92. K. Spivey, J. I. Williams and L. Wang, Chem. Phys. Lett., 2006, 432, 163–166.
93. G. Shafai, S. Hong, M. Bertino and T. S. Rahman, J. Phys. Chem. C, 2009, 113, 12072–12078.
94. M. W. Walker, L. Shao and R. A. Volz, CVGIP: Image Understanding, 1991, 54, 358–367.
95. D. Haberer, C. E. Giusca, Y. Wang, H. Sachdev, A. V. Fedorov, M. Farjam, S. A. Jafari, D. V. Vyalikh, D. Usachov, X. Liu, U. Treske, M. Grobosch, O. Vilkov, V. K. Adamchuk, S. Irle, S. R. P. Silva, M. Knupfer, B. Büchner and A. Grüneis, Adv. Mater., 2011, 23, 4497–4503.
96. J. Petroski, M. Chou and C. Creutz, J. Organomet. Chem., 2009, 694, 1138–1143.
97. D. Stellwagen, A. Weber, G. L. Bovenkamp, R. Jin, J. H. Bitter and C. S. S. R. Kumar, RSC Adv., 2012, 2, 2276–2283.

---

Footnotes

† Electronic supplementary information (ESI) available: Fig. S1–S15, Tables S1–S12 and DFTB optimized Cartesian coordinates of Au₁₀₈S₂₄(PPh₃)₁₆ cluster and xyz trajectory files of its Au₄S₄ planar ring normal mode displacements. See DOI: 10.1039/d0sc04514d

‡ These authors contributed equally to this work.

---