Comparison of S-adsorption on (111) and (100) facets of Cu nanoclusters

The favored site of a surface adsorbate, and the reasons for that site preference, are among the most fundamental types of insight into any surface chemical system. In that vein, early studies of sulfur (S) adsorption on and reconstruction of Cu(111) surfaces indicated a particular stability of structural motifs where a S adatom resides on the four-fold-hollow (4fh) site of a planar square Cu₄ unit. This, in turn, suggested an energetic preference for adsorption of S at more highly coordinated 4fh sites versus lower-coordinated 3fh sites on Cu surfaces.¹ More recent density functional theory (DFT) analysis indicated that reconstructions for the S/Cu(111) system can be stabilized by such motifs.² Along this line, a comprehensive integrated experimental and DFT analysis of step edge decoration and reconstruction for S on stepped Cu(111) surfaces consistently indicated a preference for S at 4fh sites. Specifically, (111) micro-faceted steps, which do not present natural 4fh sites, underwent a complex S-induced reconstruction in which Cu atoms shift from their original sites and thereby form a Cu atom base which enables S adsorption at 4fh-type sites.³ 

As an aside, S on other metal(111) surfaces appear to exhibit a similar behavior. S-induced reconstructions on Ni(111) have been observed to incorporate presumed stable Ni₄S units.⁴ Ag–S complexes which form on Ag(111) at low temperature, including Ag₁₆S₁₃ and larger elongated complexes, consist of overlapping units of Ag₁₆S₁₃, also incorporate prominent Ag₄S motifs.⁵ 

The determination and comparison of the binding energies for S on extended (100) and (111) surfaces of Cu is most naturally performed with plane-wave DFT analysis utilizing a slab geometry with periodic boundary conditions. Stabilities of both chemisorbed sulfur atoms and Cu–S complexes have been studied using this method.^6–8 A series of calculations with increasing lateral unit cell size with one adsorbate per unit cell enables estimation of the behavior in the limit of zero coverage (or infinite cell size). Unfortunately, such slab calculations often exhibit a surprisingly strong dependence of energetics on the choice of slab thickness, i.e., they can suffer somewhat slow convergence to a limiting behavior for infinite thickness (corresponding to the semi-infinite surface system of interest). It has been proposed that appropriate averaging of results over slab thicknesses can eliminate such quantum size effects (and also k-points and/or basis sets convergence issues).^8,9 We return to this theme below.

In this contribution, to provide a more extensive analysis of the adsorption site dependence of S bonding than the above type of slab calculations, we consider the behavior for sequences of square pyramidal nanoclusters with exposed base (100) facets, as well as tetrahedral nanoclusters with exposed (111) facets. As an aside, such analysis is potentially also relevant for characterization of chemisorption on supported metal nanoclusters. For sufficiently large clusters above ∼280 atoms, we find a consistently strong preference for binding at 4fh sites on (100) facets versus 3fh sites on (111) facets by ∼0.6 eV. However, highly accurate DFT calculations show variations in binding of around 0.4 eV for clusters as large as 200 Cu atoms. Furthermore, there is no sign of the often-anticipated¹⁰ simple exponential decay in the size dependence of the adsorption energy, even for systems with linear size as large as 3 nm. As a consequence, this brings into question a picture of the S–Cu chemical bond as being local in nature.

The above observations highlight two related challenges in understanding these adsorption systems. As emphasized above, adsorption energetics for clusters of finite size (or for slabs of finite thickness) can exhibit strong deviations from the behavior on extended surfaces. This derives in part from the lack of localization in chemical bonding which in turn complicates the characterization of such bonding, including the understanding of the difference in bonding between 3fh and 4fh adsorption sites. Actually, it has been long recognized, but perhaps under-appreciated, that locality arises from cancellation of different phases of the Bloch states in extended systems.^11–13 This type of cancellation should not be expected to occur for calculations performed on a single cluster with simple geometric shape, even with hundreds of atoms, as coherent interference can occur between electron waves scattering from the different cluster surfaces. Elimination of the strong size dependence and associated enhancement of localization should occur by introducing some type of randomization into the system, e.g., by incorporating random defects, or by introducing rough surfaces. Below, we describe two strategies to mimic such randomization which we propose will reduce the size-dependence of energetics, thus making binding strength and site preference better match those for the extended semi-infinite surface.

