Elastic polarizable environment cluster embedding approach for water adsorption on the α-Al2O3(0001) surface. A density functional study

Low coverage water adsorption on the α-Al2O3(0001) surface has been studied with a generalized-gradient density functional approach using embedded cluster and periodic slab methodologies.

An advanced cluster embedding method in an elastic polarizable environment (EPE), which enables an accurate description of the adsorption-induced substrate relaxation, has been applied systematically at various density functional levels: PW91, BP, and PBEN.

In addition, periodic slab model calculations based on the PW91 functional were carried out for varying surface supercell sizes, (2 × 2) and (3 × 3), which compare very well with the corresponding embedded-cluster results.

In agreement with two recent studies employing integrated MO + MO (IMOMO) embedded cluster and periodic Car–Parrinello BLYP methodologies, our calculations predict the 1,2-dissociative adsorption to be about 10 kcal mol−1 more favorable than molecular adsorption; however, at variance with the latter study, we predict 1,4-dissociative adsorption to be least favorable.

Analysis of adsorbate-induced relaxation renders the interaction energy with the unrelaxed substrate in the 1,4-dissociative case negative (unbound complex), thus rationalizing the smallest (by absolute value) interaction energy.

Our best estimates for binding energies, at the PBEN level, for molecular, 1,2-dissociative, and 1,4-dissociative adsorption are −22.5, −31.2, and −17.2 kcal mol−1, respectively.


Aluminum oxide (alumina) is a material studied extensively experimentally as well as theoretically.

This interest derives from many industrial applications of alumina, most importantly, as a catalyst and catalyst support,1 and as substrate for microelectronic applications.2

Alumina is also important for its role in the environment as one of the major earth mineral components.

In both respects, the surface chemistry of aluminum oxide, i.e. the adsorption/desorption phenomena on a microscopic scale, is of special interest.

Under natural conditions, mineral surfaces are most probably hydroxylated.3,4

Therefore, adsorption of water and the resultant modification of the surface should be addressed prior to studying adsorption of any other chemical onto the hydroxylated alumina surface.

Yet, while the interest in the surface hydroxylation of alumina is quite wide spread, its exact mechanisms, the influence of sample preparation, and its possible effect on measured properties are still not fully understood.5,6

Structure and reactivity of a clean α-Al2O3(0001) surface, prepared under ultrahigh-vacuum (UHV) conditions, differ markedly from the corresponding properties of a surface prepared in the presence of water.7

While the clean α-Al2O3(0001) surface is terminated by a single Al layer,8–10 the hydroxylated surface is best described as O-terminated with a double Al layer directly following the surface O layer.6,11

This finding suggests that the reaction of a clean single Al-terminated surface with water should be a complex multistage process, ultimately resulting in a surface transformation into a gibbsite-like structure, which can be formally viewed as substitution of the each surface Al by three H atoms.13

Various intermediate states of surface hydroxylation are also possible, depending on the sample history, as well as on the presence of defects, steps, and other irregularities, all of which render understanding the chemistry at the surface a difficult task.

Several experimental and theoretical studies5,7,12–16 demonstrated that the clean α-Al2O3(0001) surface readily reacts with water.

In single-crystal desorption measurements,12 adsorption energies of water were found to lie in the range 23–41 kcal mol−1.

Given the complexity of the hydroxylation of the α-Al2O3(0001) surface, we focus in the following on the incipient stage of this process.

Water adsorption on α-Al2O3(0001) has been previously studied by cluster,13,14 embedded cluster,14,15 and supercell periodic13 models.

Most studies clearly indicated a preference for dissociative adsorption over molecular adsorption and extensive changes in chemical bonding.

A recent embedded cluster IMOMO (“integrated molecular orbital + molecular orbital”)17 study15 agreed well with Car–Parrinello molecular dynamics (CPMD) calculations13 regarding the calculated binding energies for molecular adsorption and one of the dissociative adsorption modes.

Whereas large unit cells are required in supercell periodic calculations, embedded cluster models should be well suited to modeling water adsorption in the low coverage regime.

Water adsorption on the α-Al2O3(0001) surface can serve as a good example to assess cluster embedding by our new “elastic polarizable environment” (EPE) method18via a comparison with methods based on periodic slab models.

The goal of the current investigation is twofold: (i) to study single water molecule adsorption on the α-alumina (0001) surface with methods based on density functional theory (DFT), and (ii) to compare two modeling approaches, namely our new cluster-embedding technique and a conventional periodic slab (band structure) approach.

In this way, we intend to prepare for the computational investigation of more complex systems relevant to environmental chemistry.

