Effect of Doping on Rutile TiO2 Surface Stability and Crystal Shapes.

Transition-metal-(TM-)doped TiO2 has been considered as promising electrode material for the oxygen evolution reaction (OER). OER activity is expected to depend on the coordination of the surface atoms. In this study, we theoretically investigate the stability of low-index surfaces of TM-doped rutile, (110), (100), (101) and (001), with 50 % of the Ti atoms substituted by Sc, Y, V, Nb or Ta. For Sc and Y, we also consider models with O vacancies providing the most stable oxidation state of Sc and Y. Surface energies are calculated with DFT(+U). Based on the Gibbs-Wulff theorem, the shape of the single crystals is predicted. It is observed that p-doping leads to spontaneous oxygen loss and O vacancies cause surface reconstruction. The Wulff shapes of n-doped TiO2 have smaller contributions of the (110) facet and, for Nb and Ta, larger contributions of other facets. Given the higher coordinative unsaturation of the TM atoms in the latter, a higher catalytic activity is expected.

Introduction

Search for clean, affordable, reliable and also sustainable energy is a critical technological challenge of the 21st century. Hydrogen production by chemical water splitting, which includes the oxygen evolution reaction (OER) and hydrogen evolution reaction (HER), is a promising alternative to fossil fuels.[ , ]

IrO2 is regarded as the state‐of‐the‐art electrocatalyst for the OER. Besides IrO2,[ , ] other choices for an electrocatalyst are, for example, RuO2,[ , ] mixtures of both oxides or the pure metals. However, due to the scarcity and high price of Ir and Ru, other alternatives are desirable.

Titanium dioxide is a prominent object of research because of its numerous applications varying from gas sensing,[ , , , ] biosensors, solar cells[ , , , ] to optical interference applications and electrocatalytic reactions. TiO2 is used in photocatalysis,[ , , ] photocatalytic water treatment[ , ] and for selective organic oxidation reactions. One possibility to enhance the catalytic activity without increasing the cost of the catalyst in comparison to IrO2 and RuO2 is metal doping with less precious metals. Mixtures of IrO2 and TiO2 have been already synthesized as well as doped rutile catalysts.

Rutile is the most stable polymorph of single‐crystal and microcrystalline TiO2[ , ] and has therefore been chosen in the present study. Since stoichiometric TiO2 is a semiconductor,[ , ] p‐ or n‐doping with aliovalent transition metals can be used to increase its electronic conductivity. Sc, Nb and Ta were identified in previous work as possible dopants in rutile with negative segregation energies. For comparison with Sc and Nb/Ta, we also considered Y and V from the fifth and fourth period, respectively. Since these elements have three (five) valence electrons, they decrease (increase) the number of electrons in the valence band, corresponding to p‐type (n‐type) doping. For p‐type doping, additional models with oxygen vacancies were taken into account to study the most stable oxidation state for the dopants, namely +III for Sc and Y.

Aliovalent doping may lead to electronic ground states with open shells. Therefore, for all systems, diamagnetic (DM) and ferromagnetic (FM) states were considered. Ti substitution leads to defect states that may cross the Fermi level so that the doped systems have a conducting ground state. Experimentally, highest conductivity was observed for 25–33 % Nb/Ti substitution, which was also confirmed by theoretical calculations. In our models, 50 % of the Ti atoms of rutile are substituted by M=Sc, Y, V, Nb, or Ta, so that the primitive bulk unit cells have the composition MTiO4. This is above the optimal substitution range reported in the experiments, but significantly reduces the number of possible cation configurations and thus the computational effort.

The energies of the low‐index surfaces of MTiO4 were calculated at density functional theory (DFT) level and with the Hubbard U correction (DFT+U) using a plane‐wave basis and the PBE−D4[ , , ] method as described in more detail below. From the calculated surface energies, the macroscopic crystal shapes were predicted using the Gibbs‐Wulff theorem, and compared to undoped rutile.

Results and Discussion

Bulk Optimizations

