Quantum Chemistry of Solids / 24-Surface Modeling in LCAO Calculations of Metal Oxides

11.1 Diperiodic Space Groups and Slab Models of Surfaces

479

crystal calculated along the [001] direction are practically completely localized inside a domain containing just three (001) atomic planes. Furthermore, according to calculations the di erence between these WFs and the WFs generated for the slab and centered on the plane 2 (see Fig. 11.5), is negligible. This means that the 3-plane slab is su cient for modeling the (001)-surface properties of MgO. As we demonstrated above, this consideration is, in principal, confirmed by thorough calculations of MgO slabs with a varied thickness – the surface energy and other characteristics converge at the 3-plane slab.

In the case of TiO2 the bulk WFs does not fit into the 9 atomic planes shown in Fig. 11.6. As a result, the WFs, calculated from the slab states and localized in the “central” region of the slab, are di erent from the corresponding bulk WFs. Consequently, 9 planes are not su cient for correct modeling of the TiO2 crystalline surface. This result correlates with the results given in Table 11.6.

As seen from Table 11.7A the values of atomic charges for the MgO slab correspond to the purely ionic type of chemical bonding both at the surface and inside the slab. The charges in the slab model practically coincide with the bulk ones. The results, obtained by the WTAOs method and according to the traditional Mulliken and L¨owdin schemes, are nearly the same.

Table 11.7. Atomic charges for the single MgO 3-plane (001)-surface and TiO2 9-plane (110) surface slabs and bulk crystals, [780], (NWTAO is nonorthogonal WTAO, OWTAO is orthogonal WTAO)

A) MgO crystal

Atoms

Atomic charges, |e|

NWTAOs

OWTAOs

Mulliken

L¨owdin

O1, Surface

–1.98

–1.96

–1.96

–1.83

Mg1, Surface

2.00

1.97

1.96

1.83

O2

–1.99

–1.96

–1.97

–1.82

Mg2

1.97

1.95

1.97

1.82

O Bulk

–2.00

–1.98

–1.98

–1.82

B) TiO2 crystal

Planes

Atoms

Atomic charges, |e|

NWTAOs

OWTAOs

Mulliken

L¨owdin

1

O1, surface

–1.33

–1.30

–1.13

–0.88

2

Ti1, surface

2.99

2.87

2.62

1.84

Ti2

2.86

2.79

2.54

1.63

O2,3, surface

–1.48

–1.44

–1.35

–0.86

3

O4

–1.57

–1.51

–1.38

–0.91

4

O5

–1.54

–1.45

–1.28

–0.82

5

Ti3

3.06

2.95

2.66

1.72

Ti4

3.01

2.93

2.65

1.73

O6,7

–1.50

–1.46

–1.33

–0.87

Bulk

Ti

2.76

2.70

2.66

1.72

O

–1.38

–1.35

–1.33

–0.86

480 11 Surface Modeling in LCAO Calculations of Metal Oxides

For the TiO2 slab the Mulliken and L¨owdin population analyses give contradictory results. Particularly, as noted in Table 11.7B, the values of atomic charges at the surface planes of the slab, calculated according to the Mulliken scheme, are smaller in magnitude than those inside the slab. This corresponds to a lower ionicity at the crystalline surface compared to the bulk. An opposite tendency is demonstrated in the L¨owdin analysis, where the surface charges are larger or the same as in the bulk. Besides, the absolute values of the charges obtained according to the two traditional schemes are quite di erent from each other. Table 11.7B shows that the analysis performed by the orthogonal and nonorthogonal WTAOs gives approximately the same results. They show the tendency of decreasing the ionicity level at the surface, which is similar to the results of the traditional Mulliken analysis.

Let us consider the values of atomic charges, obtained for the central planes of the slab (planes 3–7 on Fig. 11.6). Due to the symmetry of the slab the planes 6 and 7 are equivalent to the planes 4 and 3, respectively. As noted above, a nine-plane slab is not su cient to model the TiO2 surface, as the electron density at the central planes of the slab is di erent from that of the bulk crystal. However, the values of atomic charges calculated according to both the Mulliken and L¨owdin schemes practically coincide with the bulk ones (Table 11.7B). In contrast to these results, the population analysis based on WTAOs reproduces this slab-bulk misfit: the values for the charges of atoms of the inner planes in the slab and the bulk crystal di er noticeably.