We will show that the embedded cluster method gives reliable energies and structures.

In the case of 1,4-dissociative water adsorption, we obtained binding energies and structures at variance with the results of a recent CPMD study,13 but we will demonstrate that the current results are not artifacts of the EPE cluster model, because we were able to reproduce them in periodic calculations using the same exchange–correlation functional.

Computational methods

Periodic versus cluster models

Surfaces of solids and chemical reactions on them present a greater challenge to computational chemistry, compared to molecular systems, because an infinite number of atoms has to be treated in principle.

Two basic approaches are commonly used in electronic structure investigations of surface systems: cluster19,20 and periodic slab (“supercell”) models.21,22

Whereas the latter are preferred for (though not limited to) modeling of ordered materials, the former are more suitable for describing spatially isolated complexes on a particular segment of a surface, where a chemically interesting phenomenon (adsorption or defect) takes place, ideally when no interaction with other adsorbates or defects besides those under direct study has to be taken into account.

These two approaches tackle the problem from different perspectives, and, to some extent, can be considered as complementary.

Both modeling strategies have advantages and limitations.23,24

The supercell approach is formally correct because it treats the solid as a laterally infinite system.

Model artifacts can be avoided by using (i) a sufficiently large surface unit cell (to eliminate undesired lateral interactions between adsorption complexes), (ii) a sufficiently thick slab model (to avoid undesirable effects of the unphysical second surface of the slab), and (iii) a sufficiently wide vacuum gap between repeated slabs (again to avoid unphysical interactions).

As a consequence of such highly accurate modeling, supercell methods can become computationally rather demanding, in particular in the limit of low coverage, because more atoms have to be included in a unit cell and, in many currently available programs, a very large number of plane-wave basis functions has to be used.

Cluster model approaches are more flexible and convenient in that periodicity is not required and one can cut a model cluster from a surface according to the specific properties of the chemically interesting site; in addition, all conventional quantum-chemical techniques used for molecules can be straightforwardly applied.25

These models, however, suffer from boundary effects, i.e. errors arising from the artificial and arbitrary termination of the infinite substrate.

Also, the cluster size often has a significant impact on the reliability of the calculation, whereas finite computational resources in many cases limit calculations to medium size clusters.19

Various schemes have been developed to reduce boundary effects in cluster models, typically by using embedding procedures,18,26–32 which incorporate at various levels corrections for otherwise neglected interactions with the remaining infinite solid.

A good embedding strategy (as will be demonstrated in this paper) can be a promising alternative to periodic supercell strategies, at least for systems with strongly ionic substrates.

(On the other hand, cluster embedding for metals is not of comparable practical importance because it poses rather challenging formal problems.33)

In general, models for ionic solids are hybrids of a cluster treated quantum-mechanically (QM) and its environment represented classically.

Most cluster embedding studies to date used the simplest way of cluster embedding in rigid arrays of point charges and pseudopotential ions.

The EPE embedded cluster method

Recently, we suggested and implemented an advanced tool for cluster embedding in an elastic polarizable environment (EPE), suitable for modeling strongly ionic substrates.18

We also presented an extension of this strategy for solids with covalent polar bonding (covEPE), e.g. zeolites.32

In the EPE method, the surface site under investigation (i.e. an adsorption complex) is described by a quantum mechanical (QM) cluster embedded in a crystalline substrate, which in turn is described by a classical force field as an elastic polarizable environment, represented by a slab model.18

The EPE approach affords a variational, mutually self-consistent description of the interaction of the QM cluster and its environment.

The environment is allowed to act on the QM cluster via short-range forces at the cluster boundary and the long-range electrostatic (Madelung) field.

QM cluster and its environment are allowed to polarize each other.

Both can undergo a relaxation due to an adsorbate, and the structure of the entire system is determined through an iterative minimization of the total energy given byEtot = Ecl + Eint + Elatwhere Ecl, Eint, and Elat are the energies of intracluster interactions, the coupling of the QM cluster and the environment, and the intralattice (environment) interactions, respectively.

The coupling energy Eint is defined in such a way that the cluster environment will not undergo any relaxation in the absence of an adsorbate; thus, an artificial distortion of the environment at the QM cluster boundary is avoided.

Details of the EPE cluster-embedding method have been described elsewhere.18

For touchstone complexes of metal atoms on MgO(001)18,34 and α-Al2O3(0001)35,36 and CO adsorption on α-Al2O3(0001),37 the results of our EPE-embedded cluster model and those of periodic calculations agreed well.