Suitably averaging over the energetics of clusters of different sizes is one way to introduce the cancellation effect described above. We find that by averaging results for a range of cluster sizes, N_Cu measured in atoms (roughly speaking in the range from N_Cu = 100 to 400), one can achieve essentially the same adsorption energies using finite clusters as those from slab geometry calculations.

A more efficient method to assess the behavior in extended surfaces is to utilize partial (fractional) occupancies, which are implemented in many DFT codes, to smear out the effect of the Fermi (HOMO) energy. In Sec. IV, we explore the effects of broadening the occupancy function and show that much faster convergence to energetics for the semi-infinite surface system can be achieved by judiciously choosing the broadening parameter. Furthermore, comparing the density of states (DOS’s) projected onto the adsorbate using the broadened occupancy, the role of antibonding orbitals is clarified, thus facilitating understanding of the difference in adsorption energy between the 3fh and 4fh sites.

Section II briefly summarizes the computational methods used in this paper. The main results comparing S binding on (111) and (100) facets of clusters of various sizes and averaging over large cluster sizes are presented in Sec. III. Results obtained by broadening the occupancy function, and the associated understanding of the difference in bonding at 3fh and 4fh sites, are presented in Sec. IV. Further discussion and conclusions are provided in Sec. V.

DFT calculations are performed using both plane-wave (VASP^14,15 version 5.4.1) and Gaussian (NWChem¹⁶) basis sets. More technical details can be found in a previous paper.¹⁷ All calculations are without spin polarization, except for the S₂ dimer in vacuum. PBE¹⁸ functionals are used in VASP and NWChem calculations. The hybrid PBE0 functional¹⁹ is also used in NWChem calculations, and its screened version (HSE06²⁰) is used in VASP calculations. For VASP calculations, the PAW potentials for Cu and S that are optimized for the PBE functional are used.²¹ The cutoff energy for the plane-wave basis set is 280 eV. For NWChem, basis sets are Los Alamos National Laboratory double zeta with effective core potential (LANL2DZ ECP) for Cu²² and 6-311++G(d,p) for S.^23–25 Some results are also checked with the larger basis sets def2-QZVP and def2-QZVPPD.^26,27

Calculations of S adsorption are performed using VASP for both slab and cluster geometries. For slab calculations, surfaces are simulated by periodic slabs of various thicknesses separated by 1.2 nm of vacuum. Supercells are chosen so that two of the basis vectors are that of superlattices of Cu(100) or Cu(111) surface, and the third is perpendicular to the slab surface. For clusters, orthorhombic supercells are used so that each supercell contains one Cu cluster, separated by 1.2 nm of vacuum in all three directions. NWChem calculations are performed for clusters only, but with open boundaries.

All clusters considered in this paper are formed by truncation of bulk fcc Cu. One can regard them as being constructed by starting from a single atom and then adding various numbers of layers with suitable structure and increasing areas. The S atom will be adsorbed near the center of the last largest layer added. Two classes of clusters are thereby constructed. To mimic adsorption on a (111) surface, we add hexagonally close-packed equilateral triangular layers with side lengths 2, 3, up to l (in atoms). This generates a series of clusters of tetrahedral (T_d) symmetry. The number of Cu atoms N_Cu in a cluster can be written as N_Cu = l(l + 1)(l + 2)/6. For l = 3m + k, where m ≥ 0 and 0 ≤ k ≤ 2 are integers, the center of the facet is a fcc site, hcp site, or top site, if k = 0, 1, or 2, respectively. For the 3fh site, we choose the center fcc and hcp sites when k = 0 and 1, and the fcc site closest to the center when k = 2.

