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

11.1 Diperiodic Space Groups and Slab Models of Surfaces

469

alternately charged planes form the repeat unit for (111) surface of MgO (Fig. 11.2), (101) surface of TiO2 (Fig. 11.3).

Table 11.3 gives the formal plane charges in the stack of di erent surfaces in the cubic perovskite ABO3 crystal with one formula unit in the primitive unit cell. The bulk atoms distribution over atomic planes is shown on Fig. 11.4.

Table 11.3. Formal and calculated (in brackets) atomic planes charges of ABO3 cubic perovskite

Surface	(001)		(110)		(111)
Termination	AO	BO2	ABO	O2	AO3	B
AIBVO3 (KNbO3)	–1	+1	+4	–4	–5	+5
AIIBIVO3 (SrTiO3)	0(+0.27)	0(–0.27)	+4	–4	–4	+4
AIIIBIIIO3 (LaMnO3)	+1(+0.95)	–1(–0.95)	+4	–4	–3	+3

Fig. 11.4. Stacking of atomic planes for cubic perovskite ABO3 surfaces

All the surfaces in cubic perovskite structure are type-3 surfaces being stacks of alternately charged planes. It is less obvious in the case of the (001) surface in AIIBIVO3 as the repeat unit consists of neutral atomic planes (see below). It is also seen that the charge of the atomic plane depends both on the oxidation states of A and B atoms (the sum of oxidation states is in all cases 6) and Miller indexes of the surface.

For the stacking sequence of alternately charged planes (producing a dipole moment perpendicular to the surface) the two-plane repeat unit produces a potential at large distances whose magnitude is given by

a
V = 2π|q|S	(11.7)

where a is the interplanar spacing. Addition of an extra neutral repeat unit of two planes on the surface of the crystal will a ect the energy of ions an infinite distance

470 11 Surface Modeling in LCAO Calculations of Metal Oxides

below the surface. Consequently, the Madelung sum cannot be truncated and must include contributions from every plane out to the surface. The potential at any ion site never reaches a constant bulk value, the energy of a neutral pair of ions never reaches the usual cohesion energy and the surface energy is infinite [772]. The stacking of planes for rutile (100), (101) and (111) type-3 surfaces is more complicated as it can not be named as the stacking sequence of alternatively charged planes.

The classification of the surfaces by Tasker is based on the formal ionicities (oxidation states) of metal atoms. As we already noted, the oxidation state and calculated atomic charges are close only in purely ionic compounds (MgO, for example). As we have seen in Chap. 9, the transition-metal atoms charges di er essentially from the oxidation states due to the covalence part of the chemical bonding with oxygen atom.

In Table 11.3 are given in brackets the atomic planes charges for (001) surface in SrTiO3 and LaMnO3 crystals calculated for the bulk crystals with the use of Wanniertype atomic functions in the population analysis, see Chap. 9 and reference [736]. As the Ti–O bond in SrTiO3 presents a non-negligible part of covalent character the actual charges of atomic planes SrO and TiO2 are nonzero. This means that SrTiO3 (001) should be considered as a polar (type-3) surface.

Type-1 or -2 surfaces have a zero dipole moment in their repeat unit and are thus potentially stable. By contrast, polar type-3 surfaces have a diverging electrostatic surface energy due to the presence of a nonzero dipole moment not only on the outer layers, but on all the repeat units throughout the material, [768]. An electrostatic instability of type-3 surfaces results from the presence of that macroscopic dipole. Type-3 surfaces can be stabilized when the macroscopic field is removed by surface reconstruction, absorption of charged specks and so on. The modification of the surface electronic structure due to the reconstruction introduces compensating charges in the outer planes, stabilizing the surface structure. The absorption of charged specks is also a very e ective mechanism to achieve the stabilization of polar orientations. Indeed, polar orientations are generally much more reactive than cleavage planes [768].

The classification of surfaces considered above was introduced for ionic or semiionic metal oxides. In covalent solids, the creation of a surface requires cutting covalent bonds, which means that dangling bonds would be present at the surface. The saturation of dangling bonds by chemisorption is important, for example, in silicates. When a surface is cut out from the bulk, unstable Si–O radicals at the surface react readily with water to give a fully hydroxylated surface with hydrophilic character [568].