A particular advantage of this advanced embedding scheme is its ability to model adsorption on strongly relaxed metal-oxide surfaces, such as α-Al2O3(0001).

Embedding studies using methodologies of similar spirit have recently been carried out to model MgO31,38 and ZnO39,40 surfaces.

EPE embedded cluster calculations

The QM cluster model employed in this work (Fig. 1) consists of a core part (the chemically important region) built of 8 Al and 12 O all-electron atoms and a boundary region given by bare ion pseudopotentials41 Al3+ without any basis functions.

Termination of a QM cluster by bare cation pseudopotentials avoids an artificial polarization of the cluster electron density by leakage into the outer region which otherwise would occur due to interaction with point charges representing cations of the surrounding.42

The environment was treated classically by a force field using a shell-model description combined with the Mott–Littleton approach.43,44

The external electrostatic field of the regular slab was represented by about 500 point charges at lattice positions around the QM cluster (the closest 240 of which were allowed to relax), whereas the total Madelung field of the remaining crystal was described via a surface charge representation of the external embedding potential.30

Following our previous strategy for modeling an α-Al2O3(0001) substrate,37 we first optimized the surface slab model at the classical atomistic level using empirical parameters.45

We employed a nine-layer slab model terminated on either side by a single Al layer.37

During this optimization, the lateral unit cell parameters were kept fixed at the calculated bulk values.

Then, an embedded cluster calculation of the unperturbed surface (i.e. without any adsorbate) was carried out.

At this step, only the QM cluster was relaxed, while the environment was kept fixed at the lattice positions determined in the preceding step.

Subsequently, the energy of the resulting structure served as reference of a clean surface.

Finally, we carried out QM/EPE calculations on the active region modified by the presence of the adsorbate.

The quantum-chemical part of the calculation was accomplished using the linear combination of Gaussian-type orbitals fitting-functions density functional method46 (LCGTO-FF-DF) as implemented in the parallel code ParaGauss.47,48

We used three gradient-corrected exchange–correlation density functionals (GGA): BP,49,50 PW91,51 and PBEN.52

The Kohn–Sham orbitals were represented by flexible Gaussian-type basis sets, contracted in generalized fashion.

O and H atoms of water were described by standard contracted basis sets, (9s,5p,1d)53,54 → [5s,4p,1d] and (6s,1p)55 → [4s,p], with the polarization p or d exponents of H and O set to 1.00 and 1.15, respectively; for the substrate we used the contracted basis sets O (9s,5p,1d) → [4s,3p,1d] and Al (12s,9p,2d)56 → [6s,4p,2d].

In the LCGTO-FF-DF approach, the classical Coulomb contribution to the electron-electron interaction is evaluated by representing the electron density with the help of an auxiliary Gaussian-type basis set.46

For every atom, the auxiliary basis was constructed as follows: the exponents of s and r2 type fitting functions were generated by scaling (all or selected) s and p orbital exponents, respectively, using a standard procedure;46 “polarization exponents” were chosen as geometric series with factors 2.5, starting with 0.1 for p and 0.2 for d exponents.

The resulting auxiliary basis sets were of the type (12s,9r2,5p,5d) for Al atoms, (9s,5r2,5p) for O atoms of the surface, (9s,5r2,5p,5d) for O atoms of water, and (6s,1r2,5p) for H atoms.

First-principle periodic calculations

First principle periodic calculations were carried out with the plane-wave based Vienna ab initio simulation package (VASP)57 using the GGA PW91 exchange–correlation density functional.51

The interaction between atomic cores and electrons was described by the projector augmentation wave (PAW) method.58,59

Brillouin-zone integrations were carried out with the help of a single k-vector (Γ point) using a generalized Gaussian smearing method60 with the standard smearing width of 0.01 eV.

An energy cut-off of 400 eV was adopted throughout.

The Al2O3 substrate was modeled using an Al-layer terminated slab of nine layers, Al–O–Al–Al–O–Al–Al–O–Al, the same model as used in the EPE calculations.

The (2 × 2) or (3 × 3) unit cells contained 24 Al and 36 O or 54 Al and 81 O centers, respectively.

The slabs were repeated periodically, with a vacuum spacing of 9 Å.

The slab surfaces were first relaxed without adsorbates.

In case of the clean surface as well as under perturbation by a H2O adsorbate, the “top” four layers of the substrate were allowed to relax fully while the rest of the slab was kept fixed at the bulk-terminated geometry obtained from a periodic PW91 optimization.

