Stability of Single Gold Atoms on Defective and Doped Diamond Surfaces

Polycrystalline boron-doped diamond (BDD) is widely used as a working electrode material in electrochemistry, and its properties, such as its stability, make it an appealing support material for nanostructures in electrocatalytic applications. Recent experiments have shown that electrodeposition can lead to the creation of stable small nanoclusters and even single gold adatoms on the BDD surfaces. We investigate the adsorption energy and kinetic stability of single gold atoms adsorbed onto an atomistic model of BDD surfaces by using density functional theory. The surface model is constructed using hybrid quantum mechanics/molecular mechanics embedding techniques and is based on an oxygen-terminated diamond (110) surface. We use the hybrid quantum mechanics/molecular mechanics method to assess the ability of different density functional approximations to predict the adsorption structure, energy, and barrier for diffusion on pristine and defective surfaces. We find that surface defects (vacancies and surface dopants) strongly anchor adatoms on vacancy sites. We further investigated the thermal stability of gold adatoms, which reveals high barriers associated with lateral diffusion away from the vacancy site. The result provides an explanation for the high stability of experimentally imaged single gold adatoms on BDD and a starting point to investigate the early stages of nucleation during metal surface deposition.

It is important to ensure that the size of the QM region embedded within the MM region is large enough to avoid any finite-size effects. The appropriate QM region size of the substrate was chosen by comparing various properties of PBE +TS /REBO-optimized embedded cluster models with varying QM region size against the parent PBE +TS -optimized periodic model. All convergence tests were conducted on an idealized surface model. First, the structural deviations of each PBE +TS /REBO-optimized embedded cluster were compared against the initially cut cluster from the PBE +TS -optimized periodic model, which can therefore be taken to be an appropriate representation of the periodic model. As can be seen from Following the structural comparison, the electronic structures of the PBE +TS /REBOoptimized embedded cluster models were compared against the PBE +TS periodic model. Figure S1(b) shows a graph of the band gap, which is the energy between the highest occupied molecular orbital and the lowest unoccupied molecular orbital, of the embedded clusters as a function of QM region size. The band gap generally decreases as the QM region size increases and tends towards the QM periodic value (2.3 eV). This shows that clusters with smaller QM region sizes do experience some finite-size effects. In contrast, as the QM region size increases, the embedded cluster becomes more structurally and energetically similar to the periodic QM structure. Despite this overall trend, over the range of QM region sizes explored, the 30-atom QM region was found to possess the closest band gap to the periodic value, while the 90-atom QM region had the next closest value (2.6 eV). The band gaps of the 10-, 20-, 60-and 70-atom QM regions, which were found to be structurally similar to the periodic model, were found to be at least 0.5 eV larger than the QM periodic value.
S-3 Figure S1. Scatter graphs showing the (a) root-mean-square deviations and (b) band gaps of a single gold atom atop PBE +TS /REBO-optimized cluster models against the initial cluster cut from the PBE +TS -optimized periodic model; and (c) adsorption energies of a single gold atom atop PBE +TS /REBO-optimized cluster models against the PBE +TS -optimized periodic model, all as a function of QM region size.
Finally, the adsorption energy of a single gold atom was evaluated as a function of the QM region size, as shown in Figure S1(c). The 70-and 90-atom QM regions (71-and 91-atoms respectively including the gold atom) resulted in the closest adsorption energy (−0.30 eV) to the QM periodic value (−0.31 eV). The 20-and 30-atom QM regions significantly overestimate the adsorption energy, while the 10-, 40-, 50-, 60-and 80-atom QM regions report similar adsorption energetics to the QM periodic model, but are not as close as the 70-and 90-atom QM regions.

S-4
Taking all results into consideration, the 90-atom QM region results in an optimized structure, band gap and gold adsorption energy most similar to the QM periodic model, and was thus chosen as the optimal QM region size within the QM/MM cluster. While a larger QM region would most likely result in a cluster with a final geometry and energetics more similar to the periodic model, convergence problems were encountered with larger QM region sizes (100 and 110 atoms). Regardless, the embedded cluster with a 90-atom QM region was found to have similar structural and energetic properties to the periodic QM cluster and was thus deemed an appropriate size.
S-5 2 Scaling of QM/MM versus Periodic QM To showcase the higher computational efficiency of the hybrid QM/MM approach, scaling graphs were constructed after conducting single-point calculations on the PBE +TS -optimized periodic model and the PBE +TS /REBO-optimized model, with a 90-atom QM region, after a gold atom was adsorbed onto the model surfaces. Figure S2 shows the computational cost of these single-point calculations as a function of number of cores, which were run on Lenovo NeXtScale nx360 M5 servers with dual Intel Xeon E5-2680 v4 (Broadwell) 14-core processors at 2.4 GHz, as available within the Orac high performance computing cluster provided by the Scientific Computing Research Technology Platform of the University of Warwick. All calculations used the Eigenvalue SoLvers for Petaflop-Applications S1 library and the ELectronic Structure Infrastructure. S2 Figure S2. Scaling graphs of single-point PBE +TS /REBO and PBE +TS calculations of the idealized oxygen-terminated diamond (110) surface. The embedded cluster model comprised 527 atoms with a 90-atom QM region, while the periodic model comprised a 92-atom unit cell.
As can be seen in Figure S2, QM/MM calculations are vastly cheaper than the periodic QM calculations. Furthermore, periodic QM calculations failed when using 16 cores or fewer due to memory issues. This confirms the superior computational efficiency of the hybrid S-6 QM/MM approach and that it can be used to access more computationally-costly methods such as MGGAs and HGGAs.