Tasker’s classification of surfaces allows some qualitative conclusions to be made about the surface stability. The quantitative calculations of the surface formation energy in slab models are considered in the next sections.

11.1.3 Singleand Periodic-slab Models of MgO and TiO2 Surfaces

As was seen above in the slab model the surface plane is supposed to be fixed by its orientation (by the set of Miller indices) relative to the bulk structure and its symmetry elements. All the atoms of the bulk primitive unit cell are distributed within one or several atomic planes with the same 2D-translation periodicity, see Figures 11.2–11.4. Let these planes of atoms form a layer, being stoichiometric and neutral as by definition it contains all the atoms of the bulk primitive unit cell. Depending on the host-crystal structure the one-layer single slab can have zero (types- 1 and type-2 surfaces) or nonzero (type-3 surfaces) dipole moment along the normal

11.1 Diperiodic Space Groups and Slab Models of Surfaces

471

to the surface plane. To ensure the zero dipole moment in the type-3 surface modeling the nonstoichiometric slabs can be used.

Figure 11.3 shows that for the type-3 (111) surface of rutile TiO2 structure O- or TiO2-terminated layer of the smallest thickness consists of 5 atomic planes: O–Ti– O–TiO2–O or TiO2–O–Ti–O–TiO2, respectively.

For the cubic perovskite ABO3 (001) surface the nonstoichiometric slabs with AO or BO2 terminations can be introduced. Figure 11.4 shows a AO-terminated slab, consisting of 7 atomic planes – 4 AO planes and 3 BO2 planes. As is seen from Fig. 11.4, the slab termination is di erent for di erent surfaces.

The symmetry of the single slab corresponds to one of 80 diperiodic (layer) space groups. Fig. 11.5 shows a 3-layer single-slab model of (001) surface of MgO crystal (each layer consists of one atomic plane). The symmetry group of this slab DG61 (P 4/mmm) belongs to a square system.

Fig. 11.5. 3-layer (3-plane) slab for MgO (001) surface

Figure 11.6 shows a 3-layer single slab of (110) rutile TiO2. In this case each layer consists of 3 atomic planes. In particular, the upper layer includes an oxygen plane (O1 atoms), a Ti2O2 plane (Ti1, Ti2, O2, O3 atoms) and an oxygen plane (O4 atoms). Oxygen atoms O5 belong to the next O–Ti2O2–O layer. The symmetry group of this slab is DG37 (P mmm).

In the multislab model (periodically repeated slab) the 3D periodicity is restored: slab of n layers (nL-slab) is supposed to be periodically repeated along the normal to the surface (let the translation vector of a single slab be c). Figures 11.7 and 11.8 show the periodic 3-layer slabs, modeling MgO (001) and TiO2 (110) surfaces, respectively. The corresponding space symmetry groups can be found from Table 11.1: G123(D41h) and G47(D21h), for MgO and TiO2 3-layer slabs, respectively.

In 3D-slab calculations the results depend not only on the number of layers in the slab (slab thickness), but also on the separation between periodically repeated

472 11 Surface Modeling in LCAO Calculations of Metal Oxides

Fig. 11.6. 3-layer (9-plane) slab for rutile TiO2 (110) surface

Fig. 11.7. Periodic slab model of (001) MgO surface

11.1 Diperiodic Space Groups and Slab Models of Surfaces

473

Fig. 11.8. Periodic slab model of rutile TiO2 (110) surface

slabs defined both by the slab thickness and the length of the translation vector along the normal to the surface. The interslab distance (vacuum gap) is equal to the di erence between the length of translation vector c and the vertical distance between the top atoms of one slab and the nearest repeated one, i.e. the di erence between the translation vector c and the single slab thickness (see Figures 11.7 and 11.8).