In our previous work, we identified Sc, Y, V, Nb and Ta as promising candidates for Ti substitution in rutile to obtain stable mixed oxides with sufficiently high conductivity. Experimentally, it was observed that 25–33 % Nb/Ti substitution leads to the highest electrical conductivity. In order to reduce the number of possible cation configurations in our models, we replaced every second Ti atom by the previously mentioned transition metals M and obtained primitive unit cells with MTiO4 stoichiometry, compare Figure 1, or, for the primitive unit cells, Figure S1 in the Supporting Information. Even though methods for modeling random cation distributions exist, for example, special quasirandom structures, we decided to apply symmetry restrictions to our surface models in order to avoid artificial dipole moments. These unphysical dipole moments result from non‐symmetrical cation arrangements in Tasker type 3 surfaces. The dipole correction implemented in VASP is not effective enough for a non‐symmetrical approach, compare Table 3 and S1 in the Supporting Information. Our approach strongly reduces the number of cation configurations and must be therefore considered a major simplification. It is, however, assumed that the most important effects of Ti substitution are taken into account. The 50 % substituted TiO2‐rutile bulk models were fully optimized with PBE−D4[ , , ] (see Computational Methods). The resulting lattice parameters a, c, their ratio c/a and the angle γ are given in Table 1. The calculated parameters of pure rutile overestimate the experimental low‐temperature values by less than 0.7 %. Ti substitution by M in rutile leads to a lowering of symmetry. The geometry optimizations were performed without symmetry constraints. Therefore, the optimized lattice parameters do not correspond to the tetragonal crystal class anymore, as can be seen from the angle γ. 50 % doped rutile solids are not yet experimentally available, so no direct comparison to experimental values is possible. However, according to Vegard's law the lattice parameters of the mixed oxides are expected to range between those of pure TiO2‐rutile and pure metal oxides MO2. The measured lattice parameters of VO2 at 360 K are smaller than those of rutile. This is also reflected in the calculated lattice parameters of VTiO4 and TiO2. Due to their larger atomic radii, the unit cells MTiO4 with M=Nb, Ta, Sc and Y have larger lattice parameters than TiO2‐rutile. This is also in agreement with the experimental lattice parameters of NbO2. The lattice parameters of NbTiO4 and TaTiO4 are very similar. YTiO4 has the largest lattice parameters, 0.3 Å and 0.2 Å larger than a and c of TiO2, respectively. The angle γ becomes slightly smaller than 90° for VTiO4, NbTiO4 and TaTiO4 and slightly larger for ScTiO4 and YTiO4. Using PBE−D4+U increases the lattice parameters, for both TiO2 and VTiO4. The lattice parameters increase with increasing U value, yielding larger deviations from experiment, compare Table 1. For pure rutile, the deviations increase from 0.6 % to 1.2 % for a and from 0.1 % to 2.2 % for c, when U is increased from 0 eV to 5 eV. Therefore, a small U value (2 eV) is used for titanium in VTiO4. But even with this small value, the lattice parameters of VTiO4 change drastically. a decreases so much that it even underestimates the lattice parameter of pure VO2, although an average value of the lattice parameters of TiO2 and VO2 is expected according to Vegard's law. c becomes larger than c of TiO2. Consequently, the c/a ratio increases from 0.629 up to 0.669. Due to this drastic change of the lattice parameters, for the following surface calculations U(Ti)=2 eV is used for TiO2 but not for VTiO4, while U(V) is set to 1 eV, minimizing the structural changes. These values are similar to previous estimates for the optimal U values, using experimental oxidation enthalpies of metal oxides as reference, U(V)=1.0 eV and U(Ti)=2.5 eV, respectively. Other studies with GGA+U also showed an increase in the lattice parameters of TiO2 with larger U.[ , , , , ] Other work applying SCAN+U, with U=2.5 eV, yielded lattice parameters of 4.62 Å and 2.99 Å for TiO2, which are very similar to our findings. One explanation for the increase in lattice parameters is the enhanced localization of d‐orbitals when the U correction is applied, leading to reduced overlap. The c/a ratio of TiO2 is relatively independent of the U value. Despite the overestimation of the lattice parameters, the structure obtained with U(Ti)=2 eV has the same ratio as the experimental structure. In a previous work the c/a ratio of TiO2 increased with U and is associated with a weaker Ti−Ti interaction. For VTiO4, all PBE−D4+U calculations result in too large c/a ratios compared to the experimental TiO2 ratio. With PBE−D4, the ratio is between the experimental values of pure TiO2 and VO2, as expected. We therefore conclude that the addition of a Hubbard U deteriorates the accuracy of calculated structure parameters for MTiO4.