Pertinent parameters of the optimized rhombohedral primitive unit cell of the bulk structure (space group Rc) were a = b = 4.80 Å and c = 13.09 Å; the spacing between Al and O layers was 0.84 Å and that between adjacent Al layers 0.49 Å.

The geometry relaxation was stopped when the total energy difference between two ionic relaxations was less than 10−3 eV and the force on each unconstrained atom was less than 0.01 eV Å−1.

For comparison, we mention that the CPMD study13 was carried out for slab models of 9 layers, separated by a vacuum gap 14 Å, a 3 × 3 surface unit cell, and an energy cut-off of 952 eV.

Results and discussion

The clean α-Al2O3 (0001) surface

The (0001) surface of α-Al2O3 is well characterized and therefore often serves as a benchmark for developing theoretical models.

The termination by a single layer of Al3+ cations has been predicted to be most stable in several atomistic61,62, and first-principles63,64 calculations and also confirmed experimentally.8–10

Relaxation was shown to be crucial for the formation of this most stable surface of α-Al2O3.

Both experiment8,9 and calculations13,64–71 demonstrated strong relaxation, by 50–90%, of the distance between the topmost aluminum and oxygen layers.

In both computational strategies, periodic slab and embedded cluster (EPE) modeling, we started by relaxing a clean surface terminated with a single Al-layer.

In the EPE model, this was done by optimizing a periodic slab model with a shell model (atomistic) force field.

In agreement with previous studies,13,68–71, a significant relaxation of the surface layers was obtained (Table 1).

Recent pseudopotential as well as all-electron slab calculations68–71 on slabs with 12–15 atomic layers, using a variety of exchange–correlation approximations (including hybrid functionals), predicted a surface relaxation behavior that is especially large for the outermost Al layer, but comprises several subsurface layers.

The interlayer distance between the first Al layer and the next O layer, Al1–O1, is substantially reduced compared to the bulk-terminated value, 0.84 Å.72

Our result from atomistic calculations, 0.29 Å, compares favorably with experimental estimates of 0.3 ± 0.1 Å8 and 0.41 ± 0.05 Å.9

Our force-field based value is also close to results of other cluster or embedded cluster calculations,14,15 which yielded 0.25 Å.

Our PW91-slab calculations either with (2 × 2) or (3 × 3) supercell, similarly to the periodic CPMD study13 where the BLYP49,73 functional had been used, predicted an even stronger decrease of the Al1–O1 interlayer distance to 0.15 Å.

Typical relaxation magnitudes obtained from periodic DFT-based calculations reported in the literature range from −79 to −94% (Table 1).

The CPMD-BLYP result13 of −82% agrees quantitatively with our PW91-slab results.

Periodic DFT-based methods significantly overestimate the experimentally determined relaxation (Table 1).

However, there is much uncertainty about these experimental values associated with difficulties in the preparation of a clean surface.

After the QM cluster was embedded and optimized in the next step of the EPE procedure, its equilibrium geometry changed somewhat compared to that optimized in periodic atomistic calculations.

This is a consequence of the fact that the force-field and the QM descriptions do not match quite well enough.18,37

Table 2 summarizes optimized interatomic distances from atomistic calculations, EPE embedded cluster model (which employed the PW91 DF for its QM part) and those from periodic slab calculations (also based on the PW91 DF).

Not unexpectedly, the geometrical parameters of the QM cluster are closer to those rendered by PW91-slab calculations than to the geometry from atomistic calculations.

The average Al1–O distance (1.72 Å) in our cluster model is 0.04 Å larger than that in the relaxed PW91-slab structure.

The interlayer spacing of the “top” Al layer to the O layer underneath, 0.23 Å, lies in between the values from periodic atomistic and DFT-based calculations.

This result corresponds to a smaller degree of surface relaxation, i.e. −72%, and is very close to other cluster model values (−70%)14,15.

Water adsorption on α-Al2O3(0001)

Adsorption of a single water molecule onto the α-Al2O3(0001) surface has previously been studied with several computational approaches.13–15

Wittbrodt et al14. carried out HF and hybrid-DF (B3LYP) calculations on medium size cluster models; they used a 6-31+G(d) basis set in the central region of the cluster and a rather small 3-21G basis in the boundary region.

These authors applied simple embedding in a field of rigid point charges to selected structures, mainly to examine the resulting changes in the energetics, which they found to be small.

They concluded that the interaction energies for molecular and “1,2-dissociative” water adsorption should roughly be the same, whereas the interaction energy for “1,4-dissociative” adsorption was calculated 3–4 kcal mol−1 smaller.

The notations “1,2-” and “1,4-” refer to different proton acceptor positions Os relative to the Al1 atom, which accepts the OH moiety; see Fig. 1.