The symmetry of the 3D-slab model is given by one of the 3D space groups G (see Table 11.1) and may depend on the slab thickness, i.e. number of layers in the slab and its termination. As is seen from Fig. 11.1, for MgO crystal (001) surface slabs of an odd number of atomic planes have inversion symmetry (relative to the central atomic plane) but slabs of an even number of atomic planes have no inversion symmetry. For the cubic perovskite ABO3 (001) surface the stoichiometric slabs (AO– BO2–AO–BO2–) consist of an even number of atomic planes and have no inversion symmetry. But the nonstoichiometric AOor BO2-terminated slabs have inversion symmetry relative to the central AO or BO2 planes, respectively.

Note that the space-symmetry group of a periodic slab di ers from the space group of the host crystal. As an example, we consider MgO crystal with symmetry of space group Oh5 . The single-slab model refers to a square lattice and has the symmetry of

474 11 Surface Modeling in LCAO Calculations of Metal Oxides

DG 61 (P 4/mmm) or DG 51( P 4mm) for slabs of odd and even number of planes, respectively (see Figures 11.2 and 11.5). As is seen from Table 11.1, the periodic slab model in these cases has the symmetry of D41h (G123) or C41v (G99) space groups, respectively.

One of the goals of slab-model studies is to investigate the dependence of the calculated surface energy on the dimensionality of the slab model and on the Hamiltonian applied. The surface energy Esn per unity of surface area for an n-layer slab is defined by the relation

	1
En	lim		[E(n, W ) − nEb]	(11.8)
s	= n,W →∞ 2S

where E(n, W ) is the total energy in a 3D-slab calculation of an n-layer slab with the vacuum gap W and Eb is the total energy of primitive unit cell in a separate bulk calculation. In the 2D-slab model W is infinite, i.e. E(n, ∞) is the n-layer single slab energy E(n). In (11.8) S is the surface area per repeating 2D-unit cell, the factor 1/2 accounts for the existence of two limiting surfaces. Esn is the energy per unit area required to form the surface from the bulk, and it is intrinsically a positive quantity.

The value of surface energy is very dependent not only on E(n, W ) but also on the bulk total energy Eb used in (11.8). In (11.8) total energies from 3D (E(n, W )) or 2D (E(n, ∞)) systems are used, creating a problem of equivalent accuracy in periodic codes for 3D and 2D systems. Boettger pointed out [773] that any di erence between Eb and the change in E(n, w) with slab thickness will cause the calculated surface energy to diverge linearly with n. Thus, increasing the slab thickness must sooner or later lead to unacceptable results because the bulk energy from a separate calculation will never exactly equal the slope of the slab energy versus n, [774]. To avoid the divergence problem, it was suggested [773], for each slab thickness n, pick as the bulk energy Eb the di erential increase in the slab total energy upon addition of one layer

	of material:
	En	lim		1
				2S [E(n, W ) − n (E(n, W ) − E(n − 1, W ))]		(11.9)
	s	= n,W	→∞

In this expression, Eb has been replaced by (E(n, W ) − E(n − 1, W )). This relation has the merit of using only slab-related quantities, making no reference to separately calculated bulk energies.

The surface formation-energy calculations impose two strong requirements to the accuracy of calculated energies: (1) numerical precision, and (2) computational consistency. While the first condition may be relatively easily satisfied by increasing the computational thresholds, the second one can lead to considerable problems due to inevitable approximations that are inherent to all quantum-mechanical computation techniques [679]. The surface modeling of real interest lies in a semi-infinite crystal i.e. the calculation of the surface energy limit when n and W tend to infinity. However, one could not a priori guarantee the convergence to definite limit in (11.8). So, the existence of such a limit is the obligatory condition for the possibility to obtain the correct value of Es using slab models.

Use of di erent surface models (2Dand 3D-slabs) and basis sets in HF (LCAO) and DFT (PW) calculations makes it di cult to understand the origin of di erences in calculated results that may be due both to the model applied and to the features of the calculation scheme chosen (basis set, integral approximations, Brillouin-zone sampling, and so on).

11.1 Diperiodic Space Groups and Slab Models of Surfaces

475

