Molecular adsorption on silicon (001): A systematic evaluation of size effects in slab and cluster models

First-principles calculations are in wide use today to describe chemical processes occurring on the silicon (001) surface. The number of atoms that can be explicitly treated is limited and hence size-constraints are invariably required; this applies to both cluster and periodic slab approaches. Using a trial set of seven molecular adsorbate conﬁgurations, we examine the dependence of calculated adsorption energies on several size parameters, namely thickness and in-plane unit cell size for slab models, as well as thickness, length, and width for cluster models. Size-converged adsorption energies are estimated by extrapolation, and are used to assess the accuracy of the more typically-sized slabs and clusters in common use today. Use of a DFT method that can be applied to both slabs and clusters allows us to assess the performance of these two approaches on an equal footing. Copyright 2013 Author(s). All arti-cle content, except where otherwise noted, is licensed under a Creative Commons Attribution 3.0 Unported License . [http://dx.doi.org/10.1063/1.4802837]


I. INTRODUCTION
The silicon (001) surface is a popular target for experimental and theoretical investigations because of its importance to semiconductor technology. As electronic devices become relentlessly smaller there is an increasing need to fully understand and then exploit the chemistry occurring on this surface. This is perhaps best illustrated by the recent demonstration of a single-atom transistor 1 which was fabricated using sophisticated lithographic techniques 2 and a detailed understanding 3,4 of the phosphine on Si(001) chemistry. Computational studies of molecular adsorption and dissociation have an important role in developing this understanding, helping to rationalise experimental observations or make predictions to guide future experiments (see e.g. topical reviews in Refs. [5][6][7][8][9][10]. Needless to say, the utility of such studies is directly related to the explanatory and predictive accuracy of the computational model used, which in turn is set by a number of methodological choices made. The objective of this present work is to examine one of the key parameters that determines the accuracy, namely the size and type of the silicon surface representation. In this we focus on calculated adsorption energies in the low-coverage/isolated-adsorbate limit as this is the typical quantity on which the construction of potential energy diagrams and reaction pathways is based. Silicon (001) is modelled using one of two principal approaches. Cluster models describe the atoms in the immediate vicinity of the adsorption site with all atoms further away being neglected. Slab models, in contrast, describe an infinitely extended surface by imposing a periodic repeat on a finite number of explicitly treated atoms. Both clusters and slabs require choices to be made in regard to the number of atoms to be included in the calculation, in order to minimise any effects due to, respectively, the cluster termination and a periodically repeated adsorbate. choices that affect accuracy, such as the quality of the exchange-correlation treatment and the basis set/pseudopotential. Size effects in cluster and slab models representing Si(001) have been explored in a number of studies, [11][12][13][14][15][16][17][18][19][20] often in the context of a specific adsorption or reconstruction problem. Direct comparisons between cluster and slab approaches are somewhat less common and tend to focus on combining the individual strengths of cluster and slab approaches (such as, more accessible high-level exchange-correlation methods in clusters and a more "natural" representation of the in-plane periodicity of the clean surface in slabs; see e.g. Refs. 15, 19, and 20).
We will be addressing here two perhaps more immediate questions, namely, to what extent commonly used cluster and slab sizes approximate a size-converged limit, and whether the two approaches even converge to the same limit. As illustrated in Ref. 19, these questions are not as trivial as they may appear at first glance. While a reasonable degree of agreement was achieved, direct comparisons between cluster and slab approaches are typically made difficult by that fact that different methods are used to solve the density functional theory (DFT) equations (Gaussiantype orbitals for clusters and planewave pseudopotentials for slabs; see Ref. 19). We avoid this problem here by using a DFT formalism that can be equivalently applied to both periodic slabs and non-periodic clusters. This ensures that calculated adsorption energy differences are directly related to the particular slab/cluster surface representation used, which affords for the first time a much needed direct comparison between Si(001) cluster and slab approaches on entirely equal footing. The size-dependencies in turn are likely to be representative for other electronic structure methods. We acknowledge in passing that such comparisons have been reported for other materials, including germanium, 21 magnesium oxide, 22 and zeolites. 23,24 Our discussion is organised as follows. In Section II, we set out our computational methodology, including a detailed description of how our slab and cluster models are constructed. In Section III A, we examine the performance of slab models when calculating adsorption energies for a trial set of seven adsorbates (see Fig. 1): acetylene (C 2 H 2 ), ethylene (C 2 H 4 ), benzene (C 6 H 6 ), formaldehyde (CH 2 O), hydrogen (H 2 ), water (H 2 O), and ammonia (NH 3 ). Specifically, we look at the effects of increasing slab thickness from 4 to 20 atomic layers, and extending in-plane unit cell sizes from (2×2) to (8×8) and (2×16). In Section III B, we turn to cluster models, and similarly observe the effects of cluster thickness (up to 20 layers), length (up to 25 dimers), and width (up to three rows) on calculated adsorption energies. Throughout Section III, we discuss the technique of compositing that allows us to economically estimate size-converged energies by combining results from several smaller calculations. The size-converged energy estimates in turn are used to assess the performance of smaller slab/cluster models that are typically used in the literature. Section IV provides a summary of our results, and suggests how they may be applied.

A. Density Functional Theory
All cluster and slab calculations reported here are performed using DFT as implemented in the DMol 3 software. 25 The electronic Kohn-Sham eigenfunctions are expanded in terms of a doublenumerical plus polarization (DNP) basis set of atom-centered functions. Exchange and correlation are treated in the generalized gradient approximation (GGA) as defined by Perdew, Burke, and Ernzerhof (PBE; Ref. 26). Orbital occupations are subject to a thermal smearing of 0.1 eV. All energies reported include a correction for scalar-relativistic effects. 27 The cubic unit cell of bulk silicon is geometry optimised using an 8×8×8 Monkhorst-Pack k-point grid. 28 This results in a calculated lattice constant of 5.473 Å, which is in good agreement with other GGA-DFT results. The small overestimate of the experimental value of 5.43 Å (Ref. 29) is typical for GGA-DFT calculations.

B. Slab models
Slab models of the Si(001) surface of n layers thickness are created by combining n atomic layers of silicon at bulk positions with an unrelaxed dihydride termination on the bottom slab surface, i.e. the side of the slab that is not used for adsorption. Figure 2(a) provides a side view illustration of a 4-layer slab with a molecular adsorbate on one side and a dihydride termination on the other. The dihydride termination is constructed by placing hydrogen atoms along the broken bulk Si-Si bonds such that the Si-H bond length is 1.49 Å. 30 During all geometry optimisations, the bottom silicon atomic layer as well as the dihydride hydrogen atoms are held fixed, whereas the n−1 top-most silicon layers and any adsorbates are allowed to fully relax. The periodic unit cell vectors are held fixed at the corresponding bulk values for the in-plane directions, and at 105.8 Å in the plane-perpendicular direction. With these settings, our slabs are separated by a vacuum gap of at least 75 Å from their periodic images. A separate dipole correction is not applied. In addition to the thickness, the size of the slab models is defined by the in-plane periodic repeat. A (2×1) unit represents a single Si-Si dimer. Larger supercells are defined using conventional notation. For example, (4×4) represents eight dimers in two rows and (6×6) represents eighteen dimers in three rows. Any exposed surface dimers are buckled into a p(2×2) pattern. Figure 2(b) shows some of the supercells used.
The irreducible Brillouin zone of our slab models is sampled using a 4×4×1 Monkhorst-Pack k-point grid 28 for a (2×2) surface unit. In trial calculations we confirm that this grid is sufficiently dense for the purpose of this work. For example, the use of a larger 6×6×1 grid resulted in adsorption energy changes of 5 meV or less. For larger supercells, we use proportionally smaller grids, rounding up to the nearest integer in the case of non-commensurate surface units such as (6×6).

C. Cluster models
In order to define our cluster models we require first a suitable surface template, to serve as a prototype atomic representation of the Si(001) surface. Briefly, our clusters are constructed by selecting a subset of silicon atoms from this template and placing hydrogen atoms along any Si-Si bonds that are broken when the cluster is extracted. In fact, defining the precise positions of these cluster-terminating hydrogen atoms is the primary purpose of the template as these atoms will remain fixed in position during geometry optimisations. This affords a simple and effective representation FIG. 2. Overview of slab models used in this work, with C 2 H 2 as an example adsorbate. (a) Side view of the 4-layer surface, with the terminating dihydride layer at the base of the slab. Buckled Si-Si dimer atoms at the surface are counted as a single atomic layer. (b) Top view of the (8 × 8) unit cell, showing only the adsorbed molecule and the two surface-nearest silicon layers. The smaller (2 × 2), (4 × 4) and (6 × 6) unit cells are indicated with red squares. Silicon, carbon, and hydrogen atoms are coloured black, grey, and white, respectively. of the strain experienced by the cluster due to the surrounding silicon atoms that are not explicitly treated. Hence, some care should be taken to define this template.
In the present work, our surface template is based on a fully optimised slab of a (2×1) H:Si(001) monohydride surface. The H:Si(001) slab has a thickness of 40 atomic layers of silicon and is symmetrically covered on both sides with a monolayer of hydrogen atoms. Repeated slabs are separated by a vacuum gap of approximately 50 Å. The slab is geometry optimised using a doubled, (2×2) in-plane unit, with k-point and fixed lattice vector settings as described above for slabs. The silicon atom positions thus obtained constitute our Si(001) surface template; the positions of the monohydride hydrogen atoms are discarded, being of no further utility. The purpose of using a monohydride slab is to ensure that all surface Si-Si dimers remain unbuckled/flat in the template. Buckled dimers are energetically preferred in both clusters and slabs; however, we do not want our cluster-termination to predispose a cluster to one dimer-buckling orientation over another. A template with symmetric, flat dimers satisfies this criterion.
With our surface template defined, cluster models are constructed by selecting silicon atoms from the template as outlined above. Cluster sizes are conveniently defined in terms of thickness (number of silicon layers represented), length (number of Si-Si dimers per dimer row), and width FIG. 3. Illustration of some of the cluster models used in this work. The notation (iD,jR,kL) denotes cluster length, width, and thickness in terms of number of dimers per row (D), number of rows (R), and number of atomic layers (L), respectively. Panels (a-i) show single row clusters of increasing thickness, while panels (j) and (k) show three-row clusters of two thicknesses. Silicon atoms and cluster-terminating hydrogen atoms are colored black and white respectively.
(number of dimer rows). Cluster-terminating hydrogen atoms are positioned along all Si-Si bonds for which only one of the silicon atoms is included in the cluster. Using the silicon atom positions from the template, a terminating hydrogen atom is placed along the broken bond axis such that the Si-H bond distance is 1.49 Å. For a given template and selection of cluster silicon atoms, this procedure uniquely defines the positions of the cluster-terminating hydrogen atoms. As stated above, these hydrogen atoms are held fixed during geometry optimisations, while all silicon atoms and any adsorbates are allowed to fully relax. Lastly, any exposed silicon dimers on a cluster surface are buckled into a p(2×2) pattern. Several cluster-sizes used in this work are illustrated in Fig. 3.F o r example, our one-dimer, single-row clusters of four, five, and six layers have stoichiometries of Si 9 H 12 ,Si 17 H 20 , and Si 21 H 28 , respectively. Larger clusters are simple length and width extensions of these prototype one-dimer clusters.

D. Molecular adsorption
This work focusses on adsorption energies as perhaps the most immediately useful quantity to delineate between various competing binding configurations and transition states of a molecule on a surface (such as when exploring the reaction paths and potential energy diagrams of molecular adsorption). In this we further focus on adsorption energies of a single molecule on an otherwise clean surface, which would correspond in experiment to measurements in a low-coverage limit on a good, defect-free surface sample. Adsorbate-adsorbate interactions such as those arising in higher-coverage situations create additional complexities that are beyond the scope of this work. 31 In order to gauge the effect of slab/cluster size parameters on adsorption energies we use a trial set of seven representative adsorption structures, shown in  Fig. 1, indicating up-and down-buckled dimer ends, respectively. All remaining dimers included in the slab and cluster models are buckled into a p(2×2) pattern. It is important to stress that we are not necessarily seeking to work with the preferred conformation/buckling configuration of an adsorbate, but rather focus on configurations that can be consistently adopted in both our slab and cluster models.
Molecular adsorption energies are calculated as the difference between the energy of the surfacemolecule system and the combined energies of the clean surface and the gas-phase molecule; i.e., Careful note should be taken of the sign convention used in this equation. Adsorption energies are understood here as reaction energies of adsorption; i.e. an exothermic adsorption process corresponds to a negative energy.

E. Estimated computational cost
In the following discussion we require a simple metric to compare the computational cost (i.e. CPU time) associated with slab or cluster models of different sizes. We are not after detailed measurements of the actual CPU time used in our specific calculations, which are highly dependent on numerous internal details of the DFT code, the parameter settings, and the computer used. What we seek instead is a simple, ad hoc measure to provide a generic estimate of relative cost. The metric used here is based on the nominal O(N 3 ) scaling of DFT, where N is a suitable size parameter such as the number of atoms, electrons, or basis functions. In our case we use the number of occupied atomic orbitals as the size parameter; i.e. nine for every silicon atom and one for every hydrogen atom. Thus our expression for the estimated computational cost is where N Si and N H are the number of silicon and hydrogen atoms, respectively, and C is a constant prefactor, that is different for slab and cluster calculations. We will use this cost estimate only to make relative comparisons between different sizes of either slabs or clusters; hence, the prefactor will cancel and needs not be further quantified.

A. Slab models
The quality of a slab model for describing adsorption processes on surfaces is primarily determined by two criteria: the thickness of the slab and the size of the in-plane periodic repeat (i.e. the surface unit cell). In all the cases considered here the in-plane repeat vectors run parallel and perpendicular to the Si-Si dimer rows. We will refer to these two directions in the following as the length and width of the slab model, respectively, which, together with the slab thickness, make for a total of three size parameters. We will discuss in turn the effects of these parameters on calculated adsorption energies.
Slab thickness. The effect of slab thickness, that is, the number of silicon atomic layers represented in the slab, is illustrated in Fig. 4. Adsorption energies are calculated for our trial set of seven molecules using a compact (2×2) in-plane repeat. For ease of comparison, energies are shown relative to those obtained for the thickest slab of 20 layers. The data is further condensed by plotting only the average relative energy over the trial set (solid curve) and the range of relative energies (errorbars). When plotted in this way convergence with slab thickness becomes apparent in the increasing proximity of these energies to zero. Convergence is clearly evident in Fig. 4: relative adsorption energies approach the zero line rapidly and the rate of convergence is similar for all seven molecules. This finding is further quantified by the inset in Fig. 4 in which the same data is plotted on a logarithmic scale. We see that adsorption energies for 4-, 6-and 8-layer slab models are converged to approximately 100, 20 and 10 meV, respectively. Energy differences between the two thickest slab models considered, namely 14 and 20 layers, are 2 meV or less.
Analogous convergence trends are observed when larger in-plane repeat units are used. This is illustrated in Fig. 5 for the case of the C 2 H 2 adsorbate on (2×2), (4×4), and (6×6) unit cells. Plotted using points connected by solid lines are calculated adsorption energies as a function of slab thickness. The data shows that the thickness dependence of the (4×4) and (6×6) results closely tracks the (2×2) data within the range of thicknesses tested. The results for larger unit cells are slightly offset to lower energies but otherwise run in parallel. This suggests that the effects of slab thickness and in-plane repeat are to some extent independent from one another. We note that the same behaviour is found for all seven adsorbates in our trial set.
The near-independence of assorted approximations (basis sets, exchange-correlation, etc.) has long been utilized in quantum chemistry in so-called composite or compound models to provide an economic alternative to more conventional all-in-one calculations, in which all approximations ) families of composite model chemistries for example manage the competing computational demands of describing electron correlation and basis-set completeness. This is achieved by combining the results of several calculations that are high-level in some approximations, and low-level in others, to estimate results that are high-level in all approximations. Similar composite models have been formulated for use in surface chemistry 19,35,36,39 where, typically, the size of the slab/cluster enters as an additional parameter. The computational benefit of compositing generally arises because of the high-order scaling associated with any DFT or ab initio calculation, which is O(N 3 ) or larger, where N is the number of atoms or basis functions (see Section II E). A suitable combination of several smaller calculations can thus be much less expensive than a single large, all-in-one calculation.
Following these ideas, we can utilise the apparently parallel convergence with slab thickness in Fig. 5 to estimate adsorption energies that would be calculated for thicker (4×4) and (6×6) slabs. We start with the thickest directly calculated slab for the given unit cell, and use the (2×2) results to provide a slab thickness correction (STC) to extrapolate to the desired thickness. For example, for the (4×4) unit cell, we can write, where we start with the directly calculated adsorption energies of a 10-layer (4×4) unit cell (i.e. four dimers in two rows each, denoted 4D,2R,10L), and extrapolate to any thickness nL via the 10to n-layer energy difference for the (2×2) unit cell (i.e. 2D, 1R). In the following we will denote the slab thickness correction term in square brackets as STC(10 → nL; 2D, 1R). The upper dotted line in Fig. 5 shows the effect of this extrapolation for the (4×4) case. A similar extrapolation for the (6×6) unit starts with the directly calculated 6-layer energy, and uses the (4×4) thickness dependence as follows: Note that our(4×4) data is limited to only 10 layers thickness; thus, in order to extrapolate beyond 10 layers we must draw on our (2×2) thickness dependence. This is illustrated in Fig. 5 where the two stages of extrapolation for the (6×6) slab are plotted using dashed and dotted lines. In order to gauge the effectiveness of these extrapolations, we compare in Table I adsorption energies calculated directly for two slab models -(4×4) 10-layer and (6×6) 6-layer -with those  obtained by extrapolation from thinner, 4-layer slabs. As the data shows, differences between the two approaches are 5 meV or less for the (4×4) 10-layer case, and 2 meV or less for the (6×6) 6-layer case. Such errors should be adequate for many applications. In terms of computational cost (see Section II E), the indirect extrapolation for the (4×4) 10-layer slab requires approximately 11 times less CPU time than the direct approach. This speed-up is due to the number of atoms involved: the direct approach requires a calculation on a huge Si 160 H 32 slab, whereas the largest slab in the indirect approach is Si 64 H 32 . Estimated CPU time savings in the case of the (6×6) 6-layer slab are somewhat more modest; the indirect approach is faster by a factor of three. Note however that in this case the range of extrapolation is much smaller, namely from four to six layers only. In the upcoming discussion we will identify several additional size extrapolation formulae. Slab length and width. In Fig. 6 we examine the dependence of the adsorption energy on the slab in-plane unit cell repeat. This dependence is illustrated in separate panels for slab thicknesses of 4, 6, 8, and 10 layers. For ease of comparison between molecular adsorbates, all energies are given relative to those obtained for the (4×4) unit: hence all (4×4) energies are zero by construction. The thinnest slab of four layers [ Fig. 6(a)] permits the largest in-plane repeat of (8×8), whereas thicker slabs are limited to smaller repeats. Overall, the data in Fig. 6(a) shows that convergence with unit cell size is somewhat less systematic than seen above with slab thickness. Nevertheless, the range of energy differences between the three largest models -(4×4), (6×6), and (8×8) -are within 15 meV. Larger energy differences occur between (2×2) and (4×4) units. At these differences, steric and strain interactions between the adsorbate and its periodic images become significant. In our data this is particularly evident for benzene, the largest molecule in our set, which is 70 meV less stable on a (2×2) slab than on a (4×4) slab.
The dependence of adsorption energy on cell size [ Fig. 6(a)-6(d)] appears to be independent of slab thickness; all four thicknesses considered exhibit similar changes in adsorption energy as unit cell size is increased. For example, increasing the unit cell from (2×2) to (4×4) leads to a stabilisation for three trial molecules (C 6 H 6 ,C 2 H 4 , and C 2 H 2 ), a destabilisation for three molecules This independence motivates another correction term and extrapolation formula. To estimate energy changes for an increase in unit cell size we use a slab length and width correction (SLWC). For example, extrapolating from a (4×4) 8-layer slab to a (6×6) 8-layer slab can be achieved using the (4×4) and (6×6) 6-layer results as follows, where the SLWC term corresponds to the simple energy difference E(6D, 3R, 6L) − E(4D, 2R, 6L).
In additional slab calculations we specifically focus on the length parameter, mindful of the fact that the Si(001) surface is electronically highly dispersed in this direction. Adsorption energies are calculated for a sequence of 4-layer (2×n) slabs with the length parameter, n, set to 2, 4, 8, and 16 dimers. Note that these slabs are only one dimer row in width, hence all adsorption energies obtained with this sequence of slabs are likely to include significant interactions between adsorbates and their inter-row periodic repeats. However, when used as a slab length correction (SLC; see below) these inter-row interactions are likely to cancel. Figure 7 plots the calculated energies as a function of slab model length, with energies given relative to those obtained with the largest, n=16 slab. The data shows the relative (2×2) energies are scattered over a range of approximately 70 meV. For the (2×4) and (2×8) slabs, the energy range is reduced to approximately 30 meV. This 30 meV quantifies the error associated with the slab calculations that are eight dimers in length.
Using the difference between (2×8) and (2×16) calculations, we can construct an SLC to account for slab extensions from eight to sixteen dimers. For instance, adsorption energies for a 4-layer (8×16) slab can be estimated from our 4-layer (8×8) results as follows, FIG. 8. Calculated adsorption energies as a function of cluster thickness for a single-row, 3-dimer cluster. All energies are given relative to those obtained for the 20-layer cluster. The error bars indicate the range of energy differences over the trial set, while the solid curve connects the average energy differences. The inset graph shows the same data on a logarithmic scale.

B. Cluster models
We now turn from the slab to the cluster approach for modelling molecular adsorption. Here we will examine the three parameters that define the size of our cluster models, namely the thickness in atomic layers, as well as the length and width in terms of number of dimers and dimer rows, respectively. Note that these three cluster parameters directly mirror those used for the slab models above, except that a periodic repeat no longer applies.
Cluster thickness. To examine the effect of cluster thickness, we show in Fig. 8 calculated adsorption energies for a sequence of single-row 3-dimer clusters extending up to 20 layers in depth. To facilitate comparison between different adsorbates, all energies are given relative to those obtained for the thickest cluster. The plotted curve in the data connects the average over the trial set of the relative energies for each thickness, while the error bars indicate the range of relative energies. As with slabs, convergence with thickness is evident in how this data approaches the zero line. In an inset to Fig. 8, the same data is plotted on a logarithmic scale. We see that 4-layer clusters considerably underestimate adsorption energies by an average 87 meV. Considerable improvements are achieved when 5-and 6-layer clusters are used with the average underestimate reduced to 20 and 3 meV, respectively. The notable flattening of the data in the logarithmic plot, in particular the error bar maxima at around 1 meV, is due to numerical errors other than cluster thickness becoming dominant. 37 We further examine thickness effects using a number of single-row clusters of different lengths, ranging from one to seven dimers. Shown in Fig. 9 using solid lines are calculated energies for the C 2 H 2 adsorbate. The data shows that 5-and 7-dimer clusters exhibit a thickness dependence that is broadly similar to the 3-dimer cluster. While adsorption energies are slightly offset from the shorter 3-dimer clusters, they run nearly in parallel with increasing thickness. In contrast, the smaller 1-dimer cluster converges much more rapidly with thickness. The behaviour described here for C 2 H 2 is replicated by all the adsorbates in our set. Overall this data suggests that cluster thickness and cluster length can be treated as independent parameters in terms of their effects on the adsorption energy, as long as the clusters are sufficiently long (i.e. three or more dimers).
FIG. 9. Calculated adsorption energies for C 2 H 2 using cluster models of varying thickness and dimer length. Points connected by solid lines indicate directly calculated energies, while dashed and dotted lines correspond to cluster thickness extrapolations based on 5-dimer and 3-dimer thickness dependencies, respectively (e.g. see Eq. (7)). Note that the directly calculated 3-dimer energies and the 5-dimer and 7-dimer extrapolations extend out to 20 atomic layers, i.e. beyond the range of the horizontal axis. In this equation, the first term on the right-hand side corresponds to the thickest (8-layer) 5-dimer cluster calculated directly. The second term in brackets describes a cluster thickness correction (CTC) from 8 to n layers based on the 3-dimer model. This extrapolation is illustrated in Fig. 9 using the dotted line extension of the 5-dimer cluster data. A similar extrapolation has been carried out for the 7-dimer cluster in Fig. 9: dashed lines indicate an extrapolation from six to eight layers using the directly calculated 5-dimer thickness dependence over this range. To extend to even thicker clusters, we use a CTC based on the 3-dimer results which is shown as a dotted line extension of the 7-dimer curve.
Cluster length. As in the case of the slab models, we measure the length of a cluster by the number of Si-Si dimers represented along the dimer row. Figure 10 plots calculated adsorption energies for our trial set of molecules as a function of cluster length. All energies are given relative to those obtained using the 25-dimer cluster. The curve in Fig. 10 connects the average relative energy over the trial set at each length, and the error bars indicate the range of relative energies. The data clearly displays convergent behavior. More specifically, we observe that 1-dimer cluster models FIG. 11. Calculated adsorption energy for acetylene, C 2 H 2 , on single-row cluster models of varying length and thickness. Points connected by solid lines represent directly calculated energies, whereas dashed and dotted lines indicate energies obtained using cluster length corrections (CLCs; e.g. see Eq. (8)). significantly (>100 meV) underestimate adsorption energies, i.e. they overestimate the strength of binding. Similarly, we find that 5-dimer clusters slightly yet consistently overestimate adsorption energies by approximately 20 meV. The convergence with dimer length is further quantified by the inset in Fig. 10 which shows the same data on a logarithmic scale. We see that 9-dimer clusters are converged to within 10 meV, while at least 17 dimers are required to achieve convergence of 1 meV.
In Fig. 11 we examine the cluster length dependence for clusters of several thicknesses in addition to the 4-layer cluster. Shown in this graph are adsorption energies obtained for the C 2 H 2 molecule, which exhibits behaviour that is representative of the full set of adsorbates. We directly calculate these energies for clusters up to nine dimers for the 5-layer cluster and up to seven dimers for the 6-layer cluster. Broadly, the length dependence of these two larger clusters follows that of the 4-layer cluster. There are differences, notably between 5 and 7 dimers where the 4-layer cluster displays a small decrease in adsorption energy (E 7D − E 5D =− 8 meV), while the 5-layer cluster displays an increase (+4 meV). However these differences are fairly small, and hence we can still use the dimer-length dependence of the 4-layer cluster to extrapolate thicker clusters to greater lengths. For example, to extend a 5-layer cluster from a 9-to an n-dimer cluster we can start with the directly calculated 9-dimer results, and extrapolate as follows: In this equation, we use a cluster length correction (CLC) which is again a simple energy difference, namely E(nD, 1R, 4L) − E (9D, 1R, 4L).
In Table II, we examine the performance of such extrapolations for two cases, a 15-dimer long 5-layer cluster and a 5-dimer 6-layer cluster. We use the 4-layer length dependence in both cases, starting the extrapolation from directly calculated clusters that are nine and three dimers long, respectively. The extrapolated adsorption energies are compared to those obtained by direct calculation. As the data in Table II shows, differences between the direct and extrapolated approaches are within 8 meV, which attests to the effectiveness of this approach.
Cluster width. The third cluster size parameter to consider is the number of dimer rows represented. It is worth noting that the vast majority of Si(001) cluster calculations in the literature are limited to single-row cluster models with the implied assumption being that inter-row effects are not important. We test this assumption using a number of cluster models that include three dimer rows. Table III compares the results of 1-and 3-row clusters of equivalent length and thickness. The data shows that extending a 3-dimer 4-layer cluster from one to three rows leads to adsorption energy increases of between 42 and 86 meV using our trial set. The same test conducted for 3-dimer 5-and 6-layer clusters reveals much smaller differences of less than 26 meV. At these thicknesses, the inter-row effect is certainly small, but perhaps should not be neglected when high accuracy adsorption energies are sought. The difference between single and triple row results can be regarded as a cluster width correction (CWC). For example, we could use such a correction to account for adjacent row effects when using a 5-dimer, 1-row cluster as follows,

C. Best estimate extrapolation
Having examined size parameters for slabs and clusters, we will now estimate for each approach adsorption energies that are converged in all size parameters. We start with a slab or cluster model that is suitably large and use a combination of size corrections to adjust the directly calculated energies for greater lengths, widths, and thicknesses.
Slab models. In our series of slab model calculations we have explored sizes that extend to 16 dimers, 4 rows, and 20 layers. We can therefore extrapolate to these dimensions as our best estimate of adsorption energies for a size-converged (sc) slab model. This is accomplished using the following expression: In this equation, the extrapolation starts from the (4×4) 10-layer results (i.e. 4D,2R,10L). This slab model is medium-sized in all three size parameters, which is advantageous for two reasons: firstly, the extrapolations in all parameters have less "ground" to cover, resulting in smaller and more accurate corrections, and secondly, the effects of the size changes in each direction are likely to be more independent from one another. Our extrapolation to the size-limit involves four steps, each corresponding to one of the four correction terms in Eq. (10). The first correction term is a slab thickness correction (STC) that accounts for an increase in thickness from 10 to 20 layers. This extrapolation uses the thickness dependence of a (2×2) slab as depicted in Fig. 4. The second and third correction terms are combined slab length and width corrections (SLWCs), extending the in-plane repeat in two stages from (4×4) to (6×6) and then from (6×6) to (8×8). The first of these stages uses the 6-layer length and width dependence, while the second stage uses the 4-layer dependence (cf.F i g .6). We adopt this particular two-stage approach here in order to make use of the thickest, most accurate slab results available across the extrapolation distance: using our 6-layer slab data we can only extend up to (6×6) and for the extension to (8×8) we must use the thinner 4-layer slab data. The fourth correction term in Eq. (10) further extends the slab length from (8×8) to (8×16), i.e. 16D,4R. This uses the slab length dependence obtained in our 4-layer (2×n) calculations (see Fig. 7). Table IV summarises the numerical details of these four steps for our set of trial molecules. We see that the first correction, the STC, systematically increases the adsorption energy by a small amount (4 to 5 meV) with very little variation between the adsorbates. The combined SLWC terms affect adsorption energies by between −3 and +4 meV. Lastly, the SLC term decreases adsorption energies by up to 14 meV. The extrapolated adsorption energies for the size-converged (16D,4R,20L) slab are obtained by summation over these terms as per Eq. (10), and are included in Table IV.
Cluster models. We have performed calculations on clusters with lengths up to 25 dimers, widths up to 3 dimer rows, and thicknesses up to 20 layers; hence we can extrapolate to a cluster of these dimensions to arrive at our best estimate adsorption energies within the cluster approach. This extrapolation is achieved using the following expression: We start with a single-row cluster of medium length and thickness (15 dimers and 5 layers). Our extrapolation to the size-converged cluster involves five correction steps. The first step is a cluster length correction from 15 to 25 dimers based on the 4-layer cluster length dependence (see Fig. 11). The second, third, and fourth steps are cluster thickness corrections that account for an increase in thickness from 5 to 6, 6 to 8, and then 8 to 20 layers, respectively, using clusters 7, 5, and 3 dimers long. In taking multiple CTC steps, we make best use of the largest cluster data available for each step (cf. Fig. 9 where a similar approach was taken for the 7-dimer cluster). The fifth step in Eq. (11) accounts for the fact that we have only used single-row cluster data in our extrapolation up to this point. We apply a cluster width correction (CWC) from 1 to 3 rows using our 3-dimer, 6-layer results (cf. Table III). Table IV gives numerical values for these cluster corrections for our seven trial molecules and the resulting best estimate adsorption energies. The CLCs from 15 to 25 dimers range between −1 and −2 meV, which is presumably a reflection of the fact that we commence our extrapolation with a relatively long cluster. The combined CTCs from 5 to 20 layers are an order of magnitude larger, systematically increasing the adsorption energies by an average 21 meV. However, note the spread of CTCs within the trial set, namely 3 meV, which implies that this correction has only a very minor effect on relative adsorption energies. The CWCs from one to three dimer rows are ranged from −12 to +25 meV over the set. This data shows that our starting cluster is not as converged in thickness and width as it is in length. With our approximate corrections we recover much of the remaining size effects to arrive at our best efforts adsorption energies for a size-converged cluster.
If our extrapolations had truly reached a size-converged limit, we should expect slab and cluster results to numerically agree on the adsorption energy for any given adsorbate. Hence, any residual differences between our extrapolated slab and cluster results provide a means to gauge the degree of overall size convergence. Differences between the slab and cluster approaches are reported in the bottom row of Table IV. They average −10 meV over the trial set and range from +1t o −19 meV, with the slab adsorption energies being generally more negative. This level of numerical agreement between two fundamentally different approaches is quite satisfying, suggesting that we have accounted for most of the size effects. This compares favorably to an earlier study of size effects in clusters and slabs where residual differences in GGA-DFT calculations were about an order of magnitude larger (up to 160 meV; see Table 10 of Ref. 19). Residual differences in our work are most likely associated with the limited widths considered in our calculations, namely, three and four dimer rows in clusters and slabs, respectively. In addition, k-point sampling accounts for errors of approximately 5 meV in slab models, as noted in Section II B.

D. Performance of smaller models.
Now that we have a good estimate of size-converged adsorption energies for slab and cluster models, E sc , we can use these energies to assess the accuracy of smaller models. We quantify the performance of a given model by calculating adsorption energy deviations, E = E−E sc , between this model and the size-converged model, as an estimate of error. We report this error as an average deviation, AVG( E), over the trial set of adsorbates, together with their range, RNG( E), calculated as the difference between the largest and smallest deviations. 38 Average errors thus provide a measure of any systematic under-or overestimates that may exist in the smaller models, while the range quantifies any stochastic errors that are attributable to size. Our trial set is admittedly small but should suffice to provide a general sense of model performance.
We consider first a variety of commonly used slab models, namely (2×2), (2×4), and (4×4) slabs, at thicknesses of four and six layers. Table V summarises the errors for these models. We see that the 4-layer (2×2), (2×4), and (4×4) slabs are associated with average errors of −74, −78, and TABLE V. Estimated errors, E (in meV), when calculating adsorption energies using selected slab models. Only the average error across the trial set of seven adsorbates and the range of errors are reported. Errors are calculated by subtracting from the adsorption energies the size-converged cluster energies of Table IV (see text for details). Also included is an indication of the relative computational costs associated with these models, with the 4-layer (2×2) case taken as a unit reference.

Slab
Method −79 meV, respectively, and error ranges of 89, 41, and 17 meV. So, while the stochastic/range error decreases with increasing unit size, all three models systematically, and nearly equally, underestimate adsorption energies. The underestimation is presumably due to the small number of atomic layers represented, as noted before (see discussion of Fig. 4). This is confirmed when we consider errors for 6-layer slabs: the systematic errors are much reduced, to +1 and −7meVfor(2×2) and (4×4) slabs, respectively. In contrast, the stochastic error is only slightly improved by this increase in thickness. With an average error of −7 meV and a range of 8 meV, the (4×4) 6-layer slab performs reasonably well and should be suitable for many applications. In applications where the focus is on relative adsorption energies between different adsorbates and/or configurations, the (4×4) 4-layer model is expected to perform similarly well, because the large systematic error will cancel out. The computational cost factors included in Table V indicate the expense involved when increasing slabs in size. The extension in unit cell size from (2×2) to (4×4) increases the computational cost by an estimated factor of 64, which is certainly considerable. Extending the thickness from four to six layers further increases the cost by a factor of about three. Depending on application, this may be altogether too costly to be practical. However, some savings can be realised by replacing the all-in-one, direct calculation, with a suitable composite model. For example, (2×2) and (4×4) 4-layer results may be combined with (2×2) 6-layer results to estimate (4×4) 6-layer energies as follows, The errors obtained with this extrapolation are included in Table V and they are are virtually the same as those obtained by direct calculation. Our error analysis suggests that the indirect compositing approach does not significantly affect the accuracy of the results. In terms of cost, the composite approach is considerably less expensive, requiring only an estimated one third of the resources required for the direct approach. A similar analysis can be conducted for cluster models. Table VI summarises the estimated errors for a variety of common single-row cluster sizes. The first entry in this table is the 1-dimer, 4-layer Si 9 H 12 cluster, perhaps the most widely used cluster representation of the Si(001) surface in the literature. With an estimated average error of −252 meV and an error range of 160 meV, this model is clearly very approximate. Extending this type of cluster to three dimers (Si 21 H 20 ) improves these errors by roughly a factor of two; nonetheless, this cluster still systematically underestimates adsorption energies by an average 127 meV. This three-dimer cluster can be improved in one of two ways: extending the length to five dimers, or extending the thickness to five layers. The data in Table VI shows that the extension in thickness provides a greater increase in accuracy (although, in terms of the error range, the 5-dimer 4-layer and the 3-dimer 5-layer clusters perform similarly well). A combined extension in length and thickness to give a 5-dimer 5-layer cluster (Si 65 H 52 ) provides a TABLE VI. Estimated errors, E (in meV), when calculating adsorption energies using selected cluster models. Only the average error across the trial set of seven adsorbates and the range of errors are reported. Errors are calculated by subtracting from the adsorption energies the size-converged cluster energies of Table IV (see text for details). Also included is an indication of the relative computational costs associated with these models, with the Si 9 H 12 cluster taken as a unit reference.

Cluster
Method reasonably well-performing cluster model with an average error of −29 meV and an error range of 47 meV. An increase of the cluster size to 6-layers for a 5-dimer Si 77 H 60 improves the systematic error by a factor of two to −13 meV, but leaves the error range effectively unchanged at 46 meV. This particular error range is in fact almost at the limit of what can be achieved using single-row clusters; inter-row effects as captured by our cluster width corrections alone account for an error range of 37 meV (cf. last row of Table III; difference between largest and smallest value). Improvements beyond these error ranges thus require that additional rows be included in the cluster.
Even when limited to a single row, 5-and 6-layer clusters at 5-dimers length may already be too large for typical, day-to-day applications, being respectively 321 and 531 times as expensive as the small Si 9 H 12 model according to our cost estimates (Table VI). However, clusters of this size can be effectively approximated using a composite model approach. For example, we combine results from 3-and 5-dimer clusters of 4 layers thickness with 3-dimer 5-layer results, as follows, and similarly for the 6-layer cluster, W es e ei nT a b l eVI that the energies obtained with these extrapolations are as accurate as those obtained by direct calculation. The cost is significantly reduced to 137 and 196, respectively, i.e. approximately 40% of the direct calculation, which recommends the size-composting technique for practical applications.

IV. SUMMARY AND CONCLUSIONS
In this paper we have closely examined size effects in slab and cluster models when calculating molecular adsorption energies. Using a trial set of seven molecular adsorbates, we separately analysed length, width, and thickness as size parameters. We extended our models to sizes significantly beyond those used in typical day-to-day applications in order to seek out their convergence limits. Our results provide evidence that these size parameters are to some degree independent from one another in their effects upon adsorption energies. We used this idea to combine results from our largest models and provide an estimate of slab and cluster energies that are size-converged in all three parameters.
Additionally, our use of DFT-GGA and localised atomic basis sets via the DMol 3 software means that our slab and cluster results were obtained on an equal footing. This facilitates direct comparison between these two approaches in that any numerical energy differences can be confidently attributed to finite size effects. Moreover, the dictum that energies from slabs and clusters should agree in a size-converged limit provides a potent measure of how close our estimates are to this limit. Residual differences between slab and cluster estimates were found to be an average 10 meV only, with none larger than 19 meV.
Needless to say, slab and cluster size are not the only parameters that govern the accuracy of calculated adsorption energies. In day-to-day application, a practitioner will have to balance the available computational resources over a range of additional approximations that are determined by the system at hand. These include basis set completeness and exchange-correlation. Accumulated evidence in the literature (e.g. Refs. 19 and 39) suggests that the effects of high-order exchangecorrelation treatments (relative to GGA-DFT) and larger basis sets on the adsorption energy are localised to the immediate adsorption site. This implies that high-order treatments can be carried out using smaller slabs/clusters, and combined with size-correction terms derived from GGA-DFT for an accurate calculation of adsorption energies. Recent applications of such composite models in surface science can be found for example in Refs. 39, 40, and 20. Analogous ideas underpin the ONIOM approach [41][42][43] in quantum chemistry to describe localised chemical reactions within a larger system. Our results reported in Section III D support such composite modelling for Si(001) in two ways. First, our accuracy assessments of smaller slab and cluster sizes relative to the size-converged limit (see Tables V and VI) can be used to match a given exchange-correlation model to an appropriately sized slab/cluster such that the estimated errors associated with each are about the same. Second, with our slab/cluster size extrapolations we provide a strategy to break down large calculations into smaller, more economic fragments.
To illustrate these ideas, consider a sophisticated G2 exchange-correlation treatment with an estimated error of approximately of 23 meV (Ref. 44). Combining this exchange-correlation model with a small Si 9 H 12 cluster would clearly result in a rather unbalanced effort, given the estimated average cluster size error of 252 meV (see Table VI). A more even match requires a much larger cluster, such as the (5D,1R,5L) Si 65 H 52 cluster with an error of 26 meV. That said, a full G2 level treatment using a cluster of this size is unlikely to be practical. This is where our size extrapolations come into play. An alternative, less costly approach is to limit the high-level G2 treatment to the small cluster, and account for the cluster size increase from Si 9 H 12 to Si 65 H 52 using GGA-DFT via the cluster length and width correction terms described above. Similarly, a hybrid-DFT B3LYP treatment with a large 6-311+G(3df,2p) Gaussian-type basis set has a reported average error of about 98 meV (Ref. 44). This level of theory is evenly matched by a (5D,1R,4L) Si 33 H 28 cluster with a size error of 93 meV (see Table VI). The costly B3LYP treatment could again be limited to a small cluster while the extension to Si 33 H 28 is handled via DFT-GGA cluster size corrections. The specific accuracy needs will inevitably be determined by the scientific problem considered. Our comprehensive survey of size effects provides a quantitative measure of accuracy to assist in the rational design of cluster and/or slab models for the task at hand.

ACKNOWLEDGMENTS
This work is supported by the Australian Research Council Centre of Excellence for Quantum Computation and Communication Technology (project number CE1100096). The computational resources used in this project were provided by the Australian National Computational Infrastructure (NCI).