Figure 1

Surface models of doped TiO2‐rutile with Ti as blue spheres, M as green spheres and O as red spheres; visualized with VESTA.

Table 3

Converged surface energies of the low‐index rutile surfaces of MTiO4 with M=Ti, V, Nb, Ta, Sc, Y in J/m2. In order to avoid artificial dipole moments, supercells (SC) with an inversion center/reflection plane were constructed for the surfaces (100) and (001). VASP‐PBE−D4 results, Hubbard U in eV.

System	(110)	(100)	(100) SC	(101)	(001)	(001) SC
TiO2	0.47	0.79	–	1.06	1.21	–
−U(Ti)=2.0	0.64	0.94	–	1.20	1.45	–
VTiO4	0.39	0.74	0.52	0.92	1.22	1.11
−U(V)=1.0	0.35	0.71	0.54	0.88	1.23	1.11
NbTiO4	0.60	0.96	0.86	1.10	1.16	1.24
TaTiO4	0.59	1.11	0.86	1.27	1.21	1.37
ScTiO4	0.73	0.93	0.63	0.67	1.19	1.33
YTiO4	0.88	0.92	0.20	0.23	0.78	1.15

Table 1

Calculated and measured (Exp.) lattice parameters a and c of the undoped and doped rutile structures, their ratio c/a and angle γ. DM denotes a diamagnetic state, FM denotes a ferromagnetic state.

	U(Ti) [eV]	U(M) [eV]	a [Å]	c [Å]	γ [ °]	c/a
TiO2 DM	–	–	4.616	2.957	90	0.641
	2.0	–	4.623	2.976	90	0.644
	2.5	–	4.625	2.981	90	0.645
	5.0	–	4.634	3.003	90	0.648
	7.0	–	4.643	3.020	90	0.650
Exp. TiO2 [a]			4.58666(4)	2.95407(3)	90	0.644
VTiO4 FM	–	–	4.601	2.892	89.5	0.629
FM	–	1.0	4.547	2.965	89.6	0.652
FM	2.0	1.0	4.546	2.981	89.6	0.669
FM	–	3.4	4.543	2.983	89.6	0.657
FM	2.0	3.4	4.547	2.993	89.6	0.658
Exp. VO2 [b]			4.5546(3)	2.8514(2)	90.0	0.626
NbTiO4 DM			4.749	2.982	88.9	0.628
Exp. NbO2 [c]			4.8463(1)	3.0315(1)	90.0	0.626
TaTiO4 DM			4.748	2.999	88.9	0.632
ScTiO4 FM			4.741	3.080	91.3	0.650
YTiO4 FM			4.917	3.222	92.3	0.655

[a] Ref. [45]; [b] Ref. [47]; [c] Ref. [49].

–

–

–

–

–

–

–

–

–

–

Ti substitution by Sc and Y may lead to the formation of oxygen vacancies due to the instability of the +IV oxidation state of M. In order to take this possibility into account, we constructed supercells and removed one oxygen atom, yielding M2Ti2O7 with M=Sc, Y. Two vacancy sites were considered, (a) between two Ti atoms and one M atom, and (b) between two M atoms and one Ti atom, see Figure S2 in the Supporting Information. The bulk optimizations were performed with CRYSTAL and PBE−D3(BJ), since it allows to maintain the initial bulk symmetry even if the symmetry of the atom positions in the cell is lowered. The use of D3(BJ) is slightly inconsistent with the VASP‐PBE−D4 surface calculations, however, the D4 correction is not yet implemented in CRYSTAL. Since the effect of using different dispersion corrections on the optimized lattice parameters is not pronounced (maximum 0.2 % for Sc2Ti2O7 and 0.3 % for Y2Ti2O7 in the bulk calculations), we do not expect significant deviations between PBE−D3(BJ) and PBE−D4. The influence of the position of the oxygen vacancy on the optimized lattice parameters is obvious from Table 2. For both M=Sc and Y substitutions, structure (a) is more stable than (b). This configuration results in negative vacancy formation energies and therefore predicts spontaneous oxygen loss of the stoichiometric bulk. The optimized bulk structure of (a) is used as basis of the following surface calculations. The vacancy formation energy of the p‐type doped rutile is calculated via Eq. 1. Here, MTiO4 is the perfect bulk, M2Ti2O7 is the bulk with oxygen vacancies and O2 is molecular oxygen in its triplet state.

