CO on Pt(111): A puzzle revisited

Today’s state-of-the-art method for calculating the interaction of atoms or small molecules with metal surfaces is considered to be density functional theory (DFT) at the generalized gradient approximation (GGA) level employing a slab or supercell representation of the surface. The method is widely used and by many assumed to be both qualitatively and quantitatively accurate. This notion has recently been challenged by Feibelman et al. [J. Phys. Chem. B 105, 4018 (2001)] who suggest that the DFT/GGA method does not correctly predict the most stable adsorption site for the CO/Pt(111) system, and they conclude that the method is not qualitatively accurate. However, using a different calculational approach we find a good agreement between the calculated potential energy surface for this system and the one inferred from experiments, indicating that the evidence supporting the view of Feibelman et al. is not yet conclusive. On the contrary, we advocate the view that the DFT/GGA method should at the moment be c...


I. INTRODUCTION
The interaction of atoms and small molecules with metal surfaces is of great academic interest, but is also extremely important in industry and society. Corrosion and heterogeneous catalysis ͑among others͒ are determined by these processes. The experimental and theoretical effort focused on understanding the interactions in order to correctly model and eventually control the reactions is therefore considerable ͑see, e.g., Refs. 1-4, and references therein͒.
Developing a fundamental understanding of the reactions on an atomic level and helping to interpret the results of increasingly complex experiments are two of the many roles theory plays. It is therefore essential to have an idea of how accurate today's state-of-the-art theoretical methods areblindly following the path set out by these theoretical calculations might leave a large number of steps to be retraced to re-establish ''firm ground.'' The state-of-the-art method of today for calculating the interactions between atoms/molecules and metal surfaces is considered to be density functional theory ͑DFT͒ at the generalized gradient approximation ͑GGA͒ level employing a slab or supercell representation of the surface ͑see, e.g., Refs. 5-11͒. One reason is that the results of a slab/supercell calculation prove rather easy to converge with respect to the number of layers used ͑see, e.g., Refs. 5, 12-15͒, whereas the cluster approach shows considerable difficulties in converging the results with respect to the size and the shape of the cluster ͑see, e.g., Refs. 12, 16 -21, but note that the very recent results of Ref. 22 indicate that it might be feasible͒. Another reason is that the introduction of the GGAs ''fixed'' the wrong site-preference for CO adsorption on Cu͑100͒, 5,12 and the lack of reaction barriers to dissociation for the H 2 /Cu system [6][7][8] as calculated within the local density approximation ͑LDA͒ of DFT.
During the last 8 years the notion that DFT/GGA calculations can accurately model the interactions between atoms/ molecules and metal surfaces both on a qualitative and quantitative level has become prevalent in our community ͑see, e.g., Refs. 4, 23-25͒. This can be understood partly on the basis that only a few studies have emphasized the difference between experiments and DFT/GGA calculations ͑see, e.g., Refs. 26 -28͒, and partly from the evidence that authors questioning the accuracy of the method at the same time continue to present convincing studies with a successful application of the same method ͑see, e.g., Refs. 25-30͒. However, a recent paper of Feibelman et al. 28 has clearly challenged the above notion: Based on an extensive set of calculations using different plane-wave based electronic structure codes they suggest that DFT/GGA clearly predicts the wrong adsorption site for the CO/Pt͑111͒ system and put this forward as a qualitative failure of the DFT/GGA method.
Conclusions on the accuracy of today's state-of-the-art theoretical methods could have far reaching consequences for our understanding and modeling of, e.g., heterogeneous catalysis, and we think that it is important to base these conclusions on as broad a basis as possible. This has prompted us to employ a linear combination of atomic orbitals ͑LCAO͒ approach to obtain DFT/GGA results for the CO/Pt͑111͒ system using a slab representation of the surface. Our results indicate that the DFT/GGA method correctly predicts the most stable adsorption site for this system. Based on this result and a number of other considerations we will offer a different, rather more positive, view than given in Ref. 28 on the qualitative and quantitative accuracy of the DFT/GGA method. Another motivation behind this study is that the convergence of previous calculations using the LCAO method 31 has been questioned in Ref. 28, and, considering the issue at hand, we think it is important to address this in some detail. We should also note that we are not the only ones being motivated by the ''puzzle'' presented in Ref. 28: The issue has already been addressed in Ref. 32. In the following section the details of the employed method will be given together with results of extensive convergence tests. Our results for the CO/Pt͑111͒ system follow in Sec. III, and these results are discussed in relation to experimental results and theoretical contributions employing different calculational methods in Sec. IV. Section V concludes.