When analyzing the Ti–O bond orders, which are given in Table 11.8, one can conclude the following.

The values, obtained by the two traditional schemes, can be hardly interpreted. Some of the bonds have the same orders in both schemes (Ti1–O5, Ti2–O1, Ti3–O4, Ti4–O5), while the others are essentially di erent (Ti1–O2, Ti2–O2, Ti3–O6, Ti4–O6), see Fig. 11.6. As to the results obtained by orthogonal and nonorthogonal WTAOs, they are close for all the bonds under consideration.

Every oxygen atom (excluding the atoms O1) has two neighboring titanium atoms at a distance of 3.687 a.u. and one at 3.727 a.u. (three-coordinated atoms). The O1 atoms are two-coordinated and have the neighbors only at 3.687 a.u. The titanium atoms are six-coordinated (Ti2, Ti3, Ti4) or five-coordinated (Ti1) and have four neighboring oxygen atoms at 3.687 a.u. and two or one – at 3.727 a.u., respectively. Only for some of the atoms (Ti2, O2, Ti4) do the Mulliken and L¨owdin schemes allow the larger orders for the shorter bonds (among the nearest-neighbor Ti-O bonds) to be obtained. For other atoms the values of the shorter bonds orders are approximately of the same magnitude as the longer ones (O6) or even significantly smaller (Ti3, O4, O5), which is unreal from the chemical point of view. The WTAOs method for all the atoms gives the expected tendency - the shorter is a bond, the larger is the value of the corresponding bond order. Also, WTAOs analysis gives the noticeably larger values for the insurface bonds. This agrees with the decreased values of the atomic charges at the surface relative to those at the inner part of the slab.

The results of the calculations of atomic covalences are presented in Table 11.9. Again, the Mulliken and L¨owdin methods demonstrate contradictory and obscure bulk–surface tendencies. For some atoms the values of the covalences, calculated according to the Mulliken scheme are close to the L¨owdin ones, for the other they di er significantly. Thus, these coincidences of the Mulliken and L¨owdin results can be regarded as casual. And since the nine-plane slab does not provide the bulk-like electron

11.1 Diperiodic Space Groups and Slab Models of Surfaces

481

Table 11.8. Orders of the bonds between the near-neighbour Ti and O atoms in the single TiO2 nine-plane (110)-surface slab and bulk crystal (NWTAO is nonorthogonal WTAO, OWTAO is orthogonal WTAO, d is a distance between neighbor atoms, s denotes in-surface atom), [780].

Atom	Neighbor atoms				Bond orders
	Atom	N	d	NWTAO	OWTAO	Mulliken	L¨owdin
O1	Ti2, s	2	3.687	0.54	0.55	0.80	0.79
Ti1	O2,3, s	4	3.687	0.37	0.39	0.42	0.62
	O5	1	3.727	0.17	0.25	0.62	0.63
Ti2	O1, s	2	3.687	0.54	0.55	0.80	0.79
	O4	2	3.687	0.27	0.29	0.47	0.50
	O2,3	2	3.727	0.11	0.13	0.17	0.40
O2,3	Ti1, s	2	3.687	0.37	0.39	0.42	0.62
	Ti2	1	3.727	0.11	0.13	0.17	0.40
O4	Ti2	2	3.687	0.27	0.29	0.47	0.50
	Ti3	1	3.727	0.17	0.21	0.52	0.54
O5	Ti4	2	3.687	0.30	0.32	0.52	0.54
	Ti1	1	3.727	0.17	0.25	0.62	0.63
Ti3	O6,7	4	3.687	0.31	0.32	0.34	0.54
	O4	2	3.727	0.17	0.21	0.52	0.54
Ti4	O5	4	3.687	0.30	0.32	0.52	0.54
	O6,7	2	3.727	0.22	0.23	0.34	0.53
O6,7	Ti3	2	3.687	0.31	0.32	0.34	0.54
	Ti4	1	3.727	0.22	0.23	0.34	0.53