The first attempt to compare 2Dand 3D-slab models within the same calculation scheme was made in [775] for HF LCAO studies of the surface properties of BaTiO3 in the cubic perovskite structure. The authors of [775] concluded that results for periodic 3D-slabs are systematically a ected by the interactions among repeated images, and possibly the fictitious field imposed by periodic boundary conditions.

In the LCAO approximation (for HF, KS and hybrid HF-DFT Hamiltonians) the systems periodic in 2 (2D-slabs) and 3 dimensions (3D-crystals) can be treated on an equal footing. The LCAO package CRYSTAL [23] gives the unique possibility to compare results of 2Dand 3D-slab models when the LCAO basis set is used. Such a comparison was made in [776] for the (001) surface of MgO crystal and in [777] for the (110) surface of rutile TiO2.

We discuss here the results for MgO crystal. In [776] the experimental value 4.21 ˚A of the f cc lattice constant for the bulk crystal was adopted in slab calculations. As was already noted, the square unit cell of one layer is stoichiometric as it consists of one Mg and one O atoms lying in one (001) plane (see Fig. 11.5). The calculations were made for single slabs of one to twenty atomic layers in a truncated-bulk geometry without relaxing the positions of atoms. The converged surface energy value was obtained for the single slab of 5 layers. The dependence of 3D-slab calculations results on the vacuum gap was investigated via calculating the 5-layer 3D-slabs with gaps up to 42 ˚A. The symmetry of slabs depends both on the 2D or 3D periodicity and on the evenness of the number n of atomic layers. For 3-D slabs the tetragonal space groups G 123 (P 4/mmm) and G 99 (P 4mm) with tetragonal lattice constant a = a/2 = 2.977 ˚A were used for odd and even n, respectively. The translation vector c was varied in the calculations.

Table 11.4 shows the convergence of results for the single slab (2D) MgO (001) surface model depending on the slab thickness.

Table 11.4. Convergence of results for the MgO (001) single slabs depending on the 2D slab thickness (Es – surface energy in J/m2, EF – Fermi energy in H), [776]

Number
of		HF		DFT (default ABS)		DFT (extra ABS)
layers,
n	Es		EF	Es	EF	Es	EF
1	1.4571		–0.3293	1.5516	–0.1035	1.5506	–0.1075
2	1.4689		–0.3397	1.5766	–0.1082	1.5344	–0.1226
3	1.4678		–0.3403	1.5341	–0.1217	1.5335	–0.1217
4	1.4678		–0.3407	1.5630	–0.1140	1.5339	–0.1208
5	1.4678		–0.3409	1.5341	–0.1204	1.5336	–0.1203
6	1.4678		–0.3410	1.5554	–0.1157	1.5337	–0.1199
9	1.4678		–0.3413	1.5342	–0.1193	1.5337	–0.1193
12	1.4678		–0.3414	1.5460	–0.1171	1.5338	–0.1190
15	1.4678		–0.3415	1.5345	–0.1188	1.5338	–0.1189
20	1.4678		–0.3415	1.5417	–0.1175	1.5338	–0.1187

476 11 Surface Modeling in LCAO Calculations of Metal Oxides

To eliminate the results dependence on Brillouin-zone sampling the dense Monk- horst–Pack k-point mesh was used. For primitive unit cells the 12×12×12 and 12×12 special-point sets have been taken for bulk and slab calculations, respectively. In the case of 3D-slabs the number of points in the third k-direction depends on the chosen value of the c translation vector in direct space. The latter was chosen to provide a similar size in all directions of the corresponding cyclic model of the crystal (the crystal is composed of equidistant supercells). The increasing of the 3D unit cell in direct space for producing the 2D supercell was accompanied by the corresponding reduction in the k-points mesh in the reciprocal space.

The same all-electron basis sets have been used both in HF and DFT calculations to allow the direct comparison of obtained energies. The choice of high precision tolerances for the Gaussian overlap criteria in truncation of Coulomb and exchange series ensured the high numerical accuracy required for the evaluation of the corresponding surface energy. In DFT calculations the exchange and correlation were treated in PWGGA approximation. In the CRYSTAL code the periodic DFT version in a basis of local Gaussian-type functions is implemented to avoid the numerical integration. This auxiliary basis set (ABS) is used for fitting the exchange-correlation potential. In DFT computations two types of ABS were used: (1) ABS1–even tempered basis set of 12 s-type Gaussian-type functions, and (2) the basis proposed with extra p, d, f , and g-functions ABS-2.