II. CALCULATIONAL DETAILS
The program ADF-BAND ͑Refs. 33-35͒ was used to solve the Kohn-Sham equations 36,37 self-consistently for a CO molecule adsorbed on a Pt͑111͒ surface, modeling the surface by a slab with translational symmetry in two directions. A combination of numerical atomic orbitals obtained from numerical Herman-Skillman-type calculations 38 and Slatertype orbitals forms a flexible basis set that has been used in the expansion of the one-electron states. The calculations can in principle be performed with all electrons included in the variational space. However, deeper lying electronic levels are kept frozen when explicit tests show that their variational inclusion has negligible effect. Pseudopotentials are never used. The calculation of the matrix elements of the Hamiltonian is performed by product Gauss numerical integration rules within the pyramids that constitute the Voronoi polyhedra surrounding the atoms. 35 ͑The Voronoi polyhedron of an atom is the part of space that is closer to that atom than to any other one; in regular crystals these become the Wigner-Seitz cells.͒ A single parameter ͑the real space integration parameter͒ governs the precision of the numerical integration; a value q corresponds to relative precision 10 Ϫq in the integrals. Finally, the k-space integration can be done accurately using the quadratic tetrahedron method. 39 The Vosko-Wilk-Nusair formulas 40 are used to calculate the exchange-correlation energy in the LDA. In this study we have employed two GGAs. The first combines the Becke correction 41 for the exchange energy with the Perdew correction 42 for the correlation energy ͑BP͒, and the second is the gradient-corrected functional of Perdew et al. 43,44 ͑PW GGA-II, which we will label PW for brevity͒. Both scalar relativistic and spin-orbit corrections are included through the zeroth-order regular approximation ͑ZORA͒. 31,45,46 Since our main goal is to decide on the relative stability of adsorption sites that are close in energy, we have performed extensive tests to ensure that our results are properly converged. The results of these tests are shown in Fig. 1: ͑a͒ Setting the real space integration parameter ͑see above and Ref. 35͒ to 5.0 results in an accuracy of 10-20 meV. ͑b͒ With a k-space integration parameter 39 set to 5 ͑corresponding to at least 15 points in the irreducible wedge of the surface Brillouin zone͒ the error is within 20-30 meV. ͑c͒ A five ͑three͒ layer slab gives results that differ by less than 10-20 ͑30-40͒ meV from a seven or eight layer slab. ͑d͒ The TZ2P basis set ͑triple zeta, i.e., three basis functions per valence orbital, plus two polarization functions͒ gives results that are within 10-20 meV of the results based on the QZMP basis set ͑quadruple zeta, i.e., four or more basis functions per valence orbital plus four or more polarization functions͒. The QZMP basis set is very large, with proven accuracy of bond energies at the 10 meV level. 47 The different basis sets tested are shown in Table I. All results remain unchanged ͑within the accuracy limit given below͒ upon unfreezing the 1s orbital of C and O, and the 4d or deeper lying orbitals of Pt. It is worth noting that full convergence with respect to the basis set can be obtained without the use of plane waves, as has been demonstrated in the past. 34 We have also carefully checked that the calculations are converged to better than 10 meV with respect to two other computational aspects, ͑i͒ the density-fitting procedure used to represent the deformation density and to evaluate a part of the periodic Coulomb potential and ͑ii͒ the evaluation of the Madelung sums. 34 Assuming that all the error sources tested above are statistically independent, the chosen settings for the calculations indicate that the results for the relative stability of the CO adsorption sites given in the next section are converged to about Ϯ50 meV of the theoretical model limit for the CO/ Pt͑111͒ system. ͑a͒ convergence with increasing real space integration parameter, ͑b͒ convergence with increasing k-space integration parameter, ͑c͒ convergence with increasing number of Pt layers, and ͑d͒ convergence with increasing size of basis sets ͑the basis sets are given in Table I͒. All results are calculated for a CO molecule adsorbed in a )ϫ)-R30°structure on one side of the Pt͑111͒ slab employing the BP approximation. Parameters used in the following are indicated by arrows.