Table 11.9. Atomic covalences for the single TiO2 nine-plane (110)-surface slab and bulk crystal, [780], (NWTAO is nonorthogonal WTAO, OWTAO is orthogonal WTAO)

Planes	Atoms		Atomic covalences
		NWTAOs		OWTAOs	Mulliken	L¨owdin
1	O1, surface	1.20		1.24	1.86	1.88
2	Ti1, surface	1.74		1.91	2.44	3.35
	Ti2	1.92		2.03	2.36	3.02
	O2,3, surface	0.95		1.02	1.28	2.01
3	O4	0.80		0.90	1.76	1.95
4	O5	0.85		0.99	1.92	2.07
5	Ti3	1.64		1.82	2.49	3.47
	Ti4	1.73		1.85	2.87	3.45
	O6,7	0.92		0.99	1.27	2.03
Bulk	Ti	2.06		2.16	2.15	3.46
	O	1.11		1.17	1.24	2.04

482 11 Surface Modeling in LCAO Calculations of Metal Oxides

density in the middle of the slab, the correspondence between the values for the covalences of the bulk and slab atoms, obtained by the traditional population analysis schemes, also cannot be considered as reliable.

The WTAOs method, in its turn, allows one to obtain the values of atomic covalences in the slab, which can be physically interpreted. First, the values calculated with the orthogonal and nonorthogonal WTAOs are alike. Secondly, the inslab atoms have a covalence di erent from the bulk one. This result can be expected due to an insu cient thickness of the nine-plane slab for the TiO2 (110) surface modeling.

And last but not least, the results of the WTAOs population analysis clearly show the increase of the covalence at the surface of the slab when compared to the inslab values. The surface atom O1 has the maximal covalence among all the oxygen atoms of the slab. This e ect is even sharper if it is remembered that this atom is only twocoordinated, while the others are three-coordinated. The large value of the covalence is also observed for the surface atoms O2/O3, but it is not as pronounced as for the atom O1. The atom Ti1, which can be regarded as a surface atom, has a smaller covalence than the atom Ti2, but taking into account that the former is only fivecoordinated, one can conclude that the e ect of higher covalence at the surface is valid for this atom as well.

The values of local characteristics of the slab electronic structure, di ering from the bulk one, may indicate the possibility of significant structural relaxation of the studied TiO2 slab. This relaxation would involve insurface and intraslab atoms, since the TiO2 nine-plane slab at its inner planes does not reproduce the bulk electronic structure. Taking into account the results of the WTAO population analysis, one can assume that the atoms would shift so that the length of the bonds would increase or decrease to compensate the excessive or deficient values of the covalence, respectively. These considerations correlate with the studies of the geometry optimization in TiO2 slabs [782].

11.2 Surface LCAO Calculations on TiO2 and SnO2

11.2.1 Cluster Models of (110) TiO2

Titanium dioxide (in rutile and anatase structures) is the most investigated crystalline system in the surface science of metal oxides. The review article [783] summarizes the results of experimental and theoretical studies of titanium dioxide (bulk and surface) made up to 2002 inclusive. The information about calculations of the surface reconstruction, surface defects and growth of metals on TiO2 is also included. The results of the later theoretical studies of rutile surfaces can be found in [784–795] and references therein. In the majority of the calculations the slab model was used for the study of periodic surface structures.

Here, we consider the molecularand cyclic-cluster models of (110) TiO2, following [770], where the connection between the slab and cluster models was considered by performing a symmetry analysis of the crystalline surface.

The cyclic-cluster (CC) model of the surface is connected to the slab-2D-supercell approach, but it is di erent due to the di erent introduction of cyclic boundary conditions (CBC). In the slab-supercell approach, these conditions are, in fact, introduced for a very large system, e.g., for the main region of the 2D periodic plane lattice, so

11.2 Surface LCAO Calculations on TiO2 and SnO2

483

that band-theory methods can be used in the slab calculations. In the CC model, the CBC are introduced for the cluster itself, so that the model turns out to be finite. The convergence of the results to the slab values must be investigated by increasing the CC. However, the CC approach o ers evident advantages: (a) the spurious boundary e ects of molecular clusters are absent; (b) the periodicity of adsorbed species or of defects is excluded because the CC is a finite system; (c) there is a one-to-one symmetry correspondence between the electronic states of the 2D CC and the band states of the corresponding slab model.

