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

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 Aun 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 PH3 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 Au108S24(PPh3)16.


Introduction
Atomically precise ligated gold clusters of nanometer dimension receive continued attention due to their unique catalytic properties 1-3 and well-dened discrete electronic energy levels 4 that potentially offer greater exibility 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][7][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. 9 Post-treatment is usually required to remove some (or all) of the ligands to allow interaction with the substrate or reactants. [10][11][12][13][14] To achieve control over such complex catalytic systems it is vitally important to understand the relationship between molecular and electronic structure, oen 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 oen employed in theoretical investigations, as they are capable to accurately describe electronic, geometrical, and vibrational structure of gold clusters and nanoparticles. [16][17][18][19][20][21][22][23][24][25][26][27][28][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 ngerprints to distinguish between bonding arrangements and orientations between gold atoms and ligands, and ligandligand interactions within clusters. 23,[30][31][32][33][34] Unfortunately, DFT calculations can become prohibitively expensive with system size, 35 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 simplied ligand models where e.g. triphenylphosphine (PPh 3 ) is replaced by phosphine (PH 3 ) or trimethylphosphine (PMe 3 ). 17,18,22,25,26,28,29,36 This, as with the complete ligand removal approach, has the additional benet that conformational searching is simplied, 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][37][38][39] and a computational truncation of the ligands will inuence the chemistry and therefore description of the catalytic properties in calculations. Integrated schemes such as ONIOM 40 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 denition 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, 41 which is oen computationally unfeasible. In practice therefore, integrated methods are oen 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 ne-tuning. Exacerbating this problem is the recent emergence of larger, so called metalloid, clusters such as Au 108 S 24 (PPh 3 ) 16 . 25 To make matters even worse, interactions of deposited gold clusters with substrate surfaces such as SiO 2 and TiO 2 may play an important role in the catalytic reaction, [10][11][12][13][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. 15 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 rst principles methods. [43][44][45][46][47][48] Among them, density-functional tightbinding (DFTB), 45,[49][50][51][52][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 twocenter 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 DFTcalculated 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 proles, (3) electronic structure, (4) vibrational normal modes of Au-P containing clusters, and (5) the adsorption of a PH 3 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 108 S 24 (PPh 3 ) 16 . 25 2 Methodology and computational details

Brief overview of DFTB2
A comprehensive review of DFTB methods can be found elsewhere 60 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. 51 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 r 0 and electron density uctuations Dr whereĤ 0 is the initial Hamiltonian constructed from the superposition of neutral atomic densities in a two-center approximation, 49 and the |J i i are occupied valence molecular orbitals (MOs), expanded as a linear combination of optimized pseudo-atomic orbitals |f m i. The optimization of these pseudoatomic valence orbitals and orbital densities for a given chemical element constitutes the determination of the electronic parameters. Dq A is a point charge 61 on atom A, and g AB Dq A Dq B represents the Coulomb interaction energy between the two point charges; 51 when A ¼ B, g AA is the chemical hardness or second derivative of the total energy with respect to the charge on atom A. g AB is dened as where 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. 62 In the framework of the two-center approximation, the Hamiltonian integrals hf m |Ĥ 0 |f n i and the overlap integrals hf m |f n i 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 rst 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 twocenter term that depends only on the chemical element type of atoms A and B and their interatomic distance r AB . 60 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. 59 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 ts 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 0 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][64][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 t 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.

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 denition of atomic connement potentials for generating pseudo-atomic orbitals |f m i as well as the connement 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. 51 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 a ), 58 and a modied version of "auorg" by Oliveira et al. (referred to as auorg c ). 59 The difference between these two parameter sets lies in the Au 6p-orbital energy; in the auorg a set it was taken as the true PBE orbital energy, while in the auorg c set it was empirically shied upward by z+0.0279 hartree. The main purpose for this orbital energy shi was to obtain improved values for cohesive energies of pure gold nanoclusters with respect to PBE. 59 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 a and auorg c 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 shied 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 shi 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 shi of the P 3d virtual orbital energy introduces signicant 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" shied value, and used our genetic algorithm (GA) optimization tool 65 to automatically generate repulsive potentials for the Au-P interaction with these orbital energy shis. 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 3 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. 59 we denote the original PBE Au 6p energy value by 'a' and its modied value by 'c'. The shied "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, ' 0 '. Hence, auorg a coupled with the original "mio" P is unchanged, while the auorg a and auorg c coupled with the new P 3d-orbital energy are referred to as auorg a 0 and auorg c 0 , respectively. Because auorg c exhibits the largest underestimation of the Au-P electronic binding energy (both Au 6p and P 3d virtual orbital energies are shied 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. The Au-P repulsive potentials were optimized on the basis of shell-resolved self-consistent charge electronic energies to t ligand binding energies and forces for the training set listed in Table 2 using our in-house genetic algorithm (GA)-based parameterization tool. 65 Aer 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 dened as the tness of the parameter sets with respect to DFTB values and reference data according to the formula are weight factors of the i th binding energy and the force of the i th structure, respectively, N eq is the number of tting 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 shied the reference ligand binding energies of PH 3 by À24 kcal mol À1 for auorg a 0 and À22 kcal mol À1 for the auorg c 0 to better reproduce ligand binding energies of larger ligands and the adsorption of PH 3 on the Au (111) surface. Surprisingly, auorg a 0 and auorg c 0 optimized Au-P repulsive potentials are almost identical, see Fig. S1 in the ESI. † In principle, the Au-P repulsive potential for auorg a 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 a 0 and auorg a .

Computational details
In the original parameterization of the auorg set, 58 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 67 because it was noted by Kepp 68 and Goel et al. 69 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). 70 No dispersion correction was employed 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)
Here we included the empirical D3 dispersion contribution 71,72 in both DFT and DFTB2 calculations. 73 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. 23 In the DFT calculations we employed an effective core potential (ECP) 74 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. 75 All non-periodic DFT calculations were carried out using the implementation of the ORCA code. 76 DFTB2 single point energy and geometry optimization calculations were performed with the DFTB+ code. 77 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.
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 6 (dppp) 4 ] 2+ , [Au 8 (PPh 3 ) 8 ] 2+ and [Au 9 (PPh 3 ) 8 ] 3+ 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 À1 full width at half maximum.
The transferability of the DFTB parameters was evaluated by comparing DFT and DFTB energy landscapes of PH 3 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 xed. Selective geometry relaxation was used by constraining the P atom of the PH 3 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, 89 the kinetic energy cutoff was 450 eV, and kpoint 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 À6 eV for achieving self-consistent eld energies, and 0.005 eVÅ À1 for the maximum force in case of geometry optimizations.
For the demonstration application, a case study on the large Au 108 S 24 (PPh 3 ) 16 metalloid was performed. The initial positions of Au, S, and P atoms were taken from the experimental crystal structure (CSD code: DAFLOO). The PH 3 and PPh 3 ligands were added based on the position of the P atoms in the CIF le. Then, these complexes went through a two-step geometry This journal is © The Royal Society of Chemistry 2020 Chem. Sci., 2020, 11, 13113-13128 | 13117 optimization: rst, only C and H atoms were optimized while Au, S and P atoms were xed, and second, all atoms were optimized. The fully relaxed complexes were used to calculate binding energy, electronic properties and the IR spectrum.

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 3 , PMe 3 and PPh 3 . 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 a is superior among the three sets with a maximum RMSD of only 0.09Å, while auorg a 0 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 a 0 and least for auorg a (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 3 ligand and smallest for the PPh 3 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, z2.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 3 , PMe 3 , PPh 3 with Au 4 cluster are À33.91, À46.44, and À48.69 kcal mol À1 respectively. Since the overbinding can be caused by low virtual orbital energies, it makes sense that the auorg a 0 is mostly affected by it. As for the ranking in performance for averaged and normalized ligand binding energies, we nd the same order as for the geometries: auorg a is the best parameter in terms of binding energy, while auorg a 0 performs the worst.
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 7 (PH 3  Aer closer inspection, we conclude that this large value results from a strongly distorted Au 7 core shape that can only be stabilized when larger ligands such as PMe 3 or PPh 3 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 3 or PMe 3 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 7 (PH 3 ) 7 ] + cluster, lower RMSDs were observed for auorg a 0 and auorg c 0 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 a 0 and auorg c 0 parameters, with energy deviations as high as À18 kcal mol À1 for the smallest PH 3 -ligated clusters (see Fig. 2 and Table S3 in the ESI †). Again we nd generally that the averaged and normalized ligand binding energy deviations become smaller when the larger PMe 3 ligands are used. Again, auorg a has smaller overbinding than the other two parameters for PH 3 ligands, and this tendency turns to underbinding when the larger PMe 3 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 nd the tendency conrmed that auorg a tends towards least overbinding and auorg a 0 towards greatest overbinding, with auorg c 0 somewhere in between but closer to auorg a 0 . Since small ligands themselves consistently tend to increase overbinding, we nd that auorg a performs best for small ligands and small Au clusters, while auorg a 0 or auorg c 0 parameters are better for larger ligands and larger gold clusters.

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 2 , PP 3 , 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 a and auorg a 0 parameters is presented and discussed in the main text. The performance of the auorg c 0 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 a and auorg a 0 , respectively. The same RMSDs of auorg c 0 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 6 (dppp) 4    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. 94 The gure 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 eld effects that are not present in our gas phase theoretical calculations, and thus we will not discuss differences in ligand geometries.
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 a and auorg a 0 , and Fig. S4 in the ESI † for auorg c 0 . All predictions by DFTB tend towards underbinding, with only three minor exceptions. From the smallest [Au 6 (dppp) 4 ] 2+ up to [Au 13 (dppm) 6 ] 5+ , the auorg a 0 binding energy deviations of #5 kcal mol À1 are in very good agreement with the DFT reference. auorg a trends towards strong underbinding in all cases. In the case of clusters with more core Au atoms (n $ 20), auorg a and auorg a 0 have both noticeable underbinding as high as 20 and 15 kcal mol À1 , 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 renement is reduced in all cases, with a maximum absolute deviation of #9 kcal mol À1 . 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 9 (PPh 3 ) 8 ] 3+ and [Au 11 (PMePh 2 ) 10 ] 3+ clusters, see Fig. S5 in the ESI † for more details.

Gold cluster-ligand bond dissociation energy curves for the [Au 8 (PPh 3 ) 8 ] 2+ 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 3 ) m (n ¼ 8, 9, 11, and 101) on titania aer 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 8 (PPh 3 ) 8 ] 2+ 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 DE of this practically dissociated geometry was dened in all methods as 0 kcal mol À1 . 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.
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 À1 (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 z5-10 kcal mol À1 . In all cases, DFTB2/auorg a 0 and DFTB2/auorg c 0 are very similar to each other and outperform DFTB2/auorg a . 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 xed electronic connement radii optimized to describe geometries and binding energies. 95 The underestimation of phosphine ligands binding to gold "surface" will be discussed further below regarding the adsorption of a PH 3 molecule on the gold (111) surface.

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 a 0 parameter, but show the corresponding performance of the other parameters in the ESI. † Fig. 7 shows the TPSS/def2-SVP and DFTB/auorg a 0 calculated energy levels of the frontier molecular orbitals for selected gold clusters. The plots showing frontier orbitals computed with auorg a 0 and auorg c 0 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 shis are appreciable and may reach z0.9 eV (see Table S6 †), with a general trend towards lower energies. Regarding HLGs, DFTB/ auorg a 0 tends to underestimate those for the singly charged   8 ] 3+ (MIVPOX-D 2h ) clusters for the TPSS and DFTB2/auorg a 0 methods. Despite some inevitable quantitative differences in the orbital topology, the amplitudes of the HOMO and LUMO orbitals obtained with the DFTB2/auorg a 0 method are very close to the ones obtained by TPSS/def2-SVP. Noticeable qualitative differences are the LUMOs of Au 7 and Au 8 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 shis of Au 7 and Au 8 . The HOMO and LUMO plots of DFTB2/auorg a and DFTB2/ auorg c 0 are shown in Fig. S8 and S9 in the ESI. † Both of these two sets of molecular orbitals have similar shapes with DFTB2/ auorg a 0 orbitals; with DFTB2/auorg a also overestimating contribution from the gold core while DFTB2/auorg c 0 underestimates the gold and overestimates the contribution from phosphine ligands in Au 7 . 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. The IR spectra of [Au 6 (dppp) 4 ] 2+ (BOTSOS), [Au 8 (PPh 3 ) 8 ] 2+ (OPAUPF), and [Au 9 (PPh 3 ) 8 ] 3+ (MIVPOX-D 2h ) clusters have been studied previously both experimentally and computationally using DFT 19 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 À1 and Au-P modes above 400 cm À1 . Both PBE/def2-SVP and DFTB2/auorg a 0 -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 8 and Au 9 is the presence of an experimentally distinct peak at 398 cm À1 . 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 À1 in the Au 8 and Au 9 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, 19 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 À1 of Au 8 are related to vibrations with weaker intensity assigned to a phenyl rocking vibration that is attached to protruding Au 1 -P 8 . In the case of Au 9 , 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