III. RESULTS
For each of the three sites top, fcc, and bridge the equilibrium geometry has been determined ͑Table II͒ for a CO molecule adsorbed in a )ϫ)-R30°structure on one side of a 3 layer slab ͑frozen at a Pt-Pt distance equal to the experimental bulk lattice constant, 5.24 bohr͒ with its molecular axis normal to the surface ͑a five layer slab gives identical results for the geometries within the accuracy of our calculations͒. The results have been obtained through a twodimensional cubic spline interpolation on a grid consisting of more than 25 points in the distance between the C and the O atom ͓denoted d(C-O)] and the distance of the C atom above the Pt͑111͒ surface plane ͓denoted d(Pt-C͔͒. The results of the interpolation has ͑through normal mode analysis͒ also been used to extract the vibrational frequency for the CO internal stretch and the vibration of the entire CO molecule against the surface ͑Table II͒.
Using the above obtained geometries the relative energy differences between the sites have been calculated for a five layer slab and the results are given in Table III. Results are given for a slab frozen at a Pt-Pt distance equal to the experimental bulk lattice constant, and a slab with a GGA optimized bulk lattice constant and a top layer relaxed along the direction of the surface normal.
Before entering a discussion of these results in relation to experimental studies or theoretical studies employing dif-ferent calculational methods, we would like to comment on previous results obtained with the ADF-BAND program. In Ref. 31 it was shown how the ZORA efficiently can be implemented for extended systems, and the adsorption of CO molecules on the ͑111͒ surface of Ni, Pd, and Pt was selected as test systems to illustrate the importance of the relativistic effects in different rows of the periodic system. The calculations were performed on a )ϫ)-R30°structure with CO on one side of a two layer slab employing a basis set of the same quality as the one labeled TZ1P in Table I, using a real space integration parameter of 4.0 and a k-space integration parameter of 5. In these calculations the top site was preferred by 0.24 eV ͑when including scalar relativistic or scalar relativistic plus spin-orbit corrections͒ above the hcp site, whereas our present calculations indicate a preference for the top site above the fcc site of only 35 meV ͑note that in Ref. 31 all reported results were given for a slab frozen at a Pt-Pt distance equal to the experimental bulk lattice constant, and we therefore compare these results to the ones we have obtained for a slab with the same Pt-Pt distance͒. From Fig. 1 we see that the difference in computational parameters can only partly account for the difference in the relative stability of the top and hollow sites found in Ref. 31 and this study. The remaining part of the difference can, however, be understood on the basis that the adsorption geometries used in the two studies are not the same. In Ref. 31 the CO bond was frozen to the gas-phase value of 2.15 bohr, and the CO distance to surface was obtain with an accuracy of about 0.2 bohr. Since the top site adsorption geometry was closer to the true equilibrium than the hollow site geometry this led to an overestimation of the preference for the top site ͑compare Table III of Ref. 31 and Table II above, and note that we have checked that the difference between the fcc and hcp sites is negligible͒. But even if the energetic preference for the top site above the hollow site was too large by about 0.2 eV ͑at the scalar relativistic or scalar relativistic plus spinorbit level͒, the main conclusion of Ref. 31 regarding the CO/Pt͑111͒ system remain valid: In Table III we have included nonrelativistic results and they indicate a preference for the fcc above the top site by 330 meV ͑at the BP level͒. Thus, it is important to include relativistic corrections for correctly predicting the most stable adsorption site for CO on Pt͑111͒ ͑more comments on this will follow in the next sec-tion͒.