Subsequently, this conclusion was dismissed in periodic CPMD-BLYP calculations of Hass et al.,13 who found dissociative adsorption to be favored by 10 kcal mol−1.

The authors of the latter supercell study attributed the discrepancy to deficiencies of the cluster approach.

To cast aside the differences in methods and basis sets, they repeated the non-embedded cluster calculations of ref. 14 at their CPMD-BLYP level and found 1,2-dissociative adsorption to be only 3 kcal mol−1 more favorable.

They also pointed to significant adsorption-induced long-range restructuring of the surface, which is not accounted for in small cluster models.

Indeed, most atoms in those cluster models had been constrained to their bulk positions to avoid unphysical relaxation at the cluster boundaries.

The more recent cluster study of Shapovalov and Truong15 employed the same Al8O12 model as Wittbrodt et al.,14 but at a higher level of theory, namely a two-layer IMOMO methodology.17

The “model system” treated at the coupled-cluster singles and doubles (CCSD) level was a non-embedded Al4O6 cluster with the adsorbate, whereas the “real system” treated at the MP2 level consisted of an Al8O12 cluster, surrounded by full ion pseudopotentials (without orbital basis functions) and embedded in an array of point charges.

A basis set of double-ζ quality with effective core potentials on Al was employed for both “real” and “model” systems.15

Also, a more elaborated embedding scheme, compared to the earlier cluster study,14 was used, still with a rigid cluster environment.

This improved model rendered binding energies for molecular and 1,2-dissociative adsorption in good agreement with those of the periodic CPMD study.13

Shapovalov and Truong did not consider the 1,4-geometry,15 although it constitutes the most interesting and controversial issue, in view of the discrepancy between the results of the cluster and the slab model, which have been attributed to long-range surface modifications not reflected by medium-size cluster models.13

In the present work, a somewhat different cluster model of the same size Al8O12 was chosen (Fig. 1).

This cluster allowed us to consider both 1,2-and 1,4-dissociative water adsorption using the same model.

The above mentioned cluster model of Wittbrodt et al. was criticized in ref. 13 especially in reference to the 1,4-dissociative path, due to limitations of the model size, i.e. under-coordination of the proton acceptor (labeled as Os4 in Fig. 2) and lack of surface Al centers neighboring the Os4 center, which suggests edge artifacts.

Our cluster model overcomes these drawbacks and takes into account long-range electrostatic effects in a consistent manner.

In addition, all degrees of freedom were allowed to relax in our model, whereas in the work of Shapovalov and Truong15 only the position of the central Al atom of the substrate and the adsorbate were optimized and the rest of the QM cluster as well as its environment were kept frozen at the substrate geometry.

As we will discuss below, energies for molecular and 1,2-dissociative adsorption obtained in both calculations are nevertheless in close agreement.

Structural aspects

In their periodic CPMD-BLYP study, Hass et al13. considered molecular as well as 1,2-, 1,4-, and 1,6-dissociative adsorption of water.

The 1,2-and 1,4-dissociative pathways were found to be most important; therefore, only these and molecular adsorption were considered in the present work.

Water adsorption noticeably affects the substrate geometry (Table 2).

As a result of adsorption, the distance between the Al1 and O1 planes increased in the EPE calculation, as was also observed in previous theoretical studies.13–15

The relaxation of the Al1–O1 interplanar distance relative to bulk geometry decreases from −69% for the bare cluster to −49% for molecular adsorption, to −15 to −26% for dissociative adsorption.

Table 3 summarizes the corresponding relaxation (in percent) from the present as well as other periodic and cluster calculations.

The relative upward displacement of the OH accepting Al atom is very similarly predicted by periodic as well as cluster methods.

For the bare surface, all periodic DFT-based calculations predict too large a relaxation compared to experiment.

Embedded cluster models generally give a smaller relaxation of the outermost layer spacing; the percentage values range from 70 to 72%.

For molecular adsorption, the calculated relaxation from periodic slab and embedded cluster models is −48 to −54%.

For 1,2-dissociation, the cluster-model values scatter over a wider range, from −11 to −20%, than those from periodic calculations, from −13 to −15%, but all methods clearly predict a significant recovery of the Al1-O1 layer spacing in the direction of the bulk value.

Of course, such comparisons of relaxation percentages have to be done with due caution, because the quoted value depends on the range around the “perturbation site” over which the interlayer spacing is averaged.

Still, because the calculated relaxation of the O1-layer was found small (when considered), at least a qualitative comparison can be made.

