On the accuracy of first-principles lateral interactions: Oxygen at Pd(100)

0.57 MB·English

by  Yongsheng Zhang

(submitted to Phys. Rev. B)
On the accuracy of first-principles lateral interactions: Oxygen at Pd(100)
Yongsheng Zhang, Volker Blum, and Karsten Reuter
Fritz-Haber-Institut der Max-Planck-Gesellschaft, Faradayweg 4-6, D-14195 Berlin, Germany
(Received 19th January 2007)
7
We employ a first-principles lattice-gas Hamiltonian (LGH) approach to determine the lateral
0
interactions between O atoms adsorbed on the Pd(100) surface. With these interactions we ob-
0
tain an ordering behavior at low coverage that is in quantitative agreement with experimental
2
data. Uncertaintiesintheapproacharisefrom thefiniteLGHexpansionandfrom theapproximate

n exchange-correlation (xc) functional underlying the employed density-functional theory energetics.
a Wecarefullyscrutinizetheseuncertaintiesandconcludethattheyprimarilyaffecttheon-siteenergy,
J
whichrationalizes theagreement withtheexperimentalcritical temperaturesfortheorder-disorder

3 transition. We also investigate the validity of the frequently applied assumption that the ordering
2 energiescanberepresentedbyasumofpairterms. RestrictingourLGHexpansiontojustpairwise
lateral interactions, we findthat thisresults in effective interactions which contain spurious contri-
] butionsthatareofequalsize,ifnotlargerthananyoftheuncertaintiese.g. duetotheapproximate
i
c xcfunctional.
s
- PACSnumbers: 68.43.Bc,68.43.De
l
r
t
m
I. INTRODUCTION performed.6,7,8,9,10,11,12,13 To make contact with a spe-
.
t cificmaterialandwithexperiment,wespecificallychoose
a
theon-surfaceadsorptionofoxygenatPd(100),forwhich
m Lateral interactions between species adsorbed at solid
detailed experimental data on the ordering behavior is
d- sthuerfatcaersgaetreocfruscuiraflamceicsrcoisecnocpeicsqtuudainetsitfieosrtahaltohnagvetibmeee.n1 available14. Sinceinthissystemhigheroxygencoverages
n above Θ 0.5monolayers [ML, defined with respect to
These interactions govern both the equilibrium, as well ∼
o the number of Pd atoms in one layer of Pd(100)] induce
as the non-equilibrium ordering behavior of the adsor-
c structures containing incorporated oxygen14,15,16,17,18,
[ bates, and thereby critically influence the surface func-
tion and properties in important applications like het- we concentrate on the low coverage regime. For this
1 regime,twoorderedstructureshavehithertobeencharac-
erogeneous catalysis. Traditionally, considerable efforts
v terizedexperimentally14,15,16,17,19,20: Ap(2 2)-Ostruc-
9 have been devoted to determine lateral interactions em- ture at 0.25ML and a c(2 2)-O structur×e at 0.5ML,
4 pirically from experimental data, e.g. from temperature ×
both with O adsorbed in the on-surface fourfold hol-
5 programmed desorption or low energy electron diffrac-
low sites. The material science motivation of our first-
1 tion (LEED) measurements. In order to simplify the in-
principlesLGH study is thento extractthe lateralinter-
0 herentlyindirectdeterminationfromsparseexperimental
7 actions operating between the adsorbed O atoms at the
data, the assumptionof exclusively pairwise interactions
0 surface and to study the ordering behavior they imply.
between the adsorbed species has often been applied. In
/ Specifically, this is to see whether we can confirm the
t recentyears,algorithmicadvancesandincreasedcompu-
a experimentally determined orderedstructures,as well as
m tationalpowerhavemadeitpossibletodeterminethelat-
thecriticaltemperaturesfortheorder-disordertransition
eralinteractionsalternativelyfromfirst-principles. Most
- in the low coverageregime.
d notably, these are approaches that parameterize lattice-
n gasHamiltonians(LGHs)with density-functionaltheory Presenting a systematic first-principles lattice-gas
o (DFT) energetics, also called cluster expansions.2,3,4,5 Hamiltonianexpansion,weindeedfindthecalculatedset
c Since the accuracyofthe determined lateralinteractions of lateral interaction energies to be fully consistent with
:
v should be of the order of kBT to properly describe the theexperimentallyreportedlowcoveragephasediagram.
Xi thermalordering,a concernwith this approachhas been Critically discussing the uncertainties of our approach,
whether the employed first-principles energetics is actu- both with respect to the employed LGH expansion and
r
a ally accurate enough. the underlying DFT energetics, we conclude that they
Within this context, the present work has a method- primarily affect the on-site energy. The lateral interac-
ologicalandamaterialssciencemotivation. Themethod- tionenergies,ontheotherhand,canbe determinedwith
ological motivation is to scrutinize both the assumption quite high accuracy, which we estimate for the present
of exclusively pairwise interactions, and the accuracy system to be around 60meV. Comparing these interac-
withwhichthefirst-principlesLGHapproachcanprovide tion energies with those determined previously empiri-
the lateral interactions. For this purpose we concentrate cally and using the pairwise interaction approximation
onasimplemodelcase,namelytheon-surfaceorderingof demonstrates that the latter assumption introduces an
atomic adsorbatesata (100)cubic surface,for whichex- error that is at least as large as this remaining uncer-
tensivestudieswithmodelinteractionshavealreadybeen taintywhencarefullydeterminingthelateralinteractions
2
from present-day first-principles calculations.
II. THEORY
A. Lattice-Gas Hamiltonian
Inordertodescribethe site-specificadsorptionofoxy-
genatomsonthePd(100)surfaceweemploytheconcept
of a lattice-gas Hamiltonian, in which any system state
is defined by the occupation of sites in a lattice, and
the total free binding energy of any configuration is ex-
panded into a sum of discrete interactions between the
latticesites(seee.g. Refs. 2,3,4,5). Foraonecomponent
system with only one site type, this energy reads (with
obvious generalizations to multi-site cases)
FIG.1: (Coloronline)TopviewofthePd(100)surface, illus-
r
NF = Fon−site n + Vp n n + tratingtheconsidered poolof 17lateral interactionsbetween
b b i u i j O atoms in on-surface hollow sites. Vp (m = 1,2,3,4,5)
Xi uX=1 (iX<j)u arethetwo-body(orpair)interactionsamtfirst,second,third,
q fourth and fifth nearest neighbor distance. Vt, Vq, Vqu are
n n n
Vt n n n +... , (1) considered compact trio, quattro and quintointeractions, re-
u i j k
uX=1 (i<Xj<k)u spectively. Light grey (yellow) spheres represent Pd atoms,
and small dark grey (red) spheres O atoms.
where the site occupation numbers n = 0 or 1 indicate
l
whether site l in the lattice is empty or occupied, with
a total of N sites occupied, and Fon−site is the free en- B. Static and vibrational average binding energy
b
ergy of an isolated species at the lattice site, including
static and vibrational contributions. There are r pair In order to generate a quantitatively accurate LGH,
interactions with two-body (or pair) interaction energies we parameterize the unknown lateral interaction ener-
Vp between species at uth nearest neighbor sites, and q giescontainedintheLGHbyfirst-principlescalculations.
u
triointeractionswithVt three-bodyinteractionenergies. The centralquantities requiredfor this parameterization
u
The sum labels (i<j) [and (i<j <k) ] indicate that are computed average free binding energies for a set of
u u
the sums run over all pairs of sites (ij) (and three sites ordered configurations of O adsorbed at Pd(100). We
(ijk)) that are separated by u lattice constants, and the write this average free binding energy as
summationis donesuch,that eachpair(ortrio)ofoccu-
F (T) = E +Fvib.(T) , (2)
piedsitescontributesexactlyoncetothelatticeenergy.21 b b b
Formally, higher and higher order interaction terms separating the total and vibrational contributions, E
b
(quattro,quinto,...) followintheinfiniteexpansionofeq. and Fvib.(T) respectively. The former is defined as
b
(1). In practice, the expansion must (and can) be trun-
catedafterafinitenumberofterms. Obviously,thejudi- 1 N
E = Etotal Etotal OEtotal .
cious choice of which interactions to consider and which b −N (cid:20) O/Pd(100)− Pd(100)− 2 O2(gas)(cid:21)
O
onesto neglectmustcriticallyaffectthe reliabilityofthe (3)
entire approach. To quantify the impact of this choice Here N is the total number of adsorbed O atoms,
O
onaccuracy,werelyontheconceptofleave-one-outcross Etotal , Etotal , and Etotal are the total ener-
validation(LOO-CV)detailedbelowtoidentifythemost O/Pd(100) Pd(100) O2(gas)
gies of the surface containing oxygen,of the correspond-
important interactions out of a larger pool of possible
ing clean Pd(100) surface, and of an isolated oxygen
interactions. Figure 1 illustrates the lateral interactions
molecule, respectively. Since a free O molecule is thus
2
containedinthispool,whichrangefrompairinteractions
used as the zero reference for E , a positive binding en-
b
up to the fifth nearest neighbor, via all trio interactions
ergy indicates that the dissociative adsorption of O is
2
up to second nearest neighbor, to several compact quat-
exothermic at T =0K.
tro and one quinto interaction. The pool focuses thus
In order to determine the vibrational contribution to
onshort-tomedium-rangedinteractions. Interactionsat
the average free binding we use the phonon density of
larger distances can be substrate-mediated elastic or of
states σ(ω) and write the vibrational free energy as
electronic origin22, but for the present system we do not
expect such interactions to play a role on the accuracy
Fvib.(T) = dω F(T,ω)σ(ω) , (4)
level of interest to this study. Z
3
where ToobtainthetotalenergyoftheisolatedO molecule,
2
we exploit the relation Etotal =2Etotal D, where
1 1 O2(gas) O(atom)−
Fvib.(T,ω) = h¯ω + Etotal is the total energy of an isolated oxygen atom,
(cid:18)2 eβh¯ω 1(cid:19) O(atom)
− andD thetheoreticalO bindingenergy. TheisolatedO
2
β¯hω
k T ln(1 e−βh¯ω (5) atomisthencalculatedspin-polarized,insidearectangu-
− B (cid:20)eβh¯ω 1 − − (cid:21) lar cell of side lengths (12 13 14)bohr, Γ-point sam-
− × ×
pling of the Brillouin zone and without spherically aver-
is the vibrational free energy of an harmonic oscillator
aging the electron density in the open valence shell. For
of frequency ω.23 k is the Boltzmann constant, and
B D we employ the previously published ultra-converged
β = 1/(k T) the inverse temperature. The vibrational
B GGA-PBE value of 6.21eV.26
contribution to the average binding energy can then be
For the calculations of the adsorbate vibrational
written in exactly the same way as eq. (3), namely
modes, the dynamical matrix is set up by displacing the
1 Oatomfromitsequilibriumpositionin0.05˚Asteps. An-
Fvib.(T) = Fvib. (T) (6)
b −NO h O/Pd(100) − ticipatingagooddecouplingofthevibrationalmodesdue
to the largemassdifference between PdandO,the posi-
N
Fvib. (T) OFvib. (T) = tions of all atoms in the substrate below the adsorption
Pd(100) − 2 O2(gas) (cid:21)
site are kept fixed in these calculations. The frequen-
1 dω Fvib.(T,ω) σ (ω) cies and normal modes are then obtained by subsequent
−N Z O/Pd(100) − diagonalizationof the dynamic matrix.
O (cid:2)
N
O
σ (ω) σ (ω) .
Pd(100) − 2 O2(gas) (cid:21)
D. Monte Carlo simulations
To evaluate this contribution in practice one must thus
determine the difference of the surface phonon density Onceareliablesetofinteractionshasbeenestablished,
of states of the adsorbate covered and of the clean sur- evaluating the LGH for any configuration on the lattice
face, σ (ω) and σ (ω) respectively, as well corresponds merely to performing an algebraic sum over
O/Pd(100) Pd(100)
as the vibrational frequencies of the gas phase molecule a finite number of terms, cf. eq. (1). Due to this sim-
contained in σ (ω). plicity, the LGH can be employed to evaluate the sys-
O2(gas)
tem partition function. Here this is done by canonical
Monte Carlo (MC) simulations for O coverages up to
C. Total Energy Calculations Θ=0.5ML.Theemployedlatticesizewas(40 40)with
×
periodic boundary conditions. Metropolis sampling used
The total energies required to evaluate eq. (3) are 2000MCpassesperlatticesiteforequilibration,followed
obtained by DFT calculations within the highly accu- by 10000 MC passes per site for averaging the thermo-
ratefull-potentialaugmentedplanewave+localorbitals dynamic functions. Increasing any of these numerical
(LAPW/APW+lo) scheme24, using the generalized gra- parameters led to identical results on the accuracy level
dient approximation (GGA-PBE)25 for the exchange- of interest to this study, i.e. here primarily critical tem-
correlationfunctional. Allsurfacestructuresaremodeled peratures that are convergedto within 5-10K.
inasupercellgeometry,employingfully-relaxedsymmet- For fixed coverage on the surface, ordered structures
ric slabs (with O adsorption on both sides of the slab) areidentifiedbyevaluatingorderparameterssensitiveto
consisting of five (100) Pdlayerswith an optimized bulk lateralperiodicities. Tocheckonthe(2 2)periodicityof
×
lattice constant of a = 3.947˚A (neglecting bulk zero- the twoexperimentally characterizedorderedstructures,
pointvibrations). Avacuumregionof 10˚Aensuresthe we divide the (100) cubic lattice into four interpenetrat-
≥
decoupling of consecutive slabs. The LAPW/APW+lo ing sub-lattices a, b, c and d in a (2 2) arrangement
×
basis set parameters are listed as follows: Muffin tin ab . This allows to separately evaluate in the MC runs
cd
spheres for Pd and O are RPd = 2.1bohr and RO = (cid:0)the(cid:1)averagenumber of occupied sites in eachsub-lattice,
MT MT
1.1bohr,respectively,thewavefunctionexpansioninside N , N , N , and N , respectively. Using Fourier theory,
a b c d
the muffin tin spheres is done up to lwf = 12, and the the order parameter for the p(2 2) structure is then
max ×
potential expansion up to lpot = 6. The energy cutoff defined as
max
for the plane wave representation in the interstitial re- 3
tghioenwbaevtewfeuenncttihoensmaunffidnEtpinotsp=he1r9e6sRisyEfomwrfatxh=e p2o0tRenytifaolr. Ψp(2×2) = 4NtotqN12+N22+N32 , (7)
max whereN =N +N +N N ,N =N +N N +N ,
Monkhorst-Pack (MP) grids are used for the Brillouin 1 a b c− d 2 a b− c d
N =N N +N +N , and N is the total number
zoneintegrations. Specifically,weusea(12 12 1)grid 3 a− b c d tot
× × ofsitesinthesimulationcell. Inthesameway,the order
for the calculation of (1 1) surface unit-cells. For the
× parameter for the c(2 2) structure is defined as
larger surface cells, care is taken to keep the reciprocal ×
spacepointsamplingidenticalbyappropriatelyreducing 1
the employed k-meshes. Ψc(2×2) = Ntotq(Na+Nd−Nb−Nc)2 . (8)
4
O atoms in Section IIID.
TABLE I: Calculated binding energies E (in meV/O atom)
b As the next step in the LGH parameterization, aver-
for O adsorption in on-surface hollow or bridge sites. The
age binding energies for different ordered configurations
nomenclatureandgeometricarrangementinthesurfaceunit-
withO atomsinon-surfacehollowsites andwithsurface
cell for thefiveordered adlayers are explained in Fig. 2.
unit-cells up to (3 3) were correspondingly computed.
(3×3)-O p(2×2)-O c(2×2)-O (2×2)-3O (1×1)-O Despite our focus o×n the lower coverage regime, the set
Θ 0.11ML 0.25ML 0.50ML 0.75ML 1.00ML doescomprisestructureswithcoveragesupto1ML,since
hollow 1249 1348 1069 643 344
thesestructuresarerequiredduringtheLGHparameteri-
bridge 1024 961 801 573 378
zationto determine in particularthe higher-ordermany-
body interactions occurring in (locally) denser adatom
arrangements. In most cases configurations that we ini-
Computing these order parameters as a function of tem- tiallypreparedwithhighcoverageofon-surfaceOatoms
perature, the critical temperature for the order-disorder weretrulymetastable,inthesensethattheyrelaxedinto
transitionisdefinedbytheinflectionpointwhereΨp(2×2) geometrieswhereallOatomsremainedinthe on-surface
orΨc(2×2) gotozero. In parallel,we alsoderive the crit- hollow sites. However, some structures directly relaxed
ical temperature from evaluating the specific heat, ob- into geometries with O incorporated below the first Pd
tainingvaluesthatareidenticaltowithin10Kwiththose layer. Since the O-O interaction in these structures does
inferred from the order parameters. not correspond to the physical situation we want to de-
scribe(andwouldthusmessupthechosenLGH),weex-
cludedtheseconfigurationsfromourset(usingastrongly
III. FIRST-PRINCIPLES LATTICE-GAS enlarged Pd first interlayer distance as criterion). This
HAMILTONIAN FOR O AT PD(100) left a set of 25 configurations with on-surface O atoms,
which was subsequently used in the parameterization of
the LGH. The binding energy data of this set is com-
A. Energetics of on-surface adsorption
piled in Fig. 3, and indicates already overall strongly
repulsive lateral interactions, which reduce the binding
Owingtothetendencyofoxygenatomstopreferhighly
energy with increasing coverage by up to 1eV.
coordinated binding sites at late transition metal sur-
faces, the high-symmetry fourfold hollow sites appear as
themostlikelyadsorptionsitesatPd(100). Ontheother
B. Lateral interactions
hand,onecannotexcludeapriorithatthetwofoldbridge
sitesarenotalsometastable,i.e. localminimaofthe po-
tential energy surface. To test this, we slightly displaced In order to determine the lateral interaction energies
abridgesiteOadatominap(2 1)configurationlaterally in the LGH, we employ eq. (1) to write down the LGH
towards a neighboring hollow×site. The resulting forces expression for each of the ordered configurations calcu-
relaxed the adatom back to the ideal bridge position, so lated by DFT (including the interactions with the pe-
that at least in this configuration the bridge site is not riodic images in the neighboring cells). Neglecting the
just a mere transition state, i.e. a saddle point of the vibrational contributions in eq. (2) for the moment, we
potential energy surface. As this might also be true for equate the right hand side of eq. (1) with NOEb for
bridgesiteadsorptioninother(local)Oadatomarrange- the corresponding configurations and arrive at a system
ments, we calculated the binding energetics of O atoms of linear equations that can be solved for the unknown
inthefourfoldhollowandinthetwofoldbridgesitesfor5 values of the lateral interaction energies. The crucial as-
different ordered overlayersspanning the coverage range pects of this procedure are therefore i) the number and
up to one ML. The periodicities of these overlayers are type of interactions included in the LGH expansion,and
explainedforthecaseofhollowsiteadsorptioninFig. 2, ii) the number and type of ordered configurations that
and Table I summarizes the calculated binding energies. arecomputedwithDFTtodeterminethe valuesofthese
interaction energies. In the following we show how i)
For lower coveragesthe fourfold hollow site is energet-
is addressed by leave-one-out cross-validation, and ii) is
icallyclearlymorestable,whichsuggestsaninsignificant
aidedby asearchfor the LGH “groundstate” structures
contribution of bridge sites to the ordering behavior at
and an iterative refinement of the input structure set.
coverages up to around 0.5ML, even if the latter are al-
waysmetastablesites. Althoughthereversaloftheener-
getic order betweenhollow and bridge sites at Θ=1ML
seen in Table I is intriguing, it clearly occurs in a cover- 1. Leave-one-out cross-validation (LOO-CV)
age range where surface oxide formation and eventually
three-dimensional oxide cluster growth takes place.17,18 In a truncated LGH expansion with finite ranged in-
Sinceourinterestliesintheon-surfaceorderingbehavior teractions,sparseconfigurationswillexhibita latticeen-
atlowcoverages,wewillthereforefocusthe LGHexpan- ergy that is simply N times the on-site energy, as soon
O
sion for the moment exclusively on adsorption into the as all adsorbed species have distances from each other
fourfoldhollowsites,andreturntothe roleofbridgesite thatexceedthelongestrangeconsideredinteraction. We
5
FIG. 2: (Color online) Top view of 5 ordered adlayers with O in on-surface hollow sites. The coverage of each configuration
from left to right panel is 1/9, 1/4, 1/2, 3/4 and 1 ML, respectively. Light grey (yellow) spheres represent Pd atoms, small
dark grey (red) spheres O atoms, and theblack lines indicate thesurface unit-cells.
m) values of the consideredlateralinteractions are obtained
400
to from least-squares fitting to the DFT energies of the re-
a
O 600 maining M 1 calculated configurations, i.e. leaving
V/ exactly the i−th configuration out of the fit. This way,
e
m 800 the CV score is intended to be a measure of the predic-
(
y tive power of a LGH expansion considering a given set
g
er1000 of lateral interactions. In general, one would expect sets
n
E containing too few interactions to be too inflexible and
g 1200
n thus leading to a high CV score,whereas sets containing
ndi1400 too many interactions as loosing their predictive power
Bi through overfitting and thereby also leading to a high
0 0.2 0.4 0.6 0.8 1
Θ (ML) CV score. The hope is thus to identify the optimum set
of considered interactions as that set that minimizes the
CV score.
FIG.3: (Coloronline)Coverage(Θ)dependenceofthecalcu-
latedDFTbindingenergiesof25orderedconfigurationswith Within this approach, we evaluate the CV score for
O atoms in on-surface hollow sites. any set of interactions out of the larger pool of 17 lat-
eral interactions shown in Fig. 1. Table II summarizes
therefore fix the on-site energy Eon−site in eq. (1) to be thesescoressubdividedintotheoptimumsetscontaining
b m lateral interactions, i.e. listed are the sets that yield
the DFT binding energy computed for 1/9ML coverage,
the lowest CV score for any arbitrary combination of m
cf. Fig. 2 and Table I. In this particular (3 3)-O con-
× lateral interactions out of the total pool of 17. For these
figuration, the minimum distance between O adatoms is
sets, we then determine the values of the considered lat-
8.37˚A,i.e. sixnearestneighborsitesaway. Thisislarger
eralinteractionsbyleast-squaresfitting tothe computed
than the farthest reaching interaction contained in our
DFT binding energies of all M = 24 ordered configura-
pool of lateral interactions, cf. Fig. 1, so that fixing
Eon−site to the 1/9ML (3 3)-O binding energy should tions and also include them in Table II. The minimum
b × CV score reached indeed decreases initially upon adding
prevent fitting noise into this parameter.
more interactions to the set, and then increases again
Togetsomeguidanceastowhichcouldbe the leading
for sets containing more than 10 interactions. Another
lateralinteractionstobeincludedintheLGHexpansion,
gratifying feature is that almost always the same inter-
weestimatethepredictivepowerofthe LGHbythecon-
actions are picked out of the pool, i.e. the optimum set
cept of leave-one-out cross-validation4,27,28,29,30. For a
for m + 1 interactions corresponds mostly to those in-
given set of LGH interactions the cross-validation (CV)
teractions already contained in the optimum set for m
score is calculated as
interactions plus one additional one. Only very rarely is
an interaction that is contained in the optimum m set
M
CV = v 1 EDFT(i) ELGH′(i) 2 . (9) not selected in the optimum m+1 set. And if this hap-
uuM Xi=1(cid:0) b − b (cid:1) pens, this concerns lateral interactions for which only
t
very small values are determined, and which are thus
Here, the sum runs over the remaining M = 24 ordered anyway not meaningful within the uncertainties of our
configurations that were calculated with DFT (apart approach. Alsoinaphysicalpicture,thedeterminedval-
from the 1/9ML coverage structure used already for ues for the lateral interactions appear quite plausible for
Eon−site), and which have DFT calculated binding en- m up to 11. The pairwise interactions decrease with in-
b
ergiesEDFT(i). ThequantityELGH′(i)ofthe ithconfig- creasing distance, and the leading higher-order trio and
b b
uration, on the other hand, is evaluated from the LGH quattro interactions are smaller in size than the most
expression for this configuration, cf. eq. (1), where the dominant nearest-neighbor pair interaction. The quinto
6
TABLE II: List of the sets containing m lateral interactions, together with their CV scores and the values determined for the
lateral interaction energies (no entry at a position in thetable means that this interaction is not contained in theset. Lateral
interactions shown in Fig. 1, but not shown here are never selected out of the pool). Negative values for the interaction
energiesindicaterepulsion,positivevaluesattraction. ThesetsshownarethosethatminimizetheCVscoreamongallpossible
sets containing m lateral interactions out of the pool of 17 shown in Fig. 1. The first line for each set corresponds to the
data obtained by fitting to 24 ordered DFT configurations, while the second line is obtained after adding 2 additional DFT
configurations to the fit (see text). Unitsfor the CV score and lateral interactions are meV.
m CV pair trio quattro
score Vp Vp Vp Vp Vp Vt Vt Vt Vt Vt Vt Vt Vq Vq
1 2 3 4 5 1 2 3 4 5 6 7 1 2
7 31 −344 −130 52 14 −120 132 −54
32 −338 −126 50 16 −114 117 −51
8 20 −324 −126 50 16 −138 120 −57 102
22 −314 −122 46 18 −138 108 −57 120
9 16 −292 −90 50 10 −168 60 −48 −51 120
16 −296 −92 50 10 −162 63 −48 −51 114
10 17 −290 −92 50 10 −168 60 −48 1 −51 120
17 −298 −92 50 10 −162 69 −51 −54 −42 120
11 18 −284 −90 46 12 6 −171 57 −51 −3 −54 132
18 −296 −92 48 10 2 −165 69 −48 −54 −42 126
12 20 −292 −70 50 8 2 −168 486 −264 414 −36 −216 132
19 −296 −92 48 10 2 −165 69 −51 −3 −54 −48 126
interaction contained in the pool of possible interactions 2. Ground state search
is never selected. In contrast, the m = 12 set shows al-
ready clear signs of linear dependencies, with some trio Before moving to the ordering behavior at finite tem-
interactions suddenly exhibiting very large values. This peratures, a crucial test for the validity of the LGH ex-
continues for expansions containing even more interac- pansionis thatit givesthe correctorderedgroundstates
tions(notlistedinTableI),whichexhibitmoreandmore at T = 0K, i.e. the lowest-energy structures at a given
obviously meaningless lateral interaction energies. coverage. Here, this refers in particular to the ground
states predicted by the DFT energetics, since the lat-
ter is the input with which the LGH expansion must be
consistent. Obviously, if the energetic order of compet-
Another equally important feature of the expansions ing configurations is wrong at the DFT level (e.g. due
up to m = 11 is the stability of the determined lateral to the employed approximate xc-functional), there is no
interaction values against adding further interactions to hopethatacorrectLGHexpansioncouldcurethisprob-
the set. In particular for the optimum sets close to the lem. To this end, itis useful to replotthe DFT database
overallCV minimum, i.e. for m equal to 9 or 10, adding showninFig. 3informofaformationenergyplot4. For-
another lateral interaction to the set changes the values mation energies ∆Ef are in general defined as an excess
forthe dominantinteractionsbylessthan2meV.Asim- energy with respect to the equivalent amounts of pure
ilar behavior is obtained for another test to which we constituents. For the present case of on-surface O ad-
subject our expansion: Having calculated another two sorption in the coverage range below 1ML we therefore
DFT configurations, we increased the set of DFT con- define
figurations employed in the fit to M = 26 and repeated
1
the entire CV score evaluation. The results are also in- ∆E = Etotal (1 Θ)Etotal
cluded in Table II and show only minimal variations for f Ntot h O/Pd(100)− − Pd(100)−
the sets up to m = 11. Almost always the same lateral ΘEtotal . (10)
(1×1)−O/Pd(100)
interactionsarepickedoutofthe poolandalsotheir val- i
ueschangebylessthan10meVcomparedtotheprevious
As in eq. (3), Etotal is the total energy for a spe-
procedureemploying24DFT configurations. Forthe set O/Pd(100)
cific adsorbate configuration with N O atoms per sur-
m = 12, adding the two DFT input energies removes O
face unit-cell (correspondingto a coverageΘ=N /N
the linear dependencies and brings the set back in line O tot
with N the number of sites per surface unit-cell),
with the other sets with fewer interactions. These find- tot
Etotal is the total energy of the clean surface, and
ings suggest that the expansions are also stable against Pd(100)
adding further DFT configurations and we finally iden- E(t1o×ta1l)−O/Pd(100) isthetotalenergyofthefullmonolayer
tify the set containing 9 lateral interactions and using (1 1)-O configuration. With this definition, ∆E re-
f
×
24 DFT configurations to determine their values as our flects the relative stability of a particular configuration
optimum LGH expansion. with respect to phase separationinto a fractionΘ of the
7
O coverages up to 1ML. The corresponding data points
are also shown in Fig. 4. If we first focus on the cover-
age range up to 0.5ML, we find the obtained LGH data
to be fully consistent with the DFT ground state line.
Namely, there is no structure predicted by the LGH ex-
pansion that would have an energy that is lower than
the DFT convex hull, and the LGH Hamiltonian yields
therefore exactly the same ordered ground states as the
DFT input data.
In fact, this is not a coincidental result demonstrating
the reliability of the achieved LGH expansion, but the
end product of an iterative procedure. We used the con-
sistencywiththeDFTgroundstatelineasanothercrite-
rion to judge whether more DFT ordered configurations
FIG. 4: (Color online) Formation energies ∆Ef as computed are required as input to the LGH expansion. Initially
with DFT for the 25 ordered configurations shown in Fig. 2
wehadstartedwithasmallernumberofDFT configura-
(big (red) circles), as well as for all configurations obtained
tionsasthesetdiscussedabove. Havinggonethroughthe
by direct enumeration and using the LGH expansion (small
sameCVscoreevaluation,wehadidentifiedanoptimum
(black) circles), see text. The red line is the convex hull for
LGHexpansion,buthadthenobtainedLGHdatapoints
theDFTdata,identifyingthreeorderedgroundstatesshown
as insets: p(2×2)-O at Θ=0.25ML (left inset), c(2×2)-O inthedirectenumerationleakingbelowtheDFTground
at Θ = 0.5ML (top inset), and (2×2)-3O at Θ = 0.75ML state line. Interpreting the corresponding structures as
(right inset). Pd = large, light (yellow) spheres, O = small, important input to the LGH expansion, we would ide-
dark (red) spheres. allycalculate themwith DFT andaddthemto the DFT
data set used for the LGH expansion. This was, unfor-
tunately, not always possible, when the structures pre-
full monolayer configuration and a fraction (1 Θ) of dicted by the LGH had surface unit-cells that exceeded
−
clean surface, and we can relate it to the binding energy our computational capabilities. In such cases, we looked
of the configuration by for other structures in smaller unit-cells, which still con-
tained what we believed were the relevant motifs and
∆Ef = Θ Eb,O/Pd(100)−Eb,(1×1)−O/Pd(100) . (11) computedthosewithDFT.This procedurewasrepeated
(cid:2) (cid:3)
several times, each time adding new structures to the
Plotting ∆E versuscoverageas done in Fig. 4 allows
f
DFT data base, until the present consistent result was
to identify the ground states, i.e. lowest-energy ordered
obtained.
phases, that are predicted by the present DFT data set,
as those lying on the convex hull (or so-called ground
state line).4 Any structure that yields a formation en-
ergythatishigherthanthisgroundstatelineisunstable In the coveragerangeabove 0.5ML, the situationis
∼
against decomposition into the two ordered configura- not that perfect. As apparent from Fig. 4 there are still
tions represented by the two closest lying corner points several LGH structures slightly below the DFT ground
on the convex hull. As apparent from Fig. 4 the convex stateline. Unfortunately,furtherimprovementalongthe
hull formed by the DFT data set exhibits three ordered sketchedlinesisinhibitedbytheabovedescribedpropen-
ground states (apart from the trivial ones at the ends of sity of structures in this coverage range to directly relax
the considered coverage range). Consistent with exist- into geometries with O incorporated below the first Pd
ing experimental data, these are the p(2 2)-O ordered layer. This renders it very tough to provide new on-
×
phase at Θ=0.25ML, and the c(2 2)-O orderedphase surface O/Pd(100) structures to the data base and im-
×
at Θ=0.5ML. A third ordered structure, (2 2)-3O at proveonthepresentLGHexpansion. Althoughnotcom-
×
Θ=0.75ML, is at best metastable, since it falls already pletelysatisfying,wethereforecontendourselveswiththe
in the coverage range above 0.5ML, for which surface achieved expansion. Particular care should therefore be
∼
oxide formation sets in. exerted, when aiming to use it to describe the higher
Using eq. (11) we can also evaluate formation ener- coverage regime, since denser adatom arrangements can
gies using the binding energies obtained from the LGH presumablynotbefullyreliablydescribed. However,due
expansion. Since the evaluation of the latter is numeri- totheoverallstronglyrepulsiveinteractions,thelocaloc-
callysignificantlylessdemanding,wecansampleamuch currence of such denser arrangementsat lower coverages
larger configuration space in this case. To this end, we is rather unlikely in the MC simulations. Correspond-
directly enumerate all combinatorially possible ordered ingly, we do expect the results obtained from our expan-
structuresinsurfaceunit-cellsofanysymmetryandwith siontobereliableforthecoveragerangebelow 0.5ML,
∼
a surface area smaller or equal to a (4 4) cell and with on which we focus in the present work.
×
8
1200 absolute temperature values and the trend of increasing
critical temperatures with coverage. The largest devia-
1000 tionofabout250KresultsatΘ=0.25ML,wheretheory
predicts a small peak in the critical temperature, which
800 is absent in the experimental data and which we dis-
)
K cuss in Section IVB below. Apart from this feature, the
(
T 600 agreementwiththeexperimentalTc(Θ)dataisquitesat-
isfying, if not quantitative.
400
D. Population of bridge sites
200
0.1 0.15 0.2 0.25 0.3 0.35
θ (ML)
TheLGH+MCsimulationsuptonowhaveexclusively
focusedonOadsorptionintothe fourfoldon-surfacehol-
FIG.5: Θ−T diagramofthecriticaltemperaturesTc forthe lowsites. Thealreadygoodagreementobtainedwithex-
order-disorder transition. Shown with stars are the experi- isting experimental data, together with the significantly
mental data from Ref. 14, and with solid circles the values lower stability of the energetically next favored high-
obtained when using the optimum LGH expansion. Addi-
symmetry bridge sites apparent in Table I, seem to sug-
tionallyshownbyemptycirclesarethevaluesobtainedwhen
gest that the on-surface low coverage ordering can in-
usingtheLDAasexchange-correlationfunctionalintheDFT
deed be understood in terms of the most stable hollow
calculations (see text).
sites only. To verify the implied negligible population of
(possiblymetastable)bridgesites,evenuptothe critical
temperaturesoftheorder-disordertransition,weproceed
C. Order-disorder transition
by including these sites into the LGH expansion. Since
the intention is at this point only to check on the influ-
Havingestablishedthegroundstateorderedstructures enceofapopulationofthesesites,weconsiderareduced
we proceed to study the ordering behavior at finite tem- poolofpossiblelateralinteractionsbetweenOatomsad-
peratures. Experimentally, this was investigated in the sorbedinbridgesites,consistingoftheequivalentstothe
coverage range up to 0.6ML by Chang and Thiel14. For hollow-hollowpair and trio interactions shownin Fig. 1.
definedinitialcoveragesatthesurface,theyidentifiedthe Due to the twofold symmetry of the bridge sites, two
presence of ordered phases at the surface by monitoring differentformsatthe sameinteratomicdistanceexistfor
LEED superstructure spots corresponding to the differ- someoftheinteractions,andinadditionthereisalateral
ent periodicities, andthe criticaltemperatures T (Θ) for interactionVp at the very shortdistance ofa/2 between
c 0
the order-disorder transition were determined by the in- O atoms sitting in immediately adjacent bridge sites co-
flectionpointofthe vanishingspotintensities atincreas- ordinated to the same Pd atom.
ing temperatures. Avoiding the O-induced reconstruc- However, when computing with DFT configurations
tion at higher coverages, we focus here on the data in containing such closely neighboring O atoms at a site
the coverage range Θ < 0.35ML, in which the p(2 2) distance of a/2, we always found them to be unstable
×
or a coexistence of p(2 2) and c(2 2) phases form against relaxation. During the geometry optimization
× ×
the ordered structures at low temperatures, cf. Fig. 4. the O atoms moved to sites further apart, indicating
For this coverage range, Chang and Thiel determined a sizable repulsive Vp. Focusing therefore on 18 or-
0
the onset of desorption in their ultra-high vacuum ex- dered configurations that do not contain O adatoms at
periments at much higher temperatures than the order- such close distance, the best sets with a varying total
disorder transition.14 From this we assume that in the number of lateral interactions are determined via LOO-
experiments,the coverageatthe surfaceremainedessen- CV in the same manner as for the hollow-hollow inter-
tiallyconstantattheinitiallypreparedcoveragevaluefor actions. Similar to the results in Table II, the differ-
all temperatures up to the critical temperatures for the entexpansionsyieldconsistentlythe sametwodominant
order-disordertransition. lateral interactions, namely the first-nearest neighbor
The experimental conditions are simulated by canoni- pair Vp(bridge bridge) and second-nearest neighbor
cal MC runs for fixed coverages and at various temper- pair V1p(bridge −bridge) interactions. Both are largely
atures. With the definitions in eqs. (7) and (8), our repulsi2ve, with−Vp(bridge bridge) 400meV and
order parameters are equivalent to LEED spot intensi- Vp(bridge bridg1e) 12−0meV. Tu≈rnin−g to the lead-
2 − ≈ −
ties, so thatthe determined criticaltemperatures for the inglateralinteractionsbetweenOatomsadsorbedinhol-
order-disordertransitioncanbedirectlycomparedtothe lowandbridgesites,weagainfoundthatstructureswith
experimentalvalues. Figure5showsthe T (Θ)curveob- O atoms in directly adjacent bridge and hollow sites at
c
tained with the optimum LGH expansion together with the very short distance of a/√8 are not stable against
a reproduction of the experimental data. Overall, we structural relaxation. As for the other interactions, this
observe very good agreement, both with respect to the time the first-nearest neighbor pair interaction and one
9
compact trio interaction turn out to be dominant in the A. Uncertainties in the LGH expansion procedure
LOO-CV based LGH expansion procedure. They are
also largely repulsive, with values of 240meV and
Table II provides detailed information about the in-
≈ −
280meV, respectively.
fluence of most approximations in the LGH expansion
−
procedure. InspectingthebasicallyindistinguishableCV
We therefore approximately consider the interactions
score for the expansions with m =9, 10, and 11, one
involving bridge site O atoms in the MC simulations by
might take the scatter in the correspondingly extracted
excluding configurations with O atoms at the above de-
lateral interactions as a rough measure for the uncer-
scribedshortesta/2anda/√8distancesbetweenbridge-
taintyintroducedbytruncatingtheLGHexpansionafter
bridge and bridge-hollow site pairs, respectively. In ad-
a finite number of terms. Concerning the finite number
dition, the two next dominant lateral interactions for
ofDFTconfigurationsemployedintheparameterization,
bridge-bridge and bridge-hollow pairs are explicitly ac-
the achieved consistency of the DFT and LGH ground
counted for using the values stated above. The consec-
stateline illustratedinFig. 4givessomeindicationasto
utive MC simulations show virtually no change in the
the minimum number of configurations that is required.
groundstateorderedstructuresandcriticaltemperatures
Correspondingly,thedifferencesinthelateralinteraction
in the coverage range up to Θ = 0.35ML. The overall
values determined when extending this minimum set by
largely repulsive interactions together with the signifi-
two further configurations (upper vs. lower line for each
cantly less stable on-site energy compared to adsorption
expansion in Table II) may be taken as reflecting the
in the hollow sites, efficiently prevents any significant
uncertainty due to employing a finite number of DFT
population of bridge sites. For all temperatures up to
configurations in the parameterization.
the order-disorder transition, we find less than 10% of
This leaves as remaining ad hoc feature of our expan-
the available bridge sites occupied with O atoms, with
sion procedure the decision to not include the on-site
the highestpopulationsobtainedforthelargercoverages
energy into the fit, but instead fix it at the value of the
in the considered range. To make sure that these results
sparsestDFT configurationcomputed, i.e. the (3 3)-O
are not affected by the uncertainty in the approximately
×
structure at 1/9ML. To this end, we redid all LGH ex-
determined lateral interactions, we varied the value of
pansions in Table II using the same LOO-CV procedure
each of the four lateral interactions by 100meV and
± described above, but this time including the (3 3)-O
each time reran the MC simulations. This had no influ-
×
structure into the set of DFT configurations and includ-
ence on the findings, so that we do not expect them to
ing the on-site energy Eon−site into the fit. The results
beinvalidatedbythecrudewayofhowthebridge-bridge b
areremarkablyconsistent,inthe sensethattheobtained
and bridge-hollow interactions are considered. Instead,
lateralinteractionsdiffer inallcasesbylessthan15meV
weconcludethatuptocoveragesofΘ 0.35MLandup
≈ from the values compiled in Table II, and for the expan-
to the critical temperatures for the order-disorder tran-
sions with m = 9, 10, 11 the on-site energy is in fact
sition, a population of (possibly metastable) bridge sites
determinedatvalues thatarewithin 15meVofthe com-
plays no role for the on-surface ordering behavior.
puted binding energy of the (3 3)-O structure. For ex-
×
pansions using fewer lateral interaction figures (m < 9)
this becomes progressively worse, and the increasing in-
flexibilityofthefew-interactionexpansionsstartstoshift
errors between the on-site energy and the lateral inter-
actions in an uncontrolled way. We therefore conclude
that for expansions with a sufficient number of interac-
IV. ACCURACY OF FIRST-PRINCIPLES
LATERAL INTERACTIONS tion figures it apparently makes very little difference to
include or not include the on-site energy into the fit; the
expansions are stable in this respect. In view of the sig-
The agreement with the experimental low coverage nificantly different inaccuracies in the determination of
phasediagram(groundstate structuresandcriticaltem- the on-site energy and lateral interactions discussed be-
peratures) suggests that the determined set of first- low, we nevertheless prefer to fix the on-site energy at
principles lateral interactions is quite reliable. In order the value of the sparsest DFT configuration computed.
to get a better understanding of the explicit uncertain- Summarizingthediscussionontheuncertaintiesinthe
ties of the different parameters in this set, we return to LGH expansion and parameterization procedure, we es-
a critical discussion of all approximations entering the timatetheuncertaintyintroducedbythevariousapprox-
LGH approach,and scrutinize their influence on the lat- imationstobeoftheorderof10-20meVinthedominant
eral interaction values. Uncertainties arise on the one lateral interactions. When using the first-principles lat-
hand side due to the truncated LGH expansion and the eral interactions to determine quantities describing the
finite number of configurations employed to parameter- mesoscopic ordering behavior, this uncertainty needs to
ize it, and on the other hand due to the approximate be taken into account. Figure 6 illustrates this for the
first-principles energetics,both with respect to total and determined critical temperatures for the order-disorder
vibrational free energy contributions. transition by comparing the results obtained for the dif-
10
800 eq. (1) solely to the N E of the different computed
O b
configurations. The formallycorrectprocedurewouldbe
to evaluate for each configuration also the vibrational
free energy contribution to the binding energy,Fvib.(T),
600 b
and consider it in the parameterization. As apparent
)
K from eq. (6), what determines this vibrational contri-
(
T bution are the changes of the vibrational modes during
400 adsorption. Since this will predominantly concern the
adsorbate-derivedvibrationalmodes,weestimatetheor-
der of magnitude of this contribution by the zero point
energy (EZPE) difference arising from the change of the
200 O stretchfrequency,ω ,totheO-substratestretch
0.1 0.15 0.2 0.25 0.3 0.35 2 O2(gas)
θ (ML) mode, i.e. we approximate eq. (6) roughly with
1 N
FIG.6: Θ−T diagramofthecriticaltemperaturesTc forthe Fbvib.(T) ∼ −N (cid:20)EOZP/PEd(100) − 2OEOZP2(Egas)(cid:21) (12)
O
order-disorder transition. Compared are the curves obtained
whenusingthefirst-principleslateralinteractionsfordifferent 1 NO ¯h 1
LGHexpansionscompiledinTableIIandeachtimeusing24 ∼ −N 2 (cid:20)ωO/Pd(100)(i) − 2ωO2(gas)(cid:21) ,
DFT configurations in the parameterization: Solid line, full O Xi=1
circles: m = 9; dashed line, empty circles: m = 10; dotted
where ω (i) is the O-substrate stretch frequency
line, empty circles: m=11. O/Pd(100)
of the ith adsorbed O atom in the corresponding con-
figuration. Within the frozen substrate approxima-
tion we calculate the latter stretch mode in good
ferent LGH expansions with m =9, 10, and 11 in Table
II, for which the extracted lateral interactions exhibit a agreement to the experimental data31 as ¯hωp(2×2) =
scatter of about the estimated magnitude. The critical 50meV in the p(2×2)-O configuration and h¯ωc(2×2) =
33meV in the c(2 2)-O configuration. Compared
temperaturesshowavariationoflessthan100Kandthe ×
to the computed stretch frequency of h¯ω =
overallformoftheT (Θ)-curveisalmostuntouched. The O2(gas)
c 190meV of the O molecule, this yields an esti-
systematicsofthepresentLGHexpansionprocedurepro- 2
mated vibrational contribution to the binding energy of
videsthusanerror-controlledlinkbetweentheelectronic
structurecalculationsandthemesoscopicstatisticalsim- −(h¯ωp(2×2)−¯hωO2(gas)/2)/2=23meVand−(h¯ωc(2×2)−
¯hω /2)/2=31meV per O atom for the two config-
ulations, which allows to assess which quantities can be O2(gas)
urations, respectively. In the LGH expansion these con-
determined with which uncertainty for a given accuracy
tributions get separated into a non-coverage dependent
of the underlying first-principles energetics.
part, which enters the on-site energy, and a coverage-
dependent part, which enters the lateral interactions.
Taking the difference between the estimated p(2 2)-O
B. Uncertainties in the first-principles energetics ×
and the c(2 2)-O vibrational free energy contributions
×
asameasureofthecoveragedependencewethusarriveat
The energeticparametersinthe LGHexpansionineq. anuncertaintyof 5meVinthelateralinteractionsand
(1) comprise total and vibrational free energy contribu- an uncertainty of∼ 30meV in the on-site energy due to
tions. Uncertaintiesinourapproacharisethereforeoutof the neglect of vibr∼ational contributions in our approach.
the treatmentofthe vibrationalfree energycontribution Theon-siteenergyisthusmuchmoreaffectedbythepre-
and the approximate DFT energetics, where the latter dominantly non-coverage dependent shift in vibrational
contains uncertainties due to the numerical setup and frequencyfromthe O stretchtothe O-substratestretch
2
due to the approximate exchange-correlationfunctional. mode.
We will discuss these three sources of uncertainties sub-
sequently.
2. Numerical uncertainties in the DFT energetics
1. Vibrational free energy contribution Turning to the effect of the approximate DFT total
energies, we first distinguish between numerical inac-
Separating each of the energetic parameters in the curacies which arise out of the finite basis set, the k-
LGH (i.e Fon−site, Vp, Vt, ...) into a total energy term point sampling or the chosen supercell geometries, and
b m m
and a term due to the vibrational free energy contribu- the more fundamental uncertainty due to the employed
tion, one arrives at a LGH expansion for the total bind- approximate xc functional. In principle, the numerical
ingenergiesandaLGHexpansionforthevibrationalfree inaccuracies can be reduced to whatever limit desired,
bindingenergies. InSectionIIIBthevibrationalpartwas although this may quickly lead to unfeasible computa-
completely neglected by equating the right hand side of tions in practice. In this respect, the finite slab and vac-