IV. DISCUSSION
From the results in Table III we see that the energy differences are small. This is in accordance with experiments which suggest that the potential energy surface governing the diffusion of CO on Pt͑111͒ is rather flat ͓the difference in adsorption energy between the top and bridge site is estimated to be 30-60 meV ͑Refs. 48, 49͔͒. Our ''best'' results indicate that the top site is preferred by 100 ͑70͒ meV above the hollow site at the BP ͑PW͒ level. We also get the same site-preference sequence as inferred from experiments, top followed by bridge and then fcc, but this could be said to be fortuitous considering the error bar of Ϯ50 meV. However, the level of agreement with experiment is certainly gratifying and the error we make seems to be well within Ϯ0.1 eV. Based on our present results and a number of other considerations we would like to offer a different, rather more positive, view on the qualitative and quantitative accuracy of the DFT/GGA method than given in Ref. 28. Support for this different view is given by the following considerations.
In Ref. 28 a preference of 0.10 eV for the fcc site above the top site was found using the full-potential linearized augmented plane wave ͑FP-LAPW͒ method at the DFT/PW level. Both the FP-LAPW and the LCAO approach are in principle able to solve the Kohn-Sham equations without approximations, and we should therefore expect to get the same results. This is not the case, but we see at least two possible sources for the discrepancy: ͑i͒ We have used the ZORA to describe the relativistic effects, treating all electrons at the scalar relativistic or the scalar relativistic plus spin-orbit level. In the FP-LAPW calculations the Dirac equation is solved for the core orbitals and the semicore and valence orbitals are treated at the scalar relativistic level. Both approaches are approximate, and since it is clear that the relativistic effects are important for the CO/Pt͑111͒ system, we think there is more work to be done at this level before the final conclusions can be reached. In this respect it is worth to note that recent cluster model calculations employing a four-component relativistic DFT method indicate that the top site is the most stable CO adsorption site on Pt͑111͒ ͑even though the authors admit to possible problems with the convergence of the cluster size͒. 50 We should also address the difference in the magnitude of the relativistic effects found here and in Ref. 28. Table III indicates that relativistic corrections stabilize the top site by almost 0.4 eV relative to the fcc site, whereas a stabilization of only 0.1 eV is found in Ref. 28. This difference is mainly due to two different ways of defining the relativistic effects: In our case the nuclear geometries are kept the same and the relativistic effects are purely electronic. In the case of Ref. 28 the nuclear geometries are optimized at each level of relativistic approximation, and the relativistic effects therefore contain changes of nuclear geometries as well as electronic effects. ͑ii͒ As shown in Sec. II, a considerable effort was spent to ensure that the calculations presented here are converged to within Ϯ50 meV. An error bar associated with the computational parameters has not been given for the FP-LAPW calculations, but if it is similar to or larger than the Ϯ50 meV we report, the discrepancy might be somewhat smaller than it appears. However, it should be conceded that also for the LCAO and FP-LAPW approaches there are many computational parameters that have to be carefully controlled in order to obtain results that are true to the DFT functional model to an accuracy better than Ϯ50 meV. As one example it is worth noting that we have used the Vosko-Wilk-Nusair formulas 40 to parameterize the LDA functional, whereas the Perdew-Zunger formulas 51 were used in Ref. 28. At the moment it is not known if part of the discrepancy could be explained by this difference, or by other differences in the computational methods, but this needs to be evaluated if in the future an accuracy of 10 meV, or eventually 1 meV, is to be attained.
With respect to the pseudopotential methods upon which many of the studies of the interaction between atoms/ molecules and metal surfaces are based, we wish to make the following cautioning comments. In a recent study we obtained the same remarkable degree of agreement for the CO/ Cu͑100͒ system as in the present study between experiment and LCAO based DFT/GGA calculations regarding the sitepreference and diffusion barrier. 30 In particular, the previously reported 52 discrepancy with experiment for the potential energy surface describing diffusion of CO on the Cu͑100͒ surface, which had been calculated with a pseudo-TABLE III. The stability of the adsorption sites relative to the top site for the two GGAs, for an unrelaxed surface with the experimental bulk lattice constant ͑5.24 bohr͒ ͑''unrelaxed''͒, and a surface with a BP optimized bulk lattice constant ͑5.32 bohr͒ with top layer relaxed along the direction of the surface normal ͑''relaxed''͒. All results are calculated for a CO molecule adsorbed in a )ϫ)-R30°structure on one side of the Pt͑111͒ slab. The energy differences are given in meV and a positive number indicates that the top site is the more stable. Calculations including no relativistic, scalar relativistic and spin-orbit corrections are labeled nr, sr, and so, respectively. All results should be understood with an error bar of Ϯ50 meV.  potential plane wave method, proved to disappear when employing the LCAO approach. In Ref. 53 a preference for the hollow site was found for the CO/Cu͑111͒ system, which is in conflict with the experimental indications that the top site is the most stable adsorption site. However, from Table III in Ref. 28 it is clear that the pseudopotential approximation can be off by as much as 0.15 eV for CO adsorption on Pt͑111͒with a similar error for the CO/Cu͑111͒ system the DFT/ GGA error would probably not be larger than Ϯ0.1 eV, even though a study avoiding the use of the pseudopotential approximation would be needed to definitively establish the magnitude of a possible discrepancy between experiments and DFT/GGA calculations for this system. It is also interesting to consider the results of Ref. 32 for the CO/Pt͑111͒ system, which were obtained with a range of pseudopotentials with varying effective radii. The results ranged from having the top and fcc sites almost equally stable ͑with pseudopotential radii chosen somewhat too large͒ to a preference for fcc by 0.08 eV with the most confined pseudopotentials. This presumably most accurate result still differs appreciably from the comparable ͑with the same PBE func-tional͒ pseudopotential result of Ref. 28 which found a preference for the fcc site by 0.24 eV. This indicates the same possible 0.15 eV error noted earlier due to the pseudopotential approximation.