Table 2

Calculated lattice parameters a, c, their ratio c/a and vacancy formation energy (Evac) of diamagnetic oxygen defective (O‐vac.) M2Ti2O7 calculated with CRYSTAL using PBE−D3(BJ) and pob‐TZVP basis sets. The γ angle was fixed to 90°.

	O‐vac. between	a [Å]	c [Å]	c/a	Evac [eV]
Sc2Ti2O7	Ti‐Ti‐Sc	4.615	3.207	0.695	−0.70
	Sc‐Sc‐Ti	4.711	3.089	0.656	−0.17
Y2Ti2O7	Ti‐Ti‐Sc	4.770	3.354	0.703	−1.25
	Sc‐Sc‐Ti	4.890	3.184	0.651	0.20

All bulk models were calculated with diamagnetic (DM) und ferromagnetic (FM) states. For TiO2 and NbTiO4, all atomic magnetizations are zero. For TaTiO4, the energy difference between DM and FM states is only about 0.03 eV and the difference in lattice parameter only 0.2 %, therefore we performed surface calculation for the DM state in order to reduce computational effort. The magnetizations of the structures with FM ground states are displayed in Table S2 in the Supporting Information. For VTiO4, the magnetizations of Ti and O are small, but for V it is close to one, as expected for the standard +V oxidation state of vanadium. For ScTiO4 and YTiO4, the metal magnetization is close to zero. In contrast, the spin density of oxygen is relatively high, especially for the surfaces, compare Table S2 in the Supporting Information. This is consistent with the calculated instability of bulk ScTiO4 and YTiO4 (see above). Consequently, models with oxygen vacancies were constructed to check if Sc and Y in their standard oxidation state +III lead to more stable models.

Surface Models

For the construction of the surface models with oxygen vacancies, the optimized CRYSTAL‐PBE−D3(BJ) bulk structures were used in order to maintain the tetragonal symmetry. All other surfaces were constructed from VASP‐PBE−D4 bulk geometries. Test calculations for NbTiO4 using VASP‐PBE−D3(BJ)[ , ] showed that the difference between VASP‐PBE−D3(BJ) and VASP‐PBE−D4 optimized bulk lattice parameters is smaller than 0.2 %. and supercells were constructed for (100) and (001) surfaces, respectively. In these models it was possible to arrange the M and Ti cations in a symmetric way, so that no dipole moment in z direction resulted which may lead to artifacts in the calculation of slab energies. The surface energies E for each surface model are calculated according to Eq. 2. Here, E Slab (n) is the energy of the relaxed surface model with n stoichiometric layers, is the energy of the optimized bulk structure and A is the area of the surface unit cell.

The resulting energies are given in Table 3. The order of surface stability of TiO2‐rutile obtained in earlier theoretical studies was verified by the present calculations. With VASP‐PBE−D4 the order of stability is (110) > (100) > (101) > (001) as well. Surface energies are also similar, exhibiting energy differences of 0.26 J/m2 for (001) and less than 0.15 J/m2 for the other surfaces, compared to our previous study at hybrid DFT level. In agreement with earlier work, the surface energies were found to increase with increasing coordinative unsaturation of Ti and O surface atoms. Convergence of within <0.1 J/m2 was obtained with 6 layers for (110), (100), (101) and with 10 layers for (001) slab models. By applying U(Ti)=2.0 eV in TiO2, the surface energies increase by 0.17, 0.15, 0.14 and 0.24 J/m2 for (110), (100), (101) and (001), respectively. However, the order of stability remains the same. Previous theoretical studies also observed an increase of the surface energy of TiO2(110) with increasing values of U.

–

–

–

–

For the surfaces (100) and (001), supercells (SC) with an inversion center or reflection plane, respectively, were constructed in order to avoid dipole moments induced by the Ti/M substitution. In the supercells each stoichiometric layer contains one Ti and one M atom and thus is dipole‐free. This is not the case in the surface models based on the primitive cell. By comparison of the (100)/(001) primitive and supercell results in Table 3, it becomes apparent that the dipole moment of the primitive surface unit cells leads to significant changes in . Consequently, only the results from dipole‐free slab calculations were taken into account.