Benchmarking MM Forcefields
It is important to ensure that the adsorption energetics of a single gold atom are not significantly affected by the choice of the embedding forcefield environment. To investigate further, the adsorption energy of a single gold atom on an idealized surface, as calculated using PBE +TS /REBO, was benchmarked against the Tersoff S3 forcefield, where PBE +TS was used as the complementary QM method. Similar to REBO, the Tersoff forcefield was developed specifically for carbon, with applications to amorphous carbon, S3 and is thus an appropriate forcefield to benchmark the REBO forcefield against. Table S1 benchmarks the REBO forcefield against the Tersoff forcefield. Both the REBO and Tersoff forcefields result in the same adsorption energy for a single gold atom on an idealized surface. Furthermore, there is a very small disparity in adsorption height (0.05Å), showing that both forcefield methods predict virtually identical adsorption energetics for the single gold atom, and that the REBO forcefield is appropriate to embed the QM region within. Table S1. Adsorption energetics and heights of a single gold atom adsorbed onto an idealized oxygen-terminated diamond (110) surface after a PBE +TS /MM optimization, using various MM forcefield methods.

MM Forcefield
Adsorption Height (Å) Adsorption Energy (eV) It is important to ensure long-range dispersion effects such as van der Waals (vdW) forces are treated appropriately, as they can have a significant effect on the adsorption structure.
The pairwise additive TS scheme does not explicitly account for beyond-pairwise vdW interactions. For this reason, the TS scheme was benchmarked against some a posteriori manybody dispersion (MBD) S6,S7 correction schemes, namely the range-separated self-consistently screened (MBD@rsSCS) S8  Dispersion corrections were benchmarked on the idealized, SCOV-defective, and delocalized triel-doped systems. The idealized and SCOV-defective systems were chosen as they would permit dispersion corrections to be benchmarked on both 'more physisorbed' and 'more chemisorbed' systems, respectively. The delocalized triel-doped system was chosen over the localized boron-doped models for three reasons: firstly, with common boron dopant densities, the probability of the dopant atom being within the bulk material is much higher than it being in the top surface layers. Secondly, the delocalized model is applicable to any triel dopant and thirdly, the predicted adsorption height and energy do not differ significantly from the localized case with the boron dopant in the third layer, as shown in the main text. Table S2 details the adsorption heights and energies of a single gold atom after TS-, S-8 Table S2. Adsorption energetics and heights of a single gold atom adsorbed onto various oxygen-terminated diamond (110) surfaces after a dispersion-corrected PBE/REBO optimization, using various a posteriori dispersion correction schemes. No data were attained using MBD@rsSCS for the SCOV-defective system due to the negative polarizabilities for some atoms after the initial FHI-aims calculation settings.  S11-S14 and highlights the importance of including a dispersion correction with such DFAs.
For the SCOV-defective surface, the TS method once again performs quite well with respect to MBD-NL. After a PBE +MBD-NL /REBO optimization, the gold atom adsorbs 0.04Å higher than with TS, and this small disparity is reflected in the adsorption energy, which is only 0.02 eV weaker than with TS. No data were attained using MBD@rsSCS for this system, due to the negative polarizabilities for some atoms after the initial FHI-aims calculation settings, which prevented the MBD@rsSCS calculation from completing. It should be noted that this is a technical limitation of the MBD@rsSCS approach that can occur in some systems under certain conditions, and is not physically meaningful. Without any dispersion correction, the gold atom again adsorbs higher than dispersion-corrected approaches, and the adsorption energy was evaluated to be at least 0.25 eV weaker, further attesting to the need for a dispersion correction.
Finally, for the delocalized triel-doped system, there does appear to be some dependency on the choice of dispersion correction. The TS correction predicts stronger adsorption than both MBD approaches, and the gold atom optimizes to a site far closer to the surface. The MBD approaches perform very similar to each other, with the gold atom predicted to adsorb 1.06Å above the surface with an adsorption energy just larger than −1.7 eV. The consistency between the MBD approaches and the relatively weaker adsorption energies are indicative of the beyond-pairwise interactions being taken into account, and these effects have a greater influence within this charged system rather than the neutral idealized and SCOV-defective systems. The lack of a dispersion correction yet again resulted in the gold atom adsorbing higher than dispersion-corrected approaches, and the adsorption energy was evaluated to be at least 0.18 eV weaker.
S-10 Figure S3 shows the binding energy curves using various dispersion-corrected PBE/REBO calculations on the idealized, SCOV-defective, and delocalized triel-doped surfaces.  Table S2 and in literature. S11-S14 Figure S3 Finally, all the curves for the delocalized triel-doped surface have a minimum at a height of 2.0Å, as can be seen in Figure S3 S-12