Table 11.4 shows the fast convergence of 2D periodic results for surface energy with increasing number of layers. Three to five MgO layers are su cient to obtain the converged values of the HF surface energy for unrelaxed slabs, in agreement with the results of other calculations, [778]. The convergence of the Fermi-level energy (highest-occupied one-electron energy level for the slab) is much slower. The limiting value seems to be reached just at n = 20. As seen from Table 11.4, the use of extended ABS2 (for the fitting of the exchange-correlation potential) proved to be necessary for obtaining the convergence of the surface energies. The DFT values of surface energy exhibit a similar fast convergence with the slab thickness as HF values, and the di erence between HF and DFT surface energy does not exceed 5%. However, the initial part of Esn for n = 3 has a di erent sign of slope for these two methods. The DFT result, which corresponds to decreasing of Esn with n, seems to be more physically reliable.

There is no large discrepancy between HF and DFT atomic charges, although the latter exhibit a slightly more covalent character [776].

The numerical value of Fermi energy EF may be used for approximate estimation (according to the Koopmans’ theorem) of the surface ionization potential that defines the adsorption energies. Due to the correlation e ects the DFT Fermi energy EF di ers essentially from that in HF calculations. Table 11.4 shows that EF is inside the valence band of a perfect crystal, i.e. the resonance surface states are predicted both in HF and DFT single-slab calculations.

Taking into account the convergence of results obtained for 2D slabs, only the 5-layer 3D models have been calculated in [776]. The results are given in Table 11.5. Table 11.5 shows that the 3D-slab with vacuum gap 6–8 ˚A gives the converged (001) surface energy of ionic MgO crystal. The DFT values of Esn exhibit a similar fast convergence as HF ones and the di erence between them is small. The convergence of the surface energy with the single-slab thickness is slower for rutile TiO2 due to

11.1 Diperiodic Space Groups and Slab Models of Surfaces

477

Table 11.5. Convergence of the results for the periodically repeated 3D MgO (001) unrelaxed slab depending on the vacuum gap for 5-layer models (primitive 2D cell), [776], (Es – surface energy in J/m2, EF – Fermi energy in a.u.)

		˚	k-point	HF		DFT
Cell dimensions, A
Gap	C		set	Es	EF	Es	EF
∞	∞		12 12	1.468	–0.3409	1.532	–0.1202
2.105	10.525		12 12 4	5.365	–0.2458	4.8937	–0.0486
4.210	12.630		12 12 3	1.503	–0.3181	1.555	–0.1164
6.315	14.735		12 12 3	1.468	–0.3215	1.533	–0.1170
8.420	16.840		12 12 2	1.468	–0.3239	1.533	–0.1174
12.630	21.050		12 12 2	1.468	–0.3273	1.533	–0.1180
21.050	29.470		12 12 1	1.468	–0.3312	1.533	–0.1186
42.100	50.520		12 12 1	1.468	–0.3352	1.532	–0.1193

the partly covalent Ti–O interaction. Table 11.6 demonstrates the convergence of the surface energy in all–electron HF LCAO single-slab calculations [777] of the (110) surface of rutile TiO2.

Table 11.6. Convergence of the results for single-slab (2D) model of (110) rutile TiO2 surface, [777], (Es – surface energy, EF – Fermi energy)

No. of atomic	No. of Ti	Total energy	Es	EF
planes (layers)	planes	(a.u)	(J/m2)	(a.u)
Perfect crystal	–	–1996.8016	–	–0.3077
9	3	–5990.211	2.193	–0.4171
15	5	–9983.815	2.183	–0.4235
21	7	–13977.418	2.183	–0.4237
27	9	–17971.021	2.183	–0.4237