The present periodic and embedded cluster results yielded larger contraction values for 1,4-dissociative adsorption than for 1,2-dissociation, from −24 to −28%, whereas the periodic CPMD study13 gave almost the same values for both dissociative modes (Table 3).

A closer look at the structures (compared in Table 2) reveals that 1,4-dissociation is indeed the only case, where the present work and the periodic CPMD study disagree notably;13 below, we will discuss this disagreement in more detail.

In Table 2, the optimized geometric parameters at the adsorption site are compared for the EPE-embedded cluster PW91 and periodic PW91 approaches employed in our study; the corresponding structures are shown in Fig. 2.

The optimized geometries within these two very different computational methodologies, cluster and slab models, agree very well; slab model geometries for (2 × 2) and (3 × 3) supercells agree within 0.01 Å.

In the EPE model, Al–Os bond lengths are predicted either 0.01–0.02 Å shorter or 0.01–0.04 Å longer.

For molecular, 1,2-, and 1,4-dissociative adsorption, the Al–Oads bond lengths differ by 0.01, 0.02 and 0.04 Å, respectively.

Molecular adsorption results in slightly elongated Al–O distances only at the adsorption site, 1.71–1.77 Å; these distances are 0.02–0.03 Å longer than in the relaxed bare QM cluster.

The periodic PW91 calculations gave significantly larger changes, from 1.68 Å for the relaxed surface to 1.73 Å in the chemisorbed complex.

The Al–Oads distance (Fig. 2) is 1.98 and 1.96 Å in cluster and slab calculations, respectively.

This is in close agreement with the value of 1.95 Å obtained in the periodic CPMD-BLYP study.13

In our embedded-cluster model for 1,2-dissociative adsorption, we did not encounter artificially elongated Al–Os2 distances of 2–2.2 Å as obtained in previous cluster calculations.14,15

Our cluster result, 1.91 Å, is identical to that from our periodic PW91 calculations, to be compared with the periodic CPMD value, 1.90 Å.

Overall, the EPE bond lengths for molecular and 1,2-dissociative adsorption agree quite well with the corresponding CPMD values.13

The geometries from our periodic PW91 calculations are in almost quantitative agreement with that former slab model study, with the exception of 1,4-dissociative adsorption, where both our EPE and periodic models agree well, but predict a structure of the adsorption complex, which differs significantly from that of the CPMD study.13

This may be a result of the different exchange–correlation functional used, PW91 here vs. BLYP.74

Hass et al13. found significant and far extending restructuring of the substrate, in particular, the top-layer Al3 center bound to Os4 descended even below the top O layer and became six-fold coordinated.

In our EPE and periodic results, we did not observe this type of restructuring, although Al3 moved 0.03 Å “down” compared to the clean substrate, i.e. the Al3-O1 spacing decreased at variance with the Al1–O1 spacing.

However, we observed another type of restructuring: Os4 is no longer bound to Al2 and thus coordinated by only two Al atoms, both in the cluster and the slab model.

This is likely a result of “bonding competition” between surface–adsorbate and intrasurface bonding.

Breaking one of the Al–O bonds possibly accounts for the fact that we calculated the energy gain of 1,4-dissociative adsorption 12–14 kcal mol−1 smaller compared to 1,2-dissociative adsorption (see below).

Energy aspects

We now turn to a discussion of the central results of our study, the interaction energies of a single water with the α-Al2O3(0001) surface for three selected adsorption modes: molecular, 1,2-dissociative, and 1,4-dissociative adsorption (Fig. 2).

The periodic CPMD slab model study of Hass et al.,13 based on the BLYP functional,49,73 found 1,2-dissociative adsorption 10 kcal mol−1 more favorable than molecular adsorption (Table 4).

As one of the main conclusions, their work rejected simple cluster models, e.g. those of ref. 14, as they failed to predict a significantly larger binding energy for dissociative than for molecular water adsorption (Table 4).

In fact, Hass et al. repeated the bare cluster model calculations of ref. 14 (i.e. without embedding) using their own CPMD-BLYP formalism.

In this way, they obtained −25.7 kcal mol−1 binding energy for molecular adsorption and −28.9 kcal mol−1 for the dissociative adsorption.

Thus, with the “free” cluster model, the preference for dissociative adsorption was notably reduced, from 10 to 3 kcal mol−1.

Similarly, without cluster embedding, close binding energies were obtained for molecular and 1,2-dissociative adsorption at the B3LYP14 and the CCSD levels;15 an MP2 cluster calculation, preferring molecular adsorption, even gave a qualitatively wrong picture.15