Conformational Isomers of SCOV-Defective Surfaces
To ensure the SCOV defect was correctly modeled with every DFA, DFA i , the carbonyl oxygen was first removed and the surface was reoptimized using DFA i /REBO. After this initial optimization, the surface structure at the defect site, centred at the former carbonyl carbon atom, should change from bent (originally trigonal planar with the carbonyl oxygen atom in the idealized system) to trigonal pyramidal. Because diamond surfaces are usually hydrogen-terminated after chemical vapor deposition growth, S16 uncoordinated carbon atoms were subsequently saturated with hydrogen species and the surface was reoptimized using DFA i /REBO, after which the surface structure at the defect site should change to tetrahedral. Based on valence shell electron pair repulsion theory, this shows a return of the ·C· atom to an sp 3 -hybridized state. This is the correct surface configuration as an oxygen atom is needed to pull the carbon atom above the diamond (110) surface plane to form a carbonyl group at the surface. S17 Without this oxygen atom, the carbon atom would remain in an sp 3 -hybridized configuration.
Only the DFAs that returned the former-carbonyl carbon atom to an sp 3 -hybridized configuration were investigated further. This was evaluated by studying the conformational isomerism of the structure centered at the former-carbonyl carbon atom. After the removal of the carbonyl oxygen atom and optimization with a given DFA, the dihedral angle between a surface ether oxygen atom and a surface carbon atom, along the bond between the former-carbonyl and corresponding ether carbon atoms, was calculated. minimum and as explained above, is not the correct physical conformation for the surface after the removal of a carbonyl oxygen atom. In order to ensure lower-rung DFAs could still be benchmarked against HGGAs, the final PBE +TS /REBO-optimized SCOV-defective structures were reoptimized using the respective HGGA +TS /REBO method. Table S4 details the difference between the energy of the synclinal conformation and the energy of the anticlinal conformation, as calculated using various TS-corrected HGGAs. As can be seen in Table S4, the energy difference between the two conformations is between 0.73-0.86 eV, indicating the greater stability of the synclinal conformation. For clarity, the Newman projections S18 of the synclinal and anticlinal conformations are also provided in Figure S4. Table S3. Dihedral angles between a surface ether oxygen atom and a surface carbon atom, along the bond between the former-carbonyl and corresponding ether carbon atoms after the removal of a carbonyl oxygen atom from the idealized oxygen-terminated diamond (110) surface and subsequent optimization with various DFAs. DFAs are ordered from low-to high-rung, along Jacob's ladder, S32 and the TS S10 dispersion correction method was applied to all GGAs and HGGAs.  Figure S4. Newman projections S18 of the synclinal (left) and anticlinal (right) conformations after optimization with various DFAs. The projection is along a bond between the formercarbonyl carbon atom and a surface ether carbon atom. The synclinal conformation is the correct model for the SCOV defect (prior to hydrogen saturation).

S-15
6 Binding Energy Curves for SCOV-Defective Surfaces Figure S5 shows the binding energy curves for the SCOV-defective surface, as calculated with various DFAs. As can be seen, LDAs and most GGAs predict similar binding energy curves.
The revPBE and PZ-LDA DFAs predict the strongest adsorption at an adsorption height of 3.0Å. PBE and PBEsol have similar curves to each other. The RPBE GGA predicts the weakest adsorption among GGAs and has a binding energy minimum at 4.0Å, which is a larger adsorption height value than all other LDAs and GGAs, much like in the idealized case. In contrast, all MGGAs result in very shallow binding energy curves for the single adatom. Furthermore, the binding energy curve for TPSS, much like RPBE, has a minimum at a value larger than all other DFAs. However, the revTPSS binding energy curve is very similar to other MGGAs' despite the strong adsorption predicted in Figure 6. Much like with the MGGAs, the PBEsol0 HGGA also has a shallow binding energy curve, while the other HGGA binding energy curves are very similar to the PBE binding energy curve. Figure S5. Unrelaxed binding energy curves showing the adsorption energy of a single gold adatom as a function of height above an oxygen-terminated diamond (110) surface with a saturated carbonyl oxygen vacancy defect. Density-functional approximations are divided according to (from left to right): local-density approximations (LDAs), Tkatchenko-Scheffler (TS)-corrected generalized gradient approximations (GGAs), meta-GGAs (MGGAs), and TS-corrected hybrid GGAs (HGGAs). S-16