Site
Regarding the data base of reliable computational results for CO adsorption and diffusion on metal surfaces we would like to make the following additional remarks. In the cases of CO/Rh͑111͒ ͑Ref. 54͒ and CO/Ru͑0001͒ ͑Refs. 24, 55͒ the correct preference for the top site has been found. For the CO/Ni͑111͒ ͑Ref. 26͒ and CO/Pd͑111͒ ͑Ref. 56͒ systems DFT/GGA calculations and experimental results agree on the fcc site being the most stable. Together with the results for CO/Pt͑111͒ and CO/Cu͑111͒ we therefore think there are enough indications that the DFT/GGA error when considering the relative stability of the CO adsorption sites falls within Ϯ0.1 eV. With this in mind it is interesting to note that the CO/Rh͑111͒ system is in Ref. 28 indicated as another case for which the DFT/GGA approach gives the wrong result-experiments indicate the top site to be the most stable, whereas the DFT/GGA calculations prefer the hcp site. However, in Ref. 54 the previously unpublished results by Mavrikakis and Nørskov quoted in Ref. 28 are contradicted, and a preference for the top site above the hcp site by 0.05 eV is found. Furthermore, as we have seen from our past and present results of the LCAO method ͑Sec. III͒, the preference for the top site for the CO/Pt͑111͒ system has been reduced from 0.24 eV to 35 meV by improving the convergence of the calculations and the geometries used. Thus, we think that the accuracy with respect to the computational parameters with which atom/molecule-surface interactions energies are obtained does not at the moment allow one to make inferences beyond the Ϯ0.1 eV level.
An additional consideration in this respect is that when the energy differences become smaller than say 0.1 eV, it becomes important to compare more quantities than the energy differences alone before reaching a conclusion on the accuracy of the applied method. 57 For CO adsorption on metal surfaces quantities that can be measured experimen-tally and also calculated theoretically are, e.g., the vibrational frequency for the internal CO stretch and the vibration of the CO molecule against the surface at the equilibrium adsorption site. In the case of CO/Pt͑111͒ these measured frequencies are 260 meV ͑Ref. 58͒ and 58 meV, 59 respectively. From the results of Table II above or Table VII in Ref. 28 we see that the calculated frequencies for CO adsorption at the top site agrees within 5 meV of the experimental result, suggesting that the adsorption site assigned by theory should be the top site. This approach would in the other cases mentioned above also have led to theory and experiment agreeing on the adsorption site assignment.
Taking all existing evidence together it seems that the DFT/GGA error falls within Ϯ0.1 eV for the relative stability of the CO adsorption sites on metal surfaces, and there are enough indications that it is difficult to converge DFT/ GGA calculations in all computational parameters ͑real space integration, k-space integration, number of layers in the slab/supercell, basis sets, spherical harmonic expansions, lattice sums, modelling of relativistic effects, pseudopotentials, etc.͒ to better than Ϯ0.1 eV. We think this provides support for the view that the DFT/GGA method should at the moment be considered qualitatively accurate for predicting the most stable CO adsorption sites on metal surfaces.
We would like to stress that the magnitude of the DFT/ GGA error suggested by the above considerations only pertains to the relative stability of the CO adsorption sites on metal surfaces. A number of other quantities that also can be considered to be relative energy differences, like, e.g., the CO adsorption energy on a metal surface and the CO molecular binding energy, show larger discrepancies between the theoretical calculated and the experimental measured ones. In this study we find a CO adsorption energy at the top site of 1.45 and 1.57 eV ͑at the BP and PW level, respectively, when including scalar relativistic corrections͒. This is 0.3-0.4 eV lower than the 1.86Ϯ0.08 eV measured in Ref. 60. For the CO molecule we find the equilibrium distance to be 2.15 bohr ͑BP, PW͒, a binding energy of 11.51 ͑11.79͒ eV at the BP ͑PW͒ level, and the CO vibrational frequency to be 264 meV ͑BP͒ and 266 meV ͑PW͒. The calculated geometry and vibrational frequency compare favorably to the experimental values, 2.13 bohr and 269 meV, respectively, but the calculated binding energy is off by 0.4 -0.7 eV compared to the measured one of 11.09 eV. 61