To mimic adsorption on the (100) surface, we instead add square layers with side lengths 2, 3, up to l (in atoms). The clusters thus generated can be viewed as octahedral clusters cut in half, thereby denoted as $Oh2$ clusters and N_Cu = l(l + 1)(2l + 1)/6. Only for l = 2m, the center of the top layer is a 4fh site, so for l = 2m + 1 we choose the 4fh site closest to the center. Examples of clusters of both 3fh and 4fh sites are shown as insets in Fig. 1.

Table I shows results of the adsorption energy, E_b determined with different methods and exchange-correlation functionals. The adsorption energy E_b is calculated by E_b = E(S + Cu_n) − E(Cu_n) − E(S₂)/2, where E(S + Cu_n) is the total energy of the Cu cluster with a single S atom adsorbed, E(Cu_n) is the total energy of the Cu cluster itself, and E(S₂) is the energy of a S₂ molecule in vacuum. For VASP calculations, the Gaussian smearing of width 0.02 eV is used. There is no smearing in NWChem calculations. Using the same PBE exchange-correlation functional, the difference between E_b obtained from plane-wave and Gaussian basis sets is generally within 0.10 eV, i.e., there is excellent agreement between the two approaches. This validation process is important, since medium to large size metal clusters are not the natural environment for either plane-wave or atomic basis set DFT codes. The agreement between the two different methods gives confidence that results presented below do not reflect numerical artifacts.

Results using the PBE0 and HSE06 functionals also generally agree well with the PBE results, the difference usually being within 0.1 eV. However, there are certain clusters (e.g., 30-atom $Oh2$⁠, 20-atom T_d) where the difference is significantly larger. Also the consistency of results for PBE0 obtained with different Gaussian basis sets is not as good as for PBE. The largest differences in the Gaussian basis sets show up in the 5-atom $Oh2$ cluster and the 20-atom T_d cluster.

Figure 1 shows the adsorption energy E_b of S on 3fh sites on T_d clusters and 4fh sites on $Oh2$ clusters of various sizes from VASP calculations. Two sets of data are calculated for each geometry. The first set, represented by solid lines in Figure 1, has the Cu atoms in the cluster fixed at their bulk positions, allowing only the S atom to relax. The second set, represented by dotted lines, allows all atoms to relax. Results are obtained using the plane-wave basis set.

The somewhat surprising result in Figure 1 is that not only is there a very large size dependence in E_b, but also the preference for 4fh over 3fh only emerges for very large clusters. For N_Cu < 100, E_b is very sensitive to the cluster size, and the variation with N_Cu dominates over any site preference. Even for N_Cu > 100, E_b can be very close for the two types of adsorption sites for clusters of similar sizes, although the preference towards 4fh sites does emerge as a trend.

Results with the fully relaxed clusters are mostly in line with the counterparts for a fixed substrate. For some of the smaller $Oh2$ clusters, however, larger deviations are observed. This can be explained by the observation that the exposed (100) surface is much less thermodynamically stable and will sometimes reconstruct from the pristine (100) structure. Also for $Oh2$ clusters, sometimes the clean and S-adsorbed clusters can relax into different shapes. For these occasions, we choose the more stable S-adsorbed configuration as the starting point and redo the calculation for the metal cluster with an S atom removed. In most cases, relaxation lowers the value of E_b slightly, although some exceptions can be found for S on 4fh sites of $Oh2$ clusters.