For both MP2 and CCSD calculations, IMOMO cluster embedding was shown to stabilize dissociative adsorption and destabilize molecular adsorption.15

When corrected for the basis set superposition error (BSSE) by the counterpoise method75 (by 4–5 kcal mol−1), the CCSD results for the embedded cluster models on molecular and 1,2-dissociative adsorption agreed almost quantitatively with those of the earlier CPMD-BLYP slab model study.13

The results for molecular adsorption, obtained in CPMD-BLYP slab model and CCSD embedded cluster calculations, were −23.3 kcal mol−1 and −23.4 kcal mol−1, respectively, and the results for 1,2-dissociative adsorption were −33.2 kcal mol−1 and −31.6 kcal mol−1, respectively.

This rather close agreement indeed is most astonishing if one considers the markedly different procedures to account for electron correlation and the completely different approaches to account for environment effects (Table 4).

To compare periodic slab models and embedded cluster models in a strict sense, it is certainly desirable to have results at hand that were obtained for the same electronic structure approach.

Thus far, this has not been the case.

Therefore, we applied our advanced embedding strategy (EPE) and periodic first principle calculations for the PW91 functional, to obtain such a benchmark.

To make this comparison more general, we also performed EPE embedded cluster calculations for two further gradient-corrected density functionals, BP and PBEN (see Computational Methods).

We estimated basis set superposition errors using the counterpoise method.75

Finally, in periodic PW91 calculations, we compared results for supercells of two sizes, (2 × 2) and (3 × 3).

As seen from Table 4, our EPE calculations give qualitatively correct energetics in all cases: molecular adsorption is 7.6–8.7 kcal mol−1 less favorable than 1,2-dissociative adsorption, depending on the functional.

Also, the three functionals consistently render the interaction for 1,4-dissociative adsorption weakest, 5.2–5.6 kcal mol−1 less (in absolute value) than molecular adsorption.

Of all density functionals, PW91 predicts the strongest interaction; BP interaction energies are 1.6–2.3 kcal mol−1 kcal mol−1 smaller (in absolute value) and PBEN energies are 4.4–5.5 kcal mol−1 smaller (Table 4).

Overall, this agreement of the interaction energies as obtained from the EPE calculations is very satisfying and the order of interaction energies follows known trends: PW91 > BP > PBEN (in absolute value).

Experience has shown that PBEN energetics are most reliable for metal–ligand interactions and such interaction energies are usually expected to be smaller than the corresponding PW91 results which tend to overestimate such interactions.74,76

The interaction energies from periodic calculations with a (3 × 3) surface unit cell and the PW91 functional are also very similar, both in trend and magnitude, to the energies from the corresponding PW91 EPE embedded cluster calculations; the cluster results are 3.5–5.7 kcal mol−1 smaller (by absolute value) than the analogous slab model results.

In fact, the interaction energies obtained with the larger (3 × 3) and smaller (2 × 2) surface unit cells bracket the EPE cluster results (Table 4).

It is remarkable that the differences between the two slab model calculations with different surface unit cells are rather substantial, 5.5–7.2 kcal mol−1, whereas the relative energetics of the three adsorption modes remain, nevertheless, essentially unaffected by the size of the surface unit cell.

Indeed, the EPE interaction energies are rather close to those of periodic (2 × 2) calculations, with differences amounting to only 1.5–2.0 kcal mol−1 (Table 4).

From this comparison it is evident that the substrate relaxation is large and extends far, as was also concluded in the CPMD-B3LYP study;13 therefore, even larger surface unit cells may be required to recover that part of the interaction energy which is masked by the artificial periodicity of the slab models.

Qualitatively, the contributions of the substrate relaxation to the total binding energies can be deduced from the binding energies computed by optimizing only adsorbate-related degrees of freedom, i.e. with the positions of substrate atoms “frozen” at the clean surface geometry.

The corresponding results, obtained for EPE embedded clusters and (2 × 2) slab models at the PW91 level, are also given in Table 4.

Only for molecular adsorption is this latter interaction energy close to the total binding energy from EPE cluster calculations, −21.7 vs. −28.0 kcal mol−1, respectively; the corresponding energies from periodic (2 × 2) slab model calculations are very similar, −20.5 (“frozen”) and −26.5 kcal mol−1, respectively.

For 1,2-dissociative adsorption, the EPE binding energy at the frozen substrate is only −1.8 kcal mol−1, and that for 1,4-dissociative adsorption is even positive, 28.2 kcal mol−1, i.e. the complex is unbound.