Chemisorption of PH 3 on the Au (111) surface
To evaluate the transferability of the new DFTB parameters, we explored the energy landscape of PH 3 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 a , DFTB2/ auorg a 0 , DFTB2/auorg c 0 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 3 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 À1 , and both auorg a 0 and auorg c 0 parameters are able to reproduce qualitatively the relative PBE energy landscape. The top site is the most preferred adsorption location of PH 3 adsorption on the Au (111) surface according to PBE, DFTB2/auorg a 0 and DFTB2/auorg c 0 . However, DFTB2/ auorg a predicts the most preferred adsorption location is somewhere in the middle of top, hollow fcc, and hollow hcp positions.
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 DE z 4.6-7.4 kcal mol À1 for auorg a 0 , the deviation DE z 7.6-9.8 kcal mol À1 for auorg c 0 , and largest with DE z 8.9-9.8 for auorg a . 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 shi, the greater the underbinding. This test shows that DFTB2/auorg a 0 is the best option to study phosphine ligands chemisorbed on Au surfaces, and that our parameter has difficulties for the interaction of phosphines with at Au surfaces, as mentioned above in Section 3.3.

3.7
Case study for Au 108 S 24 (PPh 3 ) 16 : 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 108 S 24 (PPh 3 ) 16 . 25 The core of the metalloid consists of a Russian doll motif with an Au 6 inner octahedron enclosed by a second-shell Au 38 octahedron. This core is then capped at its six tips by novel Au 4 S 4 planar rings, which are then connected by 4 more gold atoms, resulting in a (Au 2 S) 24 outer shell motif. Finally, sixteen Au(PPh 3 ) groups complete this metalloid and are bound to each side of the Au 44 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 le, 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 3 ligands. The geometry was then optimized using the DFTB2/auorg a 0 method. Fig. 10 compares the nal DFTB optimized geometry to the initial, experimental crystal structure. DFTB2/auorg a 0 predicts the averaged and normalized ligand binding energies of Au 108 S 24 (PPh 3 ) 16 to be À63.2 kcal mol À1 , 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 À1 . 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  This journal is © The Royal Society of Chemistry 2020 Chem. Sci., 2020, 11, 13113-13128 | 13123 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.
Computationally, the electronic structure of this cluster has been approximately investigated previously at the BP86/TZVPP level of theory, using a simplied model in which the triphenylphosphine ligands were replaced by PH 3 groups, and the entire structure assumed to take T d point group symmetry.
The HOMO-LUMO gap of this simplied Au 108 S 24 (PH 3 ) 16 cluster was found to be 0.68 eV. Fig. 11 shows the DFTB2/ auorg a 0 HOMO and LUMO plots of Au 108 S 24 , Au 108 S 24 (PH 3 ) 16 , and Au 108 S 24 (PPh 3 ) 16 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 108 S 24 core to fully ligated cluster denotes that the phosphine groups stabilize the core, much like in Au 70 S 20 (PPh 3 ) 12 . 26 The quantitative difference in the energy gaps between Au 108 S 24 (PH 3 ) 16 and Au 108 S 24 (PPh 3 ) 16 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 108 S 24 (PH 3 ) 16 and Au 108 S 24 (PPh 3 ) 16 clusters. In addition, both HOMO and LUMO are shied by zÀ0.80 eV when PPh 3 ligands are replaced by PH 3 ligands. In the case of the HLG, employing the simplied model Au 108 S 24 (PH 3 ) 16 yields a similar value to the HLG of the Au 108 S 24 (PPh 3 ) 16 cluster, however, the Fermi level has shied dramatically by À0.8 eV! Our calculations therefore indicate that the use of a simplied model is inadequate and can indeed signicantly affect the electronic properties of nanoclusters.
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 a 0 Hessian calculations to predict the IR spectrum for the  Au 108 S 24 (PPh 3 ) 16 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 4 (SCH 3 ) 4 , [Au 25 (SCH 3 ) 18 ] À clusters, and six Au 18 clusters protected with various types of ligands to the corresponding DFT and experimental spectra in the case of Au 18 (S-c-C 6 H 11 ) 14 . 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 108 S 24 (PPh 3 ) 16 . 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 À1 region which are attributed to PPh 3 distortions. The multiple peaks in the 200-350 cm À1 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 gure are the normal modes related to the novel Au 4 S 4 planar rings: (1) Au-S stretching at 241.1 cm À1 , (2) S-Au-S-Au-S symmetric stretching at 277.2 cm À1 and (3) S-Au-S symmetric stretch at 517.05 cm À1 . The inset of Fig. 13 shows several peaks <150 cm À1 which are mostly attributed to Au 44 core distortion.
Since the novel Au 4 S 4 planar ring motif of Au 108 S 24 (PPh 3 ) 16 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 4 (SCH 3 ) 4 , [Au 25 (SCH 3 ) 4 ] À and [Au 37 (PPh 3 ) 10 (SR) 10 Cl 2 ] + clusters. [31][32][33] The S-Au-S-Au-S asymmetric and symmetric stretching values are found to be 274 cm À1 and 277 cm À1 (see Peak 2 of the inset in Fig. 13 10 Cl 2 ] + have reported their Au-S symmetric stretches at 293 cm À1 and 239 cm À1 , respectively, which are relatively close to our DFTB-calculated Au-S stretch values; 241 cm À1 as shown as Peak 1 in the inset of Fig. 13, and 287 cm À1 as presented in Tables S11 and S12, † for the Au 108 S 24 (PPh 3 ) 16 cluster. Peak 3 in the inset of Fig. 13 shows another peak at 514 cm À1 , which is another characteristic S-Au-S symmetric stretch. One can depict the similarity of the Au-S stretches in the 200-300 cm À1 regions with that of other sulfur-capped Au clusters. 20,23,[31][32][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 108 S 24 (-PPh 3 ) 16 is reliable and provides a useful ngerprint to identify this cluster in future studies.

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 a < auorg c 0 z auorg a 0 . 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 exible 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 a is the best option to study small ligands and small Au clusters. DFTB2/auorg a 0 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 a and DFTB2/auorg a 0 are good options to be used. The performance of auorg c 0 for Au-P interactions is similar to that of auorg a 0 because the effects of shiing the Au 6p orbital energy from À0.02786 to À0.00001 hartree is less signicant than shiing the P 3d orbital energy from 0.12044 to 0.52044 hartree. However, the effect of shiing the Au 6p orbital energy has not been investigated for Au-(C, H, N, O and S) interactions, thus the use of auorg a 0 is more preferable over auorg c 0 . Besides consideration of cluster and ligand sizes, the geometrical environment of the Au-P bond also inuences 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 a 0 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 posttreatment can be done in experiment without inducing concomitant side-effects such as agglomeration. The auorg a 0 parameter set will be publicly distributed via the http:// www.db.org website. To switch from DFTB2/auorg a 0 to DFTB2/auorg a , 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 108 S 24 (PPh 3 ) 16 nanocluster 25 and investigate its molecular and electronic structure. The optimized Cartesian coordinates of Au 108 S 24 (PPh 3 ) 16 are provided in the ESI. † We found that for large ligand-protected gold clusters, using a simplied 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 108 S 24 (PPh 3 ) 16 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 specicity. 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 nonexclusive, 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-publicaccess-plan).

Conflicts of interest
There are no conicts to declare.