As indicated in Sec. I, by suitably averaging binding energies over a range of (larger) cluster sizes, one might be able to efficiently assess the adsorption behavior on extended surfaces. In general, binding energy displays quasi-periodic variation as a function of linear cluster size, which arises from the interference of the cluster boundaries and the electronic wave functions. Thus, it is natural and appropriate to average over a number of periods in order to extract a limiting large-size behavior. We note that the period depends on the cluster geometry and indeed is different for our analysis of binding at 3fh versus 4fh sites. For 3fh sites, averaging over N_Cu from 84 to 364 which corresponds to roughly two periods of oscillation yields E_b = − 1.78 ± 0.04 eV for unrelaxed substrates and −1.84 ± 0.05 eV for relaxed substrates. The errors are estimated using the standard deviations of the data divided by the number of samples, thus reflecting the general expectation that by averaging a larger range of cluster sizes, one can better approach the limiting behavior. For 4fh sites, averaging over N_Cu from 91 to 385 which corresponds to roughly one period of oscillation yields E_b = − 2.36 ± 0.03 eV for unrelaxed substrates and −2.37 ± 0.03 eV for relaxed substrates. These results are shown in Fig. 1 as horizontal solid lines running through data points that are used for the averaging.

We also calculate independently the S adsorption energy using a periodic slab geometry. For the (100) surface, large oscillations in E_b as a function of the slab thickness are found. These are due to the 2D quantum confinement effect.  Appendix B illustrates these effects through an analysis with (2 × 2) supercells (1/4 ML S coverage). To obtain bulk adsorption energies, we average over DFT results for slab thicknesses from 7 to 12 layers and obtain E_b = − 2.400 ± 0.002 eV with θ_S = 1/16 ML for an unrelaxed substrate and E_b = − 2.468 ± 0.006 eV with θ_S = 1/20 ML for a relaxed substrate. For the (111) surface, less thickness dependence is found, and we calculate the bulk adsorption energy by averaging slab thicknesses from 4 to 7 layers to obtain E_b = − 1.778 ± 0.003 eV with θ_S = 1/12 ML for an unrelaxed substrate and E_b = − 1.926 ± 0.004 eV with θ_S = 1/16 ML for a relaxed substrate. At the right side of Figure 1, we show the calculated E_b for fcc sites on Cu(111) and 4fh sites on Cu(100) with the periodic slab geometry. Consistent with the trend established for large Cu clusters, S adsorption on the 4fh site is stronger than the 3fh site in the slab geometry calculations. Note that with averaging, the cluster results are completely consistent with the slab results for unrelaxed substrates, while some deviations exist for relaxed substrates.

Note that here we focus on 3fh and 4fh sites. For S on extended Cu surfaces, other adsorption sites are significantly less favorable. DFT-PBE calculations show that the adsorption of a sulfur atom on a bridge site is 0.95 eV weaker than the 4fh site on the Cu(100) surface. Adsorption on a top site is even less favorable, being 1.54 eV weaker than the fcc site on the Cu(111) surface. Thus bridge sites and top sites have negligible population.

We conclude this subsection with some remarks about the averaging procedure. In the free electron picture, the quasi-periodic behavior of E_b arises from interference of the wave functions reflected by cluster or slab boundaries. If one can make the linear size l of the system a continuous variable, e.g., using a jellium model, then E_b and other physical quantities can be described as piece-wise continuous curves, with periodicity λ_F/2 for l → ∞,²⁸ where λ_F is the Fermi wavelength. For the averaging procedure to be effective, the phases of the data points on this oscillatory curve should be incoherent, or in other words, more or less evenly distributed among the hills and valleys of the curve. If this condition is satisfied, then the average will not be very sensitive to the range of sizes used and also approach the limiting value rather quickly. We find that this is generally true for the systems studied here. However, there are systems, e.g., (110) surfaces of coinage metals, where the phase incoherence requirement is not met.²⁹ In this case, the averaging procedure is not very effective in eliminating the quantum size effect, even averaging over slabs of up to 12 layers.³⁰ 

