Introduction
This project aims to study the (0001) surface relaxation of hematite, including determining a suitable surface size and thickness for surface energy calculation. Density Functional Theory (DFT) calculations for hematite surface geometry optimization were performed using the Perdew-Burke-Ernzerh [1] of generalized gradient approximation for the exchange-correlation functional in the DMol3 code [2]. The Monkhorst -Pack k-point grid in this calculation is determined by the convergence of the free energy of bulk hematite. The slab thickness is chosen when the atomic layer in the center of the slab remains the same after DFT calculation.
Hematite Crystal structure
This project uses energy converged hematite (α-Fe2O3) structure with lattice parameter a=b=5.05Å, c=13.81Å, α=β=γ=90 degrees, and space group R-3c [3]. The equivalent fractional coordination in bulk hematite for iron is (0, 0, 0.353694), for oxygen is (0.3059, 0, 0.25). Iron occupies two-thirds of slightly distorted octahedron sites, and oxygen locates at the tetrahedron site. The Fe-O3-Fe unit repeated along the z-axis and stacked oxygens form a hexagonal close-packed (HCP) structure (Figure 1).
Figure 1. The structure of hematite, with its space group R-3c. Oxygens (red) stacked hexagonal close-packed, and irons (blue) locate at two-thirds of the octahedron site.
Construction of (0001) Surface
Former research shows that surface (0001) in hematite react to water and metal legends actively. The (0001) surface of the hematite model was built with half-infinite slabs, ten atomic layers, and 10Å vacuum. Choosing symmetric slab reduces the generated dipole between uppermost and button layer. The thickness of the slab should be large enough to ensure the surface is enough for relaxation. In this calculation, we use ten atomic layers to relax the surface atom. The distances between O-Fe and Fe-Fe atomic layers are shown in Figure 2. Due to the distorted octahedra site in hematite structure, the irons are not in the same atomic plane and thus regarded as in two different planes with their distance of 0.607Å.
Figure 2. Ten atomic layers present with (0001) cleaved surface, blue is iron and red is oxygen.
(0001) Surface relaxation result
The slab thickness can only be chosen if the atomic positions at the center of the slabs retain the same with bulk coordinates. A DFT-optimized 1*1 surface geometry optimization was calculated initially using a 5*5*1 Monkhorst-Pack k-point grind and 4.5Å global orbital cutoff [4]. The core treatment for iron and oxygen atoms is applied to all electrons. The energy convergence and SCF tolerance are 1.0e-5 Hartree eV and 1.0e-6 eV separately. The maximum force applied on each atom is 0.002 Ha/Å, and maximum displacement for each atom is 0.05 Å. Figure 3 shows iron and oxygens coordination after relaxation. The relaxed layer-1 oxygen anions have moved 2.21% to the positive z-direction from their bulk positions; the relaxed layer-2 iron has moved 1.31% to the negative z-direction from the bulk positions. The distance of first two layers and last two layers showed a contraction from 0.847 Å to 0.688 Å (Figure3). Moreover, the oxygen in the first layer moved not only along z axis but also along x and y axis to relax the surface energy (Figure4).
Figure 3. DMol3 calculation with K-point set 5*5*1
Figure 4. A comparison of first layer oxygen original coordinates (left) and after surface relaxation (right). Oxygen in the first and last layer changed all their coordinates to relax surface energy.
The energy decreased and converged as the atoms find their preferred location (Figure 5). The total energy of using 5*5*1 and 6*6*1 k-points is -8483.9002004eV and -8483.900402eV. This geometry optimization calculation shows the coordination of the center of the atom layer (layer 5-6, Fe-Fe layer) changed swiftly than in bulk after calculation, suggesting this model required larger suitable slab size for geometry optimization. But we get similar total energy and percent realization using 5*5*1 and 6*6*1 k-points, showing that the k-point is suitable for the geometry optimization (Table 1).
Figure 5. Energy convergence in geometry optimization, giving a result of -8483.900402 eV using 6*6*1 k-point
Table 1. DMol3 geometry optimization result and percent relaxation of surface (0001) with k-point 5*5*1
Conclusion
In this project, the slab was built symmetrically with a 10Å vacuum. The Brillouin zone integrated with a 5*5*1 Monkhorst-Pack k-point grid. The global orbital cutoff is set as 4.5Å. However, the slab in this calculation is not enough to relax the (0001) surface within ten atomic layers because the center layer of the atoms moved too much after DFT calculation. A good indication of surface relaxation is the atom position of the center remaining the same after relaxation while a first and last couple of layers are relaxed.
Reference
[1] Perdew, J. P., Burke, K., & Ernzerhof, M. (1996). Generalized gradient approximation made simple. Physical review letters, 77(18), 3865.
[2] Delley, B. (2000). From molecules to solids with the DMol 3 approach. The Journal of chemical physics, 113(18), 7756-7764.
[3] Trainor, T. P., Chaka, A. M., Eng, P. J., Newville, M., Waychunas, G. A., Catalano, J. G., & Brown Jr, G. E. (2004). Structure and reactivity of the hydrated hematite (0 0 0 1) surface. Surface Science, 573(2), 204-224.
[4] Lo, C. S., Tanwar, K. S., Chaka, A. M., & Trainor, T. P. (2007). Density functional theory study of the clean and hydrated hematite (1 1¯ 02) surfaces. Physical review B, 75(7), 075425.