The CC model was applied in semiempirical calculations [796] of the adsorption geometry for organic adsorbates on a rutile (110) surface. The INDO parametrization scheme was modified to reproduce the experimental results for geometries of carbonand nitrogen-containing molecules. The CC representing the rutile (110) surface consists of 120 atoms. It includes two layers of 5 × 2 2D unit cells with six atoms in each 2D unit cell. It is seen from Fig. 11.6, that 2D primitive unit cells form a plane rectangular lattice with parameters c and a (a, c are the tetragonal lattice parameters of the bulk crystal).

The molecular-cluster (MC) model is extremely popular in surface simulations because of its simplicity and flexibility and the possibility of its use in the standard techniques of molecular quantum chemistry [320, 769]. It is also possible to study defects or adsorption reactions in the limit of low coverage. In the MC approach, a crystal with a surface is modeled by a finite (molecular) system consisting of atoms on the surface and of some atomic planes nearest to it. The 2D periodicity of the surface is disregarded so that the symmetry of such a model is described by one of the crystallographic point groups. In adsorption studies, the adsorbed species are included in the cluster-adsorbate systems. An MC is cut out of the slab to simulate a portion of the surface so that spurious e ects arise that are related to the limited cluster size and to the presence of the boundary. To eliminate these e ects, the dangling bonds are saturated either by hydrogen atoms or pseudoatoms or the whole cluster is placed in an external field simulating the influence of the bulk crystal [769].

In the MC model, the results critically depend on the choice of the cluster size, stoichiometry, and shape. Nevertheless, the MC model can provide useful information if the cluster is properly chosen and the influence of its size is investigated. A number of general rules were suggested in [770, 771] for the generation of an MC suitable for surface simulations based on the connection between the MC and the slab model. These rules have been intuitively applied in di erent studies, in most cases without explicitly considering the relationship between periodic and molecular models. The following rules allow artefacts of MC surface models to be minimized: (a) the clusters correctly reflect the stoichiometry of the bulk crystal and are electroneutral; (b) atoms on the crystal surface should be equivalent to those on the cluster surface except for corners and edges; (c) the average coordination of all cluster atoms is as close as possible to that of atoms on the crystal surface; (d) as many symmetry elements of the crystal surface as possible are retained in the cluster; and (e) an extension to rule

(a)is that each layer of the cluster should have the stoichiometry of the crystal.

In fact, such principles of the MC choice ensure the simulation of the crystal surface

as closely as possible and are directly connected with the 2D periodic slab model. By using 2D unit cells as basic units of the model clusters, it is possible to increase the cluster size in a systematic way and, thereby, to study the convergence behavior

484 11 Surface Modeling in LCAO Calculations of Metal Oxides

of calculated surface properties. As examples of MC model applications to surface studies using some of the above-mentioned rules, calculations of MgO, TiO2 (rutile, anatase structures), and Cr2O3 (corundum structure) surfaces and adsorption of small molecules on them can be considered, see [770] and references therein. However, it is also possible to find examples in the literature where these rules have not been applied and, in part, unphysical results were obtained.

In MgO and NaCl crystals, the MCs were chosen [797] so that they simulated the 2D supercell of the surface and included several layers of the bulk crystal. For these crystals, the stoichiometry of the chosen MC was ensured by this choice.

For more complicated systems, e.g., TiO2 structures, the stoichiometry of the MC is only ensured when several additional oxygen atoms are added. There are two principle ways to add these additional atoms. They can be placed on the regular surface around the cluster. This is denoted as type A in Figures 11.9 and 11.10. In some cases, this procedure leads to clusters that do not correspond to rules (d) and

(e). Another possibility is to place the additional atoms for saturation at irregular positions of the crystalline lattice so that rules (d) and (e) will be fulfilled. Examples of this second type (B) of cluster generation are presented in Figures 11.9 and 11.10.