Figure 2 shows the bond length between S and its nearest-neighbor Cu atoms from VASP for the same sets of configurations as those in Fig. 1. Unlike the adsorption energy, the respective S–Cu bond lengths for S at the 3fh and 4fh sites converge rather quickly, basically reaching their bulk limits for N_Cu > 100. Furthermore, the bond length for S on 3fh sites is about 3% shorter than on 4fh sites. The convergence to the bulk value, as plotted at the right side of the figure, is also quite apparent. The asymptotic value of 0.229 nm for Cu–S bond length at the 4fh site is slightly larger than the 0.226 nm value obtained from an experimental photoemission study.³¹ This is consistent with the general level of accuracy of DFT/PBE.

It is interesting to note that the bond length predicted by optimization of the S with a fixed substrate using the Gaussian basis sets of LANL2DZ (Cu) and 6-311++G(d,p) (S) is about 3.5% longer than the VASP prediction. This is likely due to the shortcomings of the 6-311++G(d,p) basis set for treating S. Using def2-QZVPPD for S instead predicts bond-lengths which are only 0.5% longer than the VASP values. However, a combination of LANL2DZ and def2-QZVPPD results in an unbalanced description of the system, with a much larger basis set on S than on Cu, which causes overbinding. A combination of def2-QZVP (Cu) and def2-QZVPPD (S) gives good agreement with VASP results for both adsorption energy and bond lengths (see Table I).

As shown in Sec. III, for an isolated cluster, quantum confinement of electrons introduces a correction to the large-size limit of the adsorption energy that does not decay exponentially with the system size. We also find that removing one or more atoms from the corners of a cluster can change the adsorption energy by as much as 0.4 eV for a cluster of about 100 atoms.¹⁷ As mentioned in Sec. I, these features reflect a lack of locality of chemical bonding in metallic solids. In our case, the clusters consist of a few flat surfaces (together with some edges and corners), which can create coherent interference in the wave functions. Again, localization and thus minimization of size effects come from cancellation of the phase of Bloch waves which can be produced by introducing randomness into the system. Our proposal here is that by introducing such effects to reduce size dependence, we can more efficiently assess the energetics of the semi-infinite extended surface system. Further validation of this idea is provided below.

Specifically, in this section, we explore the technique of partial (or fractional) occupancy that has been implemented in many DFT codes as a way to introduce the above-mentioned phase cancellation. In real solid systems, the probability of occupancy of energy levels for electrons approaches that of a step function, but it is often more efficient numerically in solid state electronic calculations to broaden the step function (or, more exactly, the Fermi-Dirac distribution).³² The key physics is that the position of the Fermi level, relative to the electronic band structure, is sensitive to the system size.²⁸ By adding noise to the exact position of the Fermi level, one can simulate randomness in a system. The smearing method, by broadening the occupancy function, adds uncertainty to the Fermi level and is thus a natural way to simulate “noisy” Fermi levels.

Figure 3 shows E_b calculated for unrelaxed metal substrates with Gaussian smearing but deliberately choosing a larger smearing width σ than the default value 0.2 eV used in Fig. 1. The size dependence is greatly reduced, and the convergence to the limiting large-size value of E_b = − 1.78 (−2.39) eV for 3fh (4fh) sites is more apparent. The larger the σ values, the smaller the extent of size dependence. The dramatic reduction in size dependence is consistent with the above stated proposal that enhanced smearing mimics the introduction of randomization to the system which in turn enhances localization. Ideally, the more readily assessed limiting large-size behavior evident from this analysis provides an efficient assessment of binding on a semi-infinite extended surface.

One caveat is that with large σ, the detailed form of the smearing becomes relevant. Using the Methfessel-Paxton (MP) scheme,³² for which the occupancy function approaches a step function faster than for Gaussian smearing as σ decreases, leads to somewhat different results for large σ. For example, using the first-order MP with σ = 1.0 eV, E_b on 4fh sites in $Oh2$ clusters converges to −2.25 eV versus the −2.37 to −2.40 values obtained using the other three methods (averaging different cluster sizes, slab geometries, and Gaussian smearing with σ = 1.0 eV). For 3fh sites, the MP smearing with σ = 1.0 eV yields E_b = − 1.63 eV, versus the −1.78 eV value obtained using the other methods. We conclude that Gaussian smearing is more appropriate for our purposes here.