The latter, surprising value is supported by the corresponding result of periodic calculations, 29.6 kcal mol−1.

Unfortunately, it was not possible to converge the geometry of 1,2-dissociative adsorption with the substrate fixed at the structure without perturbation by the adsorbate; several calculations converged to the energetically favorable structure of molecular adsorption.

From these calculations one concludes that the adsorbate-induced relaxation of the substrate is largest for 1,4-dissociative adsorption, ∼50 kcal mol−1, but only ∼34 kcal mol−1 for 1, 2-dissociative mode; for molecular adsorption it is only ∼6 kcal mol−1, i.e. almost negligible in comparison with the other two adsorption modes.

Our PBEN results for molecular and 1,2-dissociative adsorption from EPE-embedded cluster models agree almost quantitatively (up to 1–2 kcal mol−1) with the CCSD embedded-cluster results of ref. 15 and the CPMD-BLYP periodic slab model results of .ref. 13

However, this tight agreement does not extend to 1,4-dissociative adsorption.

Recall that, for the PW91 functional, the current EPE and periodic slab model results for 1,4-dissociative adsorption agree equally well with each other, as do the results for the other two types of adsorption complexes, molecular adsorption and 1,2-adsorption (Table 4).

Therefore, one cannot dismiss the large discrepancy between the PBEN EPE adsorption energy for 1,4-dissociative adsorption, −17.2 kcal mol−1, and the corresponding CPMD-BLYP slab model results, −32.5 kcal mol−1, as a cluster artifact (Table 4).

Unfortunately, the structure of the latter adsorption complex was not available to us to crosscheck the energy for our two computational approaches.

Finally, we return to the IMOMO cluster calculations,15 which are the cluster models that allow treatment by high-level post-HF correlation methods.

In that sense, these IMOMO cluster results can also be taken as a set of benchmark data.

In fact, they agree best of all with our PBEN energetics, corroborating the claimed higher quality of this density functional for the determination of binding energies.

The fact that two rather diverse cluster-embedding strategies gave quite similar results, with energy differences of only 0.9 and 0.4 kcal mol−1 (Table 4), supports the reliability of these modeling approaches.


This work demonstrated the performance of our new cluster embedding technique in an elastic polarizable environment (EPE), using low coverage water adsorption on the α-alumina (0001) surface as test system.

We compared embedded cluster results to periodic supercell calculations of this work and to literature data.

To render the conclusions of this study more general, we employed three density functionals, PW91, PBEN, and BP, and we also examined the convergence of the periodic slab results with respect to the size of the surface unit cell of periodic calculations.

We found the current PW91 embedded cluster results to agree well with those of periodic PW91 calculation, in terms of both structure and energetics.

Slab models suffer from artificially imposed periodicity (as was demonstrated by comparison of 2 × 2 and 3 × 3 models), which can only be overcome by using rather large surface unit cells.

Thus, the variational EPE cluster embedding method can serve as an economic alternative to periodic slab models, which should be particularly important for disordered surfaces and charged defect sites, common in oxide systems.

Still, there is room for improvement of the current embedding model by refining the force field; the present work adopted a shell model force field from the literature.45

Despite the scatter in magnitudes of calculated interaction energies by up to 5 kcal mol−1, the relative energetics of different water adsorption complexes were predicted uniformly by all three density functionals employed.

We consider the PBEN functional to yield the most accurate adsorption energies: −22.5, −31.2, and −17.2 kcal mol−1 for molecular, 1,2-dissociative and 1,4-dissociative adsorption, respectively.

The almost quantitative agreement of the present DF results with those of a high-level post-HF correlation method, IMOMO CCSD-MP2,15 confirms their reliability for the evaluation of adsorption energies on oxide surfaces.

Finally, comparison of the current results with an earlier periodic CPMD study,13 raises new questions.

For molecular and 1,2-dissociative adsorption, the CPMD study, based on the BLYP functional, reported binding energies very similar to ours, −23.3 and −33.2 kcal mol−1; however, we were not able to reproduce the rather large binding energy, −32.5 kcal mol−1, found in that study for 1,4-dissociative adsorption.

We demonstrated by comparison with our own periodic slab calculations that this discrepancy is not an artifact of a cluster model.

Our results therefore do not support the earlier conclusion13 about preferred dissociation into 1,4-mode, but rather suggest the 1,2-dissociative pathway.

Having presented this model study, we wish to emphasize its idealized nature, as a water reaction with a partially hydroxylated alumina surface would be of truly practical interest.

This more complex case has only once been addressed theoretically,13 but clearly more work is needed to gain full insight into this interesting chemical problem.