V. CONCLUSIONS
A density functional theory ͑DFT͒ calculation employing the generalized gradient approximation ͑GGA͒ and a slab or supercell representation of the surface is today's state-of-theart method for calculating the interaction of atoms or small molecules with metal surfaces. The method is widely used and it is by many assumed to be both qualitatively and quantitatively accurate. Recently Feibelman et al. ͓J. Phys. Chem. B 105, 4018 ͑2001͔͒ have challenged this view. Their extensive set of plane wave based electronic structure calculations suggest that the DFT/GGA method predicts the wrong adsorption site, the fcc site, for the CO/Pt͑111͒ system, and they see this as a qualitative failure of the DFT/GGA method.
Changes in the views on the accuracy of the state-of-theart method of today could possibly have far reaching consequences for our understanding and modeling of, e.g., heterogeneous catalysis, and should be based on as broad a basis as possible. This has been our motivation when applying a linear combination of atomic orbitals ͑LCAO͒ approach employing a slab representation of the surface to obtain DFT/ GGA results for the CO/Pt͑111͒ system. We find that the calculated potential energy surface is in good agreement with what is inferred from experiments, with the top site being the most stable adsorption site. This indicates that there is not yet conclusive evidence supporting the view that the DFT/ GGA method is qualitatively failing for the CO/Pt͑111͒ system.
On the contrary, we advocate the view that the DFT/ GGA method should at the moment be considered qualitatively accurate for predicting the most stable CO adsorption sites on a given metal surface. The precise content of this statement can be made explicit in the following observations: ͑i͒ Our results for the CO/Pt͑111͒ system together with results in the literature for the CO/Cu͑100͒, CO/Cu͑111͒, CO/ Ni͑111͒, CO/Rh͑111͒, and CO/Ru͑0001͒ systems, suggest that the DFT/GGA error for the relative stability of the adsorption sites stays within Ϯ0.1 eV when compared to experimental results. ͑ii͒ There is ample evidence to show that converging DFT/GGA calculations to better than Ϯ0.1 eV with respect to all computational parameters is a considerable task. The present data base of reliable computational results neither excludes nor proves convincingly that the DFT/GGA model for CO site-preference on metal surfaces disagrees at this accuracy level with the experimental data base. ͑iii͒ These first two statements refer to the relative stability in the adsorption energies. An additional argument supporting our favorable view of the DFT/GGA method is that when the energy differences are smaller than say 0.1 eV, other quantities with site differences beyond the computational error bars, like, e.g., vibrational frequencies and geometries, discriminate correctly between sites, i.e., the theoretical quantities at the correct site exhibit agreement with experiment.
Even though this paper indicates that the CO/Pt͑111͒ system might be less of a puzzle than recently suggested ͓J. Phys. Chem. B 105, 4018 ͑2001͔͒, there is clearly more work to be done at the theoretical level in validating the different approximations used for including relativistic effects in periodic DFT/GGA calculations.