Strictly speaking, even with Gaussian smearing, different σ values will lead to a different limiting behavior, and the physically relevant value should correspond to the limit of σ → 0. With slab geometries and a relatively small (2 × 2) supercell, we find that between σ = 0.2 and 1.0 eV, the values of E_b for S/Cu(100) do deviate, but the differences are relatively small (about 0.025 eV). For S/Cu(111), on the other hand, the change due to σ is within numerical uncertainties. The optimal choice of the form and width of the smearing function is an open question at this stage.

Perhaps more important than potentially providing a more efficient method to estimate E_b for S on extended Cu surfaces from cluster calculations, we can also use the smearing of the occupancy function to elucidate the difference between the bonding of S in 3fh and 4fh sites. One way to visualize interactions between S and a cluster is through plotting the site-projected density of states (SDOS’s) of individual atoms. Figure 4 shows the SDOS localized on the S on the center 4fh site of a 91-atom $Oh2$ cluster, obtained using Gaussian smearing of different widths σ. With a small σ, the SDOS consists of many sharp spikes, each of which corresponds to one or more molecular orbitals. (As an aside, analogous sharp spikes appear in the SDOS for slab calculations.) The highly complex SDOS, especially near the Fermi level, is directly responsible for the large size dependence of binding seen in Sec. III A. It also makes it more difficult to obtain an intuitive picture of chemical bonding. By widening the smearing, a smoother SDOS can be achieved, which facilitates interpretation of bonding. It is significant to note that Feibelman⁶ also used Gaussian-smearing of the DOS to obtain insights into Cu–S clusters on Cu(111) surfaces. In his case, the DOS was projected onto Cu atoms and his analysis used slab (rather than cluster) geometries.

The solid line in Fig. 5 shows the SDOS of a S atom on the 3fh site of a T_d cluster with 56 Cu atoms, with Gaussian smearing of 1.0 eV. Analysis of the electronic structures using the crystal orbital Hamilton population (COHP) method³³ shows that the peaks near −17 and −8 eV are mostly bonding, and the peak near −5 eV is mostly antibonding. (Note that in Fig. 5 the energy is relative to the reference configurations of individual atoms, rather than the Fermi energy as is the usual practice in solid state physics as in Fig. 4. This is done in order to make the comparison between S on different adsorption sites more transparent.) The dashed line is for an S atom on a 4fh site of the (100) face of a $Oh2$ cluster with 91 Cu atoms. Compared with S on the 3fh site, the main difference in the SDOS is that the antibonding states are more spread out. This results in a higher Fermi energy, E_F, which in turn forces the bonding state deeper below the Fermi level, thus increasing the strength of binding. Thus the difference between S adsorption on the 4fh site and the 3fh site can be understood intuitively in the following way: on a 4fh site, with more neighboring Cu atoms, the S does not have to be as close to the Cu atoms as on the 3fh to maximize the bonding coupling between the S and Cu orbitals. This in turn leads to much smaller antibonding coupling between the S and Cu atoms, which is due to the faster decay of the antibonding interactions as the separation increases. Note that the linear sizes l for the two types of clusters shown in Fig. 5 are the same, and there can be less perfect matches when choosing different clusters. Nevertheless, the qualitative picture remains the same.

Good agreement has been achieved between DFT codes employing plane-wave and Gaussian basis sets, regarding the adsorption of S on Cu clusters of various sizes. However, we find that the large size-dependence in the adsorption energies makes it challenging to estimate the limiting value of binding on an extended surface, and the associated delocalization makes it challenging to elucidate the nature of chemical bonds between the S adsorbate and the metal cluster. It has been long recognized that for small clusters (less than 50 atoms), the discreteness of the orbitals, especially the HOMO-LUMO gap, will lead to a behavior quite different from their bulk counterpart. Another issue, which is familiar in condensed matter physics, is that for an isolated cluster, interference of wave functions from the boundaries will lead to corrections that do not decay exponentially. For Cu clusters, the adsorption energy can be significantly affected (up to 0.6 eV) by what happens 1.5 nm away from the adsorption site.