If all layers of a multilayer cluster have the same stoichiometry, artificial polarization, which can a ect the calculated surface properties, is reduced. To study the e ect of the distribution of saturation atoms on the calculated properties of crystalline surfaces, the rutile (110) surface has been selected for the MSINDO calculations [770]. The water adsorbtion in the molecular and dissociated form was modeled on clusters of Ti5O10, Ti9O18, Ti14O28, and Ti18O36 (Fig. 11.9) chosen according to the rules given above. Nevertheless, even if these rules are applied, several di erent types of clusters can be constructed. Two types of possible cluster models were used. In type- A clusters additional oxygen atoms, necessary to ensure total stoichiometry (shaded in Fig. 11.9), were placed at regular lattice positions. In type-B clusters additional oxygen atoms were placed at nonlattice positions to maximize the symmetry of the clusters and to reduce polarization, as discussed in the previous section.

Table 11.10. MSINDO adsorption energies (kJ/mol) for water adsorption on rutile (110) (relaxed cluster calculations), [770]

Cluster		Type A			Type B
	molecular		dissociative	molecular		dissociative
Ti5O10	51		39	70		201
Ti9O18	139		–24	94		107
Ti14O28	66		8	92		143
Ti18O36	157		192	129		162
Ti49O98	120		151	119		157

The adsorption energies calculated with these small cluster models are presented in Table 11.10. The geometries of the clusters were optimized within the symmetry of the rutile structure. In model B, there are one or two more degrees of freedom for oxygen atoms in nonlattice positions.

11.2 Surface LCAO Calculations on TiO2 and SnO2

485

Fig. 11.9. Rutile (110) surface of size 1 ×1: 1–4 layer clusters of types A and B, [770]. Black, white and shaded spheres label Ti, O and additional O atoms.

486 11 Surface Modeling in LCAO Calculations of Metal Oxides

Fig. 11.10. Rutile (110) surfaces of size 3 × 1, 3 × 3 and 5 × 3: two-layer clusters of types A and B, [770]. The labels of atoms are the same as in Fig. 11.9.

For the cluster–water systems, the Cartesian coordinates of all water atoms, the two surface atoms close to the oxygen atom, and the hydrogen atom of water are optimized. The adsorption energy is calculated as the di erence in the total energies of the TiO2 cluster and the cluster+H2O system. Therefore, positive values indicate stabilization. A comparison is made with the results of the calculations on a much larger cluster, Ti49O98 (Fig. 11.10), which is a better representation of the surface.

From Table 11.10, it can be seen that the di erent distribution of additional oxygen atoms in clusters A and B has a substantial influence on the calculated adsorption energies. The variation of adsorption energy with an increasing number of layers is significantly reduced if model B is considered. Even for rather small systems, the adsorption energies are relatively close to that of the largest cluster, Ti49O98. The most important di erence between the two models is that only for model B is the dissociative adsorption always more stable than that of the molecular form. The relative stability of the two forms of water on the rutile (110) surface is still a matter of debate (see the next section). At the moment, we only focus on the convergence of results obtained for clusters with increasing sizes.

A comparison of models A and B of the largest cluster (Ti49O98) shows that the influence of additional atoms is negligible due to their large distance from the adsorption position. This indicates that, in this case, it is the description of the local environment near the adsorption site rather than the e ect of the global polarization of the cluster that is responsible for the di erences observed for the smaller clusters. Thus, the numerical results for the adsorption behavior of rutile (110) towards water with di erent clusters demonstrate that convergence to the methodological limit is significantly improved if the criteria mentioned above are taken into account.

Since the underlying considerations are of a general nature, they can also be applied to cluster models of surfaces in other systems. In any case, only the comparison of the results of both slab and cluster models application allow those results of calcu-

11.2 Surface LCAO Calculations on TiO2 and SnO2

487

lations that can be taken as confident to be extracted. In the next section we consider this point returning to the water adsorption on the (110) surface of rutile.

11.2.2 Adsorption of Water on the TiO2 (Rutile) (110) Surface: Comparison of Periodic LCAO-PW and Embedded-cluster LCAO Calculations