The most stable M/Ti distribution for the (110) surface was investigated using VASP‐PBE‐D3(BJ). For M=Nb and Ta, the surfaces with 5‐fold coordinated (5c) Ti and the 6‐fold coordinated (6c) transition metal M are more stable than the opposite configuration. The latter is more stable for Sc, Y and V when substituting 50 % of titanium. For Nb‐ and Ta‐doped compounds the surface energies follow the same stability order as undoped TiO2 with similar ratios. VTiO4 also follows the same trend, but the energy difference between (100) and (110) is smaller. Applying U(V)=1.0 eV in VTiO4 changes the surface energies by ±0.4 J/m2, but the order of stability remains the same. In contrast, ScTiO4 and YTiO4 have a completely different order of stability. In particular, for Y‐substituted rutile, the (100) and (101) surfaces are significantly more stable than (110).

For the oxygen vacancies in Sc2Ti2O7 and Y2Ti2O7, three different models were constructed. The O vacancies are either in the outermost surface layers, or in the first and second layer, or in the central layers of the slab. All surfaces have 8 layers with 8 doped metal atoms, resulting in four oxygen atoms being removed from each surface, compare Figures S3‐S6 in the Supporting Information. During geometry optimization all surfaces reconstructed. The surface energies of O‐deficient Sc2Ti2O7 and Y2Ti2O7 differ drastically from the energies of the models without oxygen vacancies, see Tables 3 and 4. Due to the surface reconstruction, some surface energies are even negative, which makes the application of the Gibbs‐Wulff theorem impossible. The order of stability depends on the position of the oxygen vacancies. For Sc2Ti2O7, the (110), (100) and (101) surfaces are more stable when the O vacancy is located directly at the surface, in contrast to (001), where the most stable surface has the oxygen vacancies in the first and second layer. Y2Ti2O7(100) and (101) follow the same trend as Sc2Ti2O7, whereas (110) has the reverse order. For (001), the vacancies are most stabilized when located in the middle of the slab. This can be partly explained with the spin densities of oxygen, see Table 2 in the Supporting Information. For both ScTiO4 and YTiO4, the oxygen magnetization on surfaces (110) and (101) is high for the oxygen atoms at the surface. For (100) and (001), high O magnetization is distributed among oxygen atoms at the surface and lower lying atoms, for (001) even in the third oxygen atom layer. As mentioned above, this high magnetization is consistent with the instability of ScTiO4 and YTiO4 with respect to oxygen loss, especially for the (110) and (101) surfaces. Thus, only M2Ti2O7 models with oxygen vacancies with Sc and Y in their standard oxidation state of +III should be taken into account and the ScTiO4 and YTiO4 surfaces should be seen as artifacts.

Table 4

Converged surface energies of the low‐index surfaces of Sc2Ti2O7 and Y2Ti2O7 with oxygen vacancies (O‐vac.) in different surface layers in J/m2. VASP‐PBE−D4 results.

	O‐vac. in	(110)	(100)	(101)	(001)
Sc2Ti2O7	1st	0.05	−0.31	−0.10	0.96
	1st & 2nd	0.25	0.14	0.29	0.51
	Mid	0.56	0.65	1.30	0.64
Y2Ti2O7	1st	0.25	−0.42	−0.23	0.07
	1st & 2nd	0.22	−0.11	0.01	0.31
	Mid	−0.60	0.08	0.84	−0.04

In the following, we will discuss bond lengths. The metal‐oxygen bond lengths of pure TiO2 displayed in Table 5 change by −0.16 to −0.05 Å from bulk to surfaces. When applying U(Ti)=2.0 eV, the bond length increase by up to 0.03 Å. The Ti−O lengths of VTiO4 only show minor changes of ±0.04 Å when using U(V)=1.0 eV. This is surprising considering the observed changes of the lattice parameters. For NbTiO4 and TaTiO4, the M−O and Ti−O bond lengths are similar. ScTiO4 and YTiO4 exhibit considerably longer M−O bonds than Ti−O.

Table 5