It follows from Table 11.6 that surface energy is su ciently well reproduced for the 5 layers (15 atomic planes) slab, in accordance with the results of HF all–electron calculations [779]).

As in the case of MgO crystal the convergence of the Fermi energy with slab thickness is much slower than that of surface energy (see the last column of Table

11.6). Calculations of 3D-slabs for rutile TiO2 [777] show that the periodic supercells
√		√
in the plane with extensions 2c × 2a	2	and 2c × a	2	give the same surface energies
˚		˚
for the vacuum gap values 41 A and 15 A , respectively.

The LCAO approximation allows comparison of the electron-charge distribution in the bulk crystal and on the surface. The application of Wannier functions for this comparison can be found in [781], where the Boys localization criteria was applied and only the valence-band states were included for WFs generation (for bulk crystals such an approach is discussed in Chap. 9). As was seen in Chap. 9 for bulk crystals, the population analysis using Wannier-type atomic orbitals (WTAO) gives more adequate

478 11 Surface Modeling in LCAO Calculations of Metal Oxides

results for charge distribution than traditional Mulliken population analysis or the projection techniques applied in PW calculations.

The WTAOs for systems with 2D periodicity were generated in [780] for the first time and used for the surface charge distribution analysis in a single-slab model of (001)MgO and (110)TiO2. First, WFs of the bulk crystal, corresponding to valence bands, were used as a tool allowing one to estimate the values for the slab parameters needed to adequately model the surface electronic structure. Secondly, a minimal valence basis of WTAOs, which are the Wannier functions, generated from occupied and vacant Bloch states and centered on atoms, were used for the analysis of chemical bonding at the surfaces under study. HF LCAO calculations were performed for 3- layer single slabs of (001)MgO (3 atomic planes) and rutile (110) TiO2 (9 atomic planes), see Figures 11.5 and 11.6.

For construction of the WTAOs for the minimal valence basis, the corresponding energy bands are to be chosen. The s- and p-bands of oxygen atoms form the higher valence bands in both slabs considered. In the case of the TiO2 slab, the d-bands of the titanium atoms are the lower conduction bands, while the s-bands of the metal atoms in both slabs are located among vacant states high in energy. The latter makes the contribution of s-WTAOs of the metal atoms into the covalence of these atoms negligible.

One can estimate the slab thickness needed for surface modeling by a preliminary qualitative analysis of the behavior of WFs, corresponding to valence bands of a bulk crystal. Indeed, the o diagonal elements of the density matrix of a bulk crystal fall o with distance as the Wannier functions of valence bands or faster (exponentially in the case of insulators), see Chap. 4. Thus, beyond a su ciently large domain of the crystal (basic domain) the values of Wannier functions, assigned to the central unit cell of the domain, are negligible. Consequently, truncation of the rest of the crystal and imposition of the periodical boundary conditions on the domain-edge atoms, practically do not a ect the electron-density distribution inside the domain chosen. The basic domain of the crystal corresponds to a special set of k-points, which provides a convergence of the calculation results relative to the extension of this set.

A similar ideology can be used for the slab model. Let us consider the slab as consisting of the finite number of atomic planes. Let, at the surface planes of the slab chosen, the numerical values of bulk valence-band WFs be close to zero. The corresponding number of atomic planes determines the minimal thickness of a slab for the studied surface, which can be expected to provide the convergence of the results relative to increasing its thickness. However, the convergence might not occur for the chosen thickness of the slab due to an absence of periodical boundary conditions at the slab surface. But, if the slab WFs, localized in its central part, are close to the bulk ones and the localization of the slab WFs, centered near the slab edge, is to a certain extent the same as in the bulk crystal, the electron density in the middle part of the slab do not di er much from the bulk electron density, and would not change when the thickness of the slab is increased. Thus, the analysis of bulk WFs can help in estimating the slab thickness, which allows the bulk-like results in the central planes of the slab to be obtained and thus corresponds to more or less adequate surface modeling. Still, this is only a qualitative and estimative analysis.

Let us illustrate the aforesaid by the results obtained for MgO and TiO2 crystals. In [780] it is shown that Wannier functions for the valence bands of the bulk MgO