The adsorption of water on TiO2 surfaces has been extensively investigated using both experimental and theoretical methods, see [790] and references therein. Nevertheless, the adsorption of H2O on the TiO2 (110) surface is still a matter of controversy. From experiments, it has been proposed that H2O adsorbs mainly associatively and dissociates at defect sites, see [798] and references therein. If dissociation does occur, it is only at low coverages (< 15%) that may be associated with surface defects. In contrast, most DFT 3D periodic calculations with the plane-wave (PW) basis predict dissociation at all coverages or an equivalent amount of dissociative and associative mechanisms. In contrast, Hartree–Fock (HF) embedded-cluster calculations with an atomic (Gaussian) basis set [801] predict that the associative mechanism should be favored due to overestimation of H-bonding in the dissociated configuration by DFTPW studies.

The first two-periodic all-electron HF LCAO calculations of the rutile relaxed surfaces, made in [779], gave atomic displacements of surface atoms that did not di er significantly from the later results of DFT-PW investigations. Further periodic LCAO studies of TiO2 bare surfaces have been made in [777, 799, 800]. For studies of H2O adsorption on TiO2 the single-slab periodic HF-LCAO and DFT-LCAO methods were first applied in [790] and compared with PW-DFT results to test various methods with cyclicand embedded-cluster calculations and resolve discrepancies between the methods.

In the discussion of numerical results of this study we use the following labels for atoms of the (110)TiO2 shown on Fig. 11.6: Ti1, Ti2-Ti5f ,Ti6f (fiveand sixfold coordinated titanium atoms), O1, O2-Obr, O3f (bridging and threefold coordinated oxygen atoms).

Both in 2D (single-slab) LCAO and 3D (periodic-slab) PW calculations the same 3-layer (9 atomic planes) slab was used with the fixed positions of the middle-layer atoms. Both associative (molecular) and dissociative water monolayer adsorption has been investigated.

DFT-PW calculations were performed on di erent H2O−TiO2 structures to determine which one corresponds to the most stable arrangement on the rutile (110) surface. To this end, a 3D-supercell consisting of 1 × 1 or 2 × 1 surface unit cells was

used to model the (110) surface geometry. The smallest surface unit cell was chosen
˚	√		˚
for bare-surface calculations, having dimensions of |c| (2.959 A) and		2|a| (6.497 A)

in the (001) and (–110) directions, respectively, where a and c are translation vectors for the bulk rutile unit cell. This surface unit cell is doubled in the (001) direction for the hydroxylated or hydrated surfaces. The calculations were performed for 3-layer slabs with a total cell thickness ≈ 19 ˚A, i.e. slab thickness + vacuum gap 10 ˚A. A model with the fixed atomic positions of the central layer is assumed to be more appropriate to the real 2D surface relaxation because real crystals are not thin films and the bulk crystal structure probably exists a few atomic layers beneath the mineral

488 11 Surface Modeling in LCAO Calculations of Metal Oxides

surface. Use of the constrained central layer was also accompanied by imposition of the inversion symmetry. This symmetry saves computational time and minimizes any possible dipole moment of the slab.

Adsorption was simulated with a pair of H2O molecules (one at each side of the 3-layer TiO2 slab) using a 2 × 1 unit cell. The use of one H2O molecule on each side allows a cell with inversion symmetry to be used and creates a system with zero dipole moment associated with the slab.

In a PW study of water associative adsorption (with structure optimization) it was found that the lowest energy corresponds to the following structure: the water O is directly bonded to the five-coordinate Ti5f - atom, a single H-bond between an H in H2O and a bridging oxygen atom (Obr) and the torsion angle H − O − H · · · Obr 117◦, see Fig. 11.11. In the structure with an angle of 180◦ all atoms of the water molecule and bridging oxygen are in the same plane, but this structure appears to be less favorable.

Fig. 11.11. Associative (molecular) adsorption at half-monolayer coverage

The dissociative adsorption of H2O onto the (110) surface was modeled assuming that one H atom was bonded to the Obr next to the Ti atom with a terminal OH group. (Note that this Ti atom was originally 5-coordinated on the bare (110) surface.) Due to the half-monolayer configuration, the neighboring OH groups can be adsorbed either in a zigzag geometry or inline. In the latter case, there exists the possibility of H- bonding between neighboring bridging and terminal OH groups: Obr − H · · · Oter − H. Consequently, the inline configuration only was considered because this H-bonding should lower the adsorption energy. As in the case of molecular H2O adsorption,