Calculated M−O bond lengths of the nearest neighbors of the undoped and doped rutile structures of the bulk and low‐index surfaces after optimization in Å. For the doped structures both Ti−O and M−O are listed. DM denotes a diamagnetic state, FM denotes ferromagnetic states.

	M	Bulk	(110)	(100)	(101)	(001)
TiO2 DM		1.95	1.79–1.88	1.83–1.90	1.81–1.88	1.78–1.84
−U(Ti)=2.0		1.96	1.82–1.91	1.85–1.93	1.83–1.91	1.80–1.87
VTiO4 FM	Ti	1.94	1.81–1.92	1.79–1.88	1.79–1.88	1.77–1.91
	V	1.92	1.78–1.87	1.74–1.90	1.67–1.90	1.77–1.95
−U(V)=1.0	Ti	1.96	1.79–1.89	1.77–1.92	1.80–1.87	1.77–1.88
	V	1.93	1.80–1.90	1.73–1.91	1.69–1.92	1.78–1.92
NbTiO4 DM	Ti	1.98	1.83–1.99	1.82–1.97	1.82–1.91	1.80–2.01
	Nb	2.03	1.91–2.04	1.83–1.99	1.83–2.01	1.89–2.04
TaTiO4 DM	Ti	2.00	1.85–2.02	1.83–1.96	1.88–1.97	1.85–2.01
	Ta	2.02	1.92–2.02	1.85–1.98	1.89–1.97	1.88–2.03
ScTiO4 FM	Ti	1.94	1.76–1.91	1.79–1.90	1.80–1.92	1.78–1.87
	Sc	2.09	2.03–2.08	1.97–2.07	1.94–2.07	1.95–2.08
YTiO4 FM	Ti	1.93	1.75–1.88	1.78–1.89	1.80–1.90	1.75–1.87
	Y	2.23	2.17–2.22	2.13–2.20	2.14–2.21	2.03–2.22

In order to predict the equilibrium shapes of single crystals of the oxides, the Gibbs‐Wulff theorem was applied using the open‐source program VESTA and the calculated surface energies. The Wulff shape of pure TiO2 (Figure 2) is a rectangular prism consisting of 77 % (110) and a smaller amount of (001) surface planes with truncated corners of (101). This is in good agreement with experimental results regarding the overall crystal shape as well as the percentage of the (110) surface.[ , , ] The relative contribution of all surfaces is given in Table 6. The corners consist of approximately 22 % (101), whereas (100) is not present. The rutile particles investigated by Ohno et al. have a similar shape compared to the present results. Other works suggested dominating (110) and (100) surfaces or (110) and (111). Theoretical Wulff constructions were presented by Ramamoorthy et al. showing appreciable contributions of the (100) surface but no (001) facets. Novell‐Leruth et al. report (110) as the dominating facet and some (101) contributions.

Figure 2

Wulff constructions of (a) TiO2, (b) TiO2 with U(Ti)=2.0 eV, (c) NbTiO4, (d) TaTiO4 (e) VTiO4, (f) VTiO4 with U(V)=1.0 eV, (g) ScTiO4 and (h) YTiO4; visualized with VESTA.

Table 6

Contribution of the low‐index surfaces to the crystal shape according to the Wulff scheme in % (VESTA results).

System	(110)	(100)	(101)	(001)
TiO2	77	0	22	<1
−U(Ti)=2.0	73	0	27	0
VTiO4	72	6	22	0
−U(V)=1.0	79	0	21	0
NbTiO4	72	0	28	<1
TaTiO4	76	0	21	3
ScTiO4	15	35	50	0
YTiO4	0	55	45	0

The Wulff shapes of n‐doped NbTiO4 and TaTiO4 resemble that of rutile TiO2. Approximately 72(76)% of the overall surface area comprises (110) terminations, 28(21)% are (101) and less than 1(3)% comprises (001), for the Nb‐(Ta‐)doped compound. When V‐doped, the rectangular prism of rutile has additional truncated edges consisting of 6 % (100) and the (001) surface contributions disappear, but still bares similarities to pure TiO2.

The Hubbard U correction changes the Wulff shapes of TiO2 and VTiO4 only slightly, as illustrated by Figure 2. The small contributions of surfaces (001) and (100) vanish for both TiO2 and VTiO4. Changes of the surface contributions are smaller than 8 %, compare Table 6.

In contrast, p‐doping with Sc or Y without oxygen vacancies completely changes the predicted crystal shape of rutile. As (110) is not the most stable surface anymore, the (100) and (101) facets become much more pronounced. In the order VTiO4, ScTiO4, YTiO4 the crystal shapes become more cubic capped by square pyramids on top and bottom. The (110) surface facets vanish completely for YTiO4 and the surfaces (100) and (101) have almost the same percentage, compare Table 6, whereas for ScTiO4, (101) is the most pronounced surface with approximately 50 %. With 15 %, the (110) surface is not the most relevant surface anymore. These changes have to be seen as artifacts because high spin densities of oxygen in MTiO4 surface models indicate oxygen loss. This is also evident from the O vacancy formation energy in the surfaces, compare Table 4, and in the bulk, compare Table 2. For the reconstructed Sc2Ti2O7 and Y2Ti2O7 surfaces with oxygen vacancies, calculated surface energies are negative, which prevents the application of the Gibbs‐Wulff theorem.

The rutile (110) surface has alternating 5‐fold (5c) and 6‐fold coordinated (6c) metal atoms at the surface. The oxygen atoms also expose two different coordinations, namely 2c and 3c. Thus, this surface has undercoordinated 5c‐Ti and 2c‐O. The (100) and (101) surfaces also feature 5c metal atoms and 2c and 3c oxygen atoms, but lack 6‐fold coordinated metal atoms. For (101), the 3c oxygen atoms are not on the surface and for (100), they are below the 2c‐O. The (100) surfaces exposes MO5‐square pyramids that are tilted towards the surface normal. Due to the tilted MO5‐square pyramids the topmost oxygen atoms of the surface are geometrically blocked, which could hinder adsorption processes. On the (001) surface, only 4c metal sites and 2c oxygen atoms are present.

Gunasooriya and Nørskov state that significant contributions to the OER activity come from less stable surfaces compared to the most dominant one. In their work, the surfaces IrO2(101) and (100) are more active than the most stable surface (110). Experimentally, for the state‐of‐the‐art OER catalysts IrO2 and RuO2 the (100) surface is more active than (110), which can be partially attributed to high atomic densities of coordinatively unsaturated sites of Ir/Ru, which are electrochemically active. Other experiments state that the (001) and (101) surfaces of RuO2 are more active than (110). It has been suggested that the OER activity correlates with the number of active undercoordinated metal sites. For RuO2‐rutile, the most active surfaces are (100) and (101) with higher density of undercoordinated Ru compared to (110) and (111). Also theoretical calculations of the overpotential for RuO2 report a higher activity for (001) compared to (110). The least stable surface, (001), has highly unsaturated metal atoms and should therefore be the most reactive surface from the described ones due to the excess of bonding sites. (100) as well as (101) should also be more catalytically active compared to (110). Taking this into account, all transition metal doped rutile compounds should be more reactive than pure rutile TiO2 because of the smaller contribution of the stable (110) surface. According to these calculations, the least stable but most reactive (001) surface is only present for NbTiO4 and TaTiO4 of the doped systems.

Conclusion

We analyzed the surface stabilities of pure and 50 % doped rutile TiO2 with substituents M=V, Nb, Ta, Sc, Y which form solid solutions with rutile to MTiO4. For pure TiO2‐rutile, (110) is the most stable surface, therefore the most relevant and also most studied for catalysis and surface science. The predicted crystal shape is a rectangular prism with truncated corners, similar to that observed experimentally. Substitution with aliovalent transition metals change the order of surface stability and thus the crystal shape in comparison to pure TiO2. Nb‐ and Ta‐doped compounds resemble TiO2 with only minor changes, but in favor of (001). For VTiO4, the (100) surface becomes more prominent which should affect the catalytic activity. The higher percentage of surfaces with more unsaturated metal atoms is expected to increase the catalytic activity of the n‐doped compounds compared to pure rutile TiO2. The use of the Hubbard U correction leads to larger deviations of lattice parameters in comparison to available experimental references.

p‐doping with Sc or Y facilitates formation of oxygen vacancies in ScTiO4 and YTiO4. All oxygen‐deficient Sc2Ti2O7 and Y2Ti2O7 surfaces undergo strong reconstruction. Negative surface energies indicate that the oxygen atoms on top of the surfaces are eliminated under equilibrium conditions. Due to the negative surface energies no Wulff construction is possible.