A natural question is then, how can calculations on small to medium size clusters be relevant to adsorption on extended single-crystal surfaces? A simple but effective method is to average over results for clusters over a suitable range of sizes (as described in Sec. III). One could anticipate similar results from suitably averaging over different shapes, or by performing analysis for clusters with rough side surfaces. Another strategy which is particularly efficient for plane-wave methods is to utilize the partial occupancy technique which was originally developed for numerical efficiency. By choosing an appropriate smearing function (e.g., Gaussian), we can reliably assess binding on extended surfaces from calculations on medium size clusters.

By averaging contributions from different orbitals, we can understand the adsorption of S on metal clusters in a way that is both intuitive and also rests on firm quantitative grounds. We suggest that the stronger binding of S to 4fh sites is due to the weaker antibonding interactions compared with 3fh, while having similar bonding interactions. This interpretation of chemical bonds as a competition between bonding and antibonding interactions through interference energies, as advocated a long time ago by Ruedenberg,³⁴ is key to understanding the site preference of simple adsorbates.

J.S.B., J.L., T.L.W., J.W.E., and D.J.L. performed the computations in this work and were supported for this work by the U.S. Department of Energy (USDOE), Office of Basic Energy Sciences, Division of Chemical Sciences, Geosciences, and Biosciences through the Ames Laboratory Chemical Physics program. Use of NERSC computational resources is acknowledged. PAT provided experimental insights into key issues for S adsorption systems and was supported for this work by NSF Grant No. CHE-1507223. The work was performed at Ames Laboratory which is operated for the USDOE by Iowa State University under Contract No. DE-AC02-07CH11358.

On fcc(111) surfaces, the so-called A-step creates microfacets resembling the (100) surface locally. Thus adsorption of S along an A-step may be akin to adsorption on a 4fh site. In order to study this via the cluster approach, we create steps on top of a cluster by adding an incomplete layer, or an island, on one face of the cluster. In Fig. 6, we consider two types of A-steps, one formed by an island that has its boundary as close as possible to the edge of the cluster, thus exposing a step edge with length l − 2 on a cluster with side length l. Note that the larger island with side length l − 1 consists of Cu atoms on hcp sites, rather than fcc sites. DFT-PBE results for S adsorption along this kind of step edge are shown in Fig. 6 as the black pluses. The average result for clusters with l = 8–12 is −2.52 eV, which is slightly lower than the equivalent value of −2.36 eV for the 4fh site on the (100) surface (Sec. III and Fig. 1). The other type of step has one row of Cu atoms removed from the island in the first type, thus one of the step edges is further removed from the edge of the cluster. See insets of Fig. 6 for illustrations. Results for S adsorption on these types of steps are shown in Fig. 6 as red asterisks. The average value for l = 8–12 is −2.09 eV, which lies between −1.77 eV (3fh) and −2.36 eV (4fh) obtained in Sec. III. Therefore, the expectation that A-steps on Cu(111) are more favorable adsorption sites than flat terraces are met, although some differences are found depending on configurations further away from the step edges.

Here, we quantify how the S adsorption energy depends on the thickness of the slab in calculations with semi-infinite slab geometries. Table II lists the adsorption energy E_b for S on Cu(100) and Cu(111) surfaces, calculated using slabs of different thicknesses measured by the number of layers L. All atoms are allowed to relax except for the bottom layer of Cu atoms. Also listed are the average value 〈E_b〉 and the standard deviation δE_b for each L calculated using data up to L. For example, for L = 12, we use data from 7 to 12. While the extent of variations using slabs is much smaller than results using clusters, the convergence to the bulk limit is also slow. Also note that variations of a few meV can be due to numerical errors.