Computational Methods

All calculations, except for the bulk geometry optimizations of oxygen‐deficient Sc2Ti2O7 and Y2Ti2O7, were carried out with the plane wave Vienna ab initio simulation package (VASP),[ , , , ] version 6.1.2, using the Perdew‐Burke‐Ernzerhof (PBE) functional and the projector‐augmented wave (PAW) approach.[ , ] The respective PAW parameters were extracted from the VASP POTCAR files library, Ti_sv_GW (12 valence electrons (VE)), V_sv_GW (13 VE), Nb_sv_GW (13 VE), Ta_sv_GW (13 VE), Sc_sv_GW (11 VE), Y_sv_GW (11 VE) and O_GW (6 VE). In order to take long‐range dispersion into account, the D4 correction[ , ] was applied except for test calculations with PBE−D3(BJ).[ , ] The energy cutoff was set to 900 eV. The spin state was the same in corresponding bulk and slab calculations, either diamagnetic or ferromagnetic. For spin polarized calculations with V, Sc and Y, the total magnetic moment was not fixed during calculations. Full geometry optimizations were performed for the bulk structures using 4×4×6 k‐points.

For the surfaces, the primitive unit cells of the 50 % doped bulk were cut along the different hkl planes. During surface optimization, the lattice constants were fixed at the optimized bulk values but with γ=90°. This was either approximated from the slightly distorted VASP‐PBE−D4 bulk structures, or based on CRYSTAL‐PBE−D3(BJ) (see below) bulk optimizations with symmetry restriction, when the structural distortion was too large. In the surface models, only the atomic positions were relaxed. Convergence criteria were set to 10−6 eV for the electronic self‐consistent cycle and 0.01 eV/Å for structure optimization. The Monkhorst‐Pack grids were set to 8×4×1 for (110) slabs, 6×4×1 for (100) slabs, 4×4×1 for (101) slabs and 8×8×1 for (001) slabs. Denser k‐point meshes were tested for TiO2 and the energy difference was 0.01 eV, which was considered as converged. Convergence of the calculated surface energy within <0.1 J/m2 was obtained with 6 layers for (110), (100), (101) and with 10 layers for (001) slab models. All presented calculations were performed with these numbers of layers except for the 8‐layered M2Ti2O7 surfaces. The converged vacuum distance was identified as 12 Å. The surface calculations were performed for primitive unit cells for (110) and (101), supercells had to be generated for (100) and (001) surfaces to avoid a dipole moment, namely a 2×1 supercell for (100) and a supercell for (001). In these supercells, metal layers contain both M and Ti, see Figure S1 in the Supporting Information.

For DFT+U, the implementation by Dudarev et al. with , simply called U here, with and being averaged Coulomb interactions, was used. Different Hubbard U values for the transition metals were tested, among them literature values of U=3.4 eV for V from VO2[ , , ] and U=2.0 eV for titanium from TiO2. Molecular oxygen was calculated with VASP‐PBE−D4 in a triplet state and a large non‐cubic simulation box. During the VASP bulk geometry optimizations of the oxygen‐deficient M2Ti2O7 bulk, the angle γ changed significantly. These calculations were performed with CRYSTAL17, version 1.0.2 and pob‐TZVP basis sets. Different from VASP, CRYSTAL allows symmetry‐restricted bulk optimizations, so the tetragonal symmetry could be maintained. Due to the lack of D4‐dispersion implementation, PBE−D3(BJ)[ , ] was used. The truncation criteria for bielectronic integrals were set to 10−7, 10−7, 10−7, 10−14, 10−42, the Gilat net was , the Monkhorst‐Pack net . Anderson mixing of the Fock/Kohn‐Sham matrix and the direct inversion of the iterative subspace convergence accelerator (DIIS) were applied. The positions of the oxygen defects in M2Ti2O7 are displayed in Figures S2‐S6 in the Supporting Information.

Conflict of interest

The authors declare no conflict of interest.

As a service to our authors and readers, this journal provides supporting information supplied by the authors. Such materials are peer reviewed and may be re‐organized for online delivery, but are not copy‐edited or typeset. Technical support issues arising from supporting information (other than missing files) should be addressed to the authors.

Supporting Information

Click here for additional data file.