2020 Volume 61 Issue 1 Pages 78-87
Modelling of nonpolar surfaces is indispensable for first principles calculations to understand key surface properties, such as the surface energy and band edge positions including the ionization potential and electron affinity. An algorithm to generate nonpolar slab-and-vacuum models that require surface modification to attain stoichiometry was developed and outlined in this paper. Removal of half of the atoms in the outermost layers in a stripe pattern made possible auto-generation of many important nonpolar and stoichiometric slab-and-vacuum models including all spinel and perovskite aristotype surfaces. In addition, a computational procedure that assists the generation of step models was proposed that facilitates investigation of step edges that are often active reaction sites.
Fig. 2 Examples of generated ZnAl2O4 surface slabs. Red, blue, and grey balls indicate Zn, Al, and O, respectively. The surface orientation is (a–c) (111), (d) (100), (e) (110), and (f) (210), respectively.
Understanding the interactions of atoms and molecules at a surface is a vital issue in many areas including catalysis, photocatalysis, gas sensing, and crystal growth. Some electronic properties are surface dependent, where an important example is band edge positions that make possible prediction of the Schottky barrier heights of metal-semiconductor interfaces, band offsets of semiconductor heterointerfaces, and doping limits.1–4) First principles calculations complement experiments because the surface is described precisely at the atomic level, allowing investigation of both realistic and hypothetical situations. Slab-and-vacuum models (hereafter simply “slab models”) in three-dimensional (3D) periodic boundary conditions, where infinitely extending 2D thin films are separated from their images by finite vacuum, are commonly used to model surfaces. Automated generation of such models are highly desirable in the world of high-throughput calculations,5,6) and there are efforts to address this issue.7–10)
An unavoidable problem when discussing surfaces is polarity. A slab where the two surfaces are not equivalent has a macroscopic dipole moment perpendicular to the surface plane that diverges when the thickness is increased. This dipole moment can be cancelled by charge modification or reconstruction including desorption or adsorption of atoms, but finding the surface-dependent compensation mechanism is not trivial.11,12) Tasker’s categorization of surfaces into three distinct types13) is a well-known classification of polarity. Planes in a Tasker’s type 1 surface are neutral with both anions and cations, whereas those in a Tasker’s type 2 surface are charged and arranged symmetrically such that there is no dipole moment perpendicular to the unit cell. In contrast, a Tasker’s type 3 surface is charged with a perpendicular dipole moment and considered as a polar surface.
On the other hand, Hinuma et al. classified nonpolar slabs into three categories.9) After simply cleaving slabs from bulk, the surface is nonpolar type A if it is not polar and each layer of atoms is stoichiometric, and nonpolar type B if it is not polar and stoichiometric as a whole. A nonpolar type C surface requires modification of the surface to obtain a nonpolar and stoichiometric slab. Prototypical examples of the three types are rocksalt (100), fluorite (111), and rocksalt (111) surfaces, respectively.
This paper expands the algorithm by Hinuma et al.9,10) that automatically generates nonpolar and stoichiometric nonpolar models by simply cleaving bulk, if any exist. The procedure to obtain, if possible, a nonpolar and stoichiometric type C surface by removing half of the atoms in the outermost layers is outlined. All surfaces of the spinel (space group type Fd$\bar{3}$m, number 227) and perovskite (Pm$\bar{3}$m, number 221) aristotypes are Tasker type 3 surfaces but are actually nonpolar type C (see section 2.3). These crystal structures have a wide variety of important applications: the diverse properties of spinel structure materials regarding magnetism, optics, electricity, and catalysis are playing significant roles in applications of data storage, biotechnology, electronics, laser, sensor, conversion reaction, and energy storage/conversion.14) On the other hand, perovskite structures are used in (electro)catalysis, catalysis, superconductivity, ferroelectricity, piezoelectricity, and thermoelectricity.15,16) Automated generation of surfaces of these important crystal structures could therefore significantly accelerate research of (electro)catalysis, catalysis, superconductivity, ferroelectricity, piezoelectricity, and thermoelectricity.
In addition, a procedure to assist generation of surfaces with steps is also given. Steps are well-defined atomic size defects and can be used to tune the functionality of a material by changing its concentration through changing the inclination (vicinal) angle.17) Experimental observation of the spacing has been reported using atomic force microscopy (AFM) tomography together with high resolution transmission electron microscopy (HRTEM).18) Step edge sites are less coordinated than terrace sites and often acts as active reaction sites in catalysts.19–21) As an example in simple elementary substance metals, the {211} surface has been used as a surrogate surface to model under-coordinated step edge sites in face centred cubic (fcc) metals because this surface can be viewed to consist of {111} and {100} facets. Examples of studies employing the fcc {211} surface include investigations on H adsorption on Pd,22) O adsorption on Pt,23) and methane reforming on Ni.24) A more complicated case is steps on a compound. The c-plane of sapphire, or (0001) surface of α-Al2O3, become nucleation sites for nanostructures such as WSe2,25) carbon nanotubes,26) and nanowires of GaN27) and ZnO.28,29) The charge distribution near the surface is an additional issue in ionic compounds. The excess charge approach by Finnis,30) which counts the charge in a fixed volume near a step, has been used to evaluate the stability of steps of ionic compounds.31) Investigation of stability without expensive calculations is a useful tool but is not sufficient when evaluating the formation energy of steps. A well-defined step formation energy is obtained only by using crystallographically nonpolar models, which is a concept independent from stability. Note that slab models with steps of not only compounds but of elementary substances will be polar in the crystallographic definition if the two surfaces of the slab are not identical. Therefore, a computer assisted procedure to generate, if it is possible, nonpolar step surfaces for any system would greatly advance the study of steps on a surface, be it an elementary substance, a complex oxide, or whatever crystal that is available. Such a process is proposed in this paper.
Indices i = {1, 2, 3}, m is an integer, and n is a positive integer unless indicated otherwise. * is used for numbers that are not relevant to the discussion, and * is not necessarily the same number when there are multiple appearances. I is the identity matrix. Cartesian coordinates are denoted as (x, y, z).
2.1 Basis vectorsBasis vectors are expressed using column vectors and are those of a cell under 3D periodic boundary conditions. Any point x can be represented by a column vector of coordinate triplets as xT = (x1, x2, x3). The position of this point x in a supercell with basis vectors (a, b, c) is denoted by a column vector x, and this point is in the supercell when 0 ≤ xi < 1 for all i. A transformation matrix converts one set of basis vectors into another. For instance, a 3D transformation matrix M converts basis vectors from (a1, b1, c1) to (a2, b2, c2) as
| \begin{align} (\mathbf{a}_{2},\mathbf{b}_{2},\mathbf{c}_{2}) &= (\mathbf{a}_{1},\mathbf{b}_{1},\mathbf{c}_{1})\boldsymbol{{M}} \\ &= (\mathbf{a}_{1},\mathbf{b}_{1},\mathbf{c}_{1}) \begin{pmatrix} M_{11} & M_{12} & M_{13}\\ M_{21} & M_{22} & M_{23}\\ M_{31} & M_{32} & M_{33} \end{pmatrix} \end{align} | (1) |
An isometry [see section 1.2.1 of the International Tables of Crystallography A, 6th Edition (ITA)32)] is a mapping assigning a unique image point to each point of the space while keeping all distances invariant. An isometry can be represented using matrix formulation as
| \begin{equation} \begin{pmatrix} \tilde{x}_{1}\\ \tilde{x}_{2}\\ \tilde{x}_{3} \end{pmatrix} = \begin{pmatrix} W_{11} & W_{12} & W_{13}\\ W_{21} & W_{22} & W_{23}\\ W_{31} & W_{32} & W_{33} \end{pmatrix} \begin{pmatrix} x_{1}\\ x_{2}\\ x_{3} \end{pmatrix} + \begin{pmatrix} w_{1}\\ w_{2}\\ w_{3} \end{pmatrix} \end{equation} | (2) |
This isometry is also denoted as $\skew3\tilde{\boldsymbol{{x}}} = \boldsymbol{{Wx}} + \boldsymbol{{w}}$ or $\skew2\tilde{\boldsymbol{{x}}} = (\boldsymbol{{W}},\boldsymbol{{w}})\boldsymbol{{x}}$, where W is the matrix part and w is the column translation part.
The geometric element is the set of fixed points of a symmetry operation after removal of the intrinsic translation part (section 1.2.2.4 of the ITA). The geometric elements of inversion, two-fold rotation/screw, and mirror/glide isometries are discussed here. The coordinate triplets of a fixed point is denoted by xF. For an inversion isometry (−I, w), the fixed point is defined by (−I, w)xF = xF, or xF = w/2, and hence the geometric element is this single point. For a two-fold rotation/screw and mirror/glide isometry (W, w), the intrinsic translation part is t/2 = (1/2)(W + I)w by definition and after some algebra we obtain (W, w − t/2)xF = xF hence (W − I)(xF − w/2) = 0. An intrinsic translation part that is exactly a lattice vector results in a rotation or mirror isometry. Note that (W − I)−1 cannot be defined because (W − I) is a singular matrix and the set of xF that satisfies this relation is not a point, which is the case when (W − I) is not singular, but is the rotation/screw axis or mirror/glide plane.
Explicit forms of isometries are outlined below. Unconventional basis vector choices are frequently used in this paper, thus the most general form is given.
2.2.1 InversionThe matrix part should be the negative of the identity matrix, and there are no restrictions on the vector part. The form is
| \begin{equation} (\boldsymbol{{W}},\boldsymbol{{w}})_{\bar{1}} = \begin{pmatrix} \bar{1} & 0 & 0 & w_{1}\\ 0 & \bar{1} & 0 & w_{2}\\ 0 & 0 & \bar{1} & w_{3} \end{pmatrix} \end{equation} | (3) |
Three-fold and higher order rotation/screw symmetry are not considered because we are interested in symmetry operations that maps the crystal onto itself after applying it twice. A two-fold rotation/screw operation along the a axis has the form
| \begin{equation} (\boldsymbol{{W}},\boldsymbol{{w}})_{2x} = \begin{pmatrix} 1 & 0 & w_{13} & w_{1}\\ 0 & \bar{1} & 0 & w_{2}\\ 0 & 0 & \bar{1} & w_{3} \end{pmatrix} . \end{equation} | (4) |
| \begin{equation} (\boldsymbol{{W}},\boldsymbol{{w}})_{2x}(\boldsymbol{{W}},\boldsymbol{{w}})_{2x} = \begin{pmatrix} 1 & 0 & 0 & 2w_{1} + w_{3}w_{13}\\ 0 & 1 & 0 & 0\\ 0 & 0 & 1 & 0 \end{pmatrix} , \end{equation} | (5) |
| \begin{equation} (\boldsymbol{{W}},\boldsymbol{{w}})_{2y} = \begin{pmatrix} \bar{1} & 0 & 0 & w_{1}\\ 0 & 1 & w_{23} & w_{2}\\ 0 & 0 & \bar{1} & w_{3} \end{pmatrix} \end{equation} | (6) |
The matrix form of two-fold rotation/screw operations along a + b is
| \begin{equation} (\boldsymbol{{W}},\boldsymbol{{w}})_{2 +} = \begin{pmatrix} 0 & 1 & w_{13} & w_{1}\\ 1 & 0 & w_{13} & w_{2}\\ 0 & 0 & \bar{1} & w_{3} \end{pmatrix} . \end{equation} | (7) |
| \begin{equation} (\boldsymbol{{W}},\boldsymbol{{w}})_{2 -} = \begin{pmatrix} 0 & \bar{1} & w_{13} & w_{1}\\ \bar{1} & 0 & - w_{13} & w_{2}\\ 0 & 0 & \bar{1} & w_{3} \end{pmatrix} . \end{equation} | (8) |
The matrix form of mirror and glide operations is
| \begin{equation} (\boldsymbol{{W}},\boldsymbol{{w}})_{mz} = \begin{pmatrix} 1 & 0 & w_{13} & w_{1}\\ 0 & 1 & w_{23} & w_{2}\\ 0 & 0 & \bar{1} & w_{3} \end{pmatrix} . \end{equation} | (9) |
| \begin{equation} (\boldsymbol{{W}},\boldsymbol{{w}})_{mz}(\boldsymbol{{W}},\boldsymbol{{w}})_{mz} = \begin{pmatrix} 1 & 0 & w_{13} & 2w_{1} + w_{3}w_{13}\\ 0 & 1 & w_{23} & 2w_{2} + w_{3}w_{23}\\ 0 & 0 & \bar{1} & 0 \end{pmatrix} . \end{equation} | (10) |
The conventional unit cell as summarized in Hinuma et al.33) can be obtained using, for instance, the spglib code.34) Its basis vectors are denoted as (a, b, c), its basis vector lengths as a, b, and c, and interaxial angles as α, β, and γ. An out-of-plane vector of the (hkl) surface, cOP, is defined as cOP = ha + kb + lc, and an “in-plane vector” h′a + k′b + l′c must satisfy (h′, k′, l′) · (h, k, l)T = 0. An in-plane vector may or may not be perpendicular to the out-of-plane vector.
According to the definitions and proposed algorithms in Hinuma et al.,9) an (hkl) primitive cell, which is constructed from the conventional unit cell, is defined as a primitive cell with basis vectors (aP, bP, c1P) where aP and bP are in-plane vectors that are Gaussian reduced. An (hkl) n-supercell is defined as an 1 × 1 × n supercell of the (hkl) primitive cell and has basis vectors (aP, bP, cnP). Different (hkl) n-supercells are related by cnP = nc1P, and the (hkl) 1-supercell is the same as the (hkl) primitive cell. The (hkl) n-supercell of a polar surface does not contain an isometry of the form
| \begin{equation} (\boldsymbol{{W}},\boldsymbol{{w}}) = \begin{pmatrix} * & * & * & *\\ * & * & * & *\\ 0 & 0 & \bar{1} & 2z_{c} \end{pmatrix} . \end{equation} | (11) |
Nonpolar type A and B slab models are derived by employing a (hkl) n-supercell and removing atoms with z-coordinates outside ranges z− + m ≤ z ≤ z+ + m. The slab centre is defined as zc = (z+ + z−)/2 with 0 < z− < z+ < 1, and the (hkl) n-supercell contains an isometry of the form in eq. (11).
The proposed algorithm derives nonpolar type C slabs by taking a supercell with double in-plane area from an (hkl) n-supercell, which is hereafter denoted as the (hkl) n-double cell, removing atoms with z-coordinates outside ranges z− + m ≤ z ≤ z+ + m, and then removing half of the atoms with z-coordinates z+ + m (top layer) and z− + m (bottom layer). The two in-plane basis vectors of the (hkl) n-double cell are denoted as aD and bD, and these are related to aP and bP through two 2D transformation matrices such that
| \begin{equation} (\mathbf{a}_{\text{D}},\mathbf{b}_{\text{D}}) = (\mathbf{a}_{\text{T}},\mathbf{b}_{\text{T}})\boldsymbol{{T}}_{2} = (\mathbf{a}_{\text{P}},\mathbf{b}_{\text{P}})\boldsymbol{{T}}_{1}\boldsymbol{{T}}_{2} \end{equation} | (12) |
Nonpolar type C slabs can be obtained with this scheme if the final slab contains at least one of five in-plane isometry types: $\bar{1}$, 2y, 21y, m, and b. Isometries are prioritized in the following order when there are multiple such isometries:
$\bar{1}$ penetrated by 2y, $\bar{1}$ penetrated by 21y, other $\bar{1}$, 2y, 21y, m, and b.
Note that if 2x, 21x, or a exist, basis vectors can be transformed to convert the symmetry operation to be 2y, 21y, and b, respectively. 2y, 21y, and b are preferentially used over 2x, 21x, and a to align all principal isometry rotations along bD and/or avoid an non-zero intrinsic translation part along aD. Although the order is arbitrary, $\bar{1}$ is given preference because of ease in finding the symmetry-related site opposite the slab. 2y and 21y penetrating $\bar{1}$ is automatically preserved if such $\bar{1}$ is used as the principal isometry, hence $\bar{1}$ penetrated by 2y is given highest priority followed by $\bar{1}$ penetrated by 21y. $\bar{1}$ not penetrated by 2y nor 21y comes next, followed by 2y, 21y, m, and b. The relative order of m and b does not matter because these cannot appear simultaneously. The translation part of a principal isometry has the restriction 0 ≤ w1 < 1 and 0 ≤ w3 < 1. If there are two choices for w1 in principal isometries $\bar{1}$, 2y, and 21y, or in other words, there is one choice where 0 ≤ w1 < 0.5 and another where 0.5 ≤ w1 < 1, the former is used. Two choices appear when aD/2 is a lattice vector.
3.2 Symmetry of slabsThe symmetry of a single slab can be described with layer group symmetry with five Bravais lattices and 80 layer groups (see International Tables of Crystallography E (ITE)35) for details). However, the models that are discussed are entirely based on a 3D model under 3D periodic boundary conditions. Symmetry reduction caused by cleaving can result in specialized cell metrics (enhancement of the Euclidean normalizer, section 3.5.2 of the ITA) of the slab. Therefore, naïve application of the layer groups should be avoided.
The Bravais lattices of layer groups are primitive hexagonal (hp), primitive square (tp), primitive rectangular (op), centred rectangular (oc), and primitive oblique (mp). There is a constraint on the interaxial angle γ of hp, op, and oc conventional cells. As a result, aT · bT = 0 is imposed for cells with metrics allowed in hp, op, and oc lattices (denoted as potential hp, op, and oc regardless of the actual layer group). T1 of a potential hp cell is chosen such that det T1 = 2 and $\sqrt{3} |\textbf{a}_{\text{T}}| = |\textbf{b}_{\text{T}}|$, whereas T1 = I and |aT| = |bT| is adopted for a potential tp cell. A potential oc cell has det T1 = 2, and aT/2 + bT/2 is a basis vector. Whether the slab has potential hp, tp, or oc symmetry is determined by the flowchart in Fig. 1.
Flowchart to determine T1.
Gaussian reduction in the (hkl)-primitive cell forces aP = bP but γP can be either 60° or 120°. Tentative in-plane basis vectors are designed such that $a_{\text{T}} = \sqrt{3} b_{\text{T}} = \sqrt{3} a_{\text{P}}$, γT = 90°, and det(T1) = 2. This type of slab is obtained only when cleaved perpendicular to a three- or six-fold rotation, screw, or rotoinversion axis. The only possible orientations are {0001} of the hexagonal crystal family and {111} of the cubic crystal family. Observation of symmetry operations of space group types in the ITA reveals that relevant point groups and possible in-plane isometries are $\bar{3}$ (in-plane isometry is $\bar{1}$), 32 (two-fold rotation and screw, denoted as 2 + 21 for short), $\bar{3}$m (either $\bar{1}$ + 2 + 21, $\bar{1}$, or 2 + 21), $\bar{6}$ (m), 6/m ($\bar{1}$ + m, $\bar{1}$, or m), 622 (2 + 21), $\bar{6}$m2 (either 2 + 21 + m, 2 + 21, or m), 6/mmm (either $\bar{1}$ + 2 + 21 + m, $\bar{1}$ + 2 + 21, $\bar{1}$ + m, 2 + 21 + m, or 2 + 21), m$\bar{3}$ ($\bar{1}$), 432 (2 + 21), and m$\bar{3}$m ($\bar{1}$ + 2 + 21, $\bar{1}$, or 2 + 21). An example of a slab categorized as potential hp but is actually not hp is a space group type P3112 (number 151) crystal with {0001} orientation where the slab centre plane includes two-fold rotation and screw axes. The three-fold screw symmetry is removed when cleaving the slab and the three-fold rotation symmetry operation, which exists in all hp slabs, is not found in the resulting slab.
The following statements holds for hp slabs based on observation of the ITA: “if there is a two-fold rotation axis as an in-plane isometry, there is also a two-fold screw rotation axis as an in-plane isometry, and vice versa” and “if there is an inversion centre in the slab centre plane, it is penetrated by at least one two-fold rotation axis”. The flow to determine the principal isometry is as follows:
Tentative in-plane basis vectors are designed such that aT = bT and γT = 90°. T1 = I is used. This type of slab is obtained only when cleaved perpendicular to a four-fold rotation, screw, or rotoinversion axis. The only possible orientations are {001} of the tetragonal and cubic crystal families. Observation of symmetry operations of space group types in the ITA reveals that relevant point groups are 4/m, 422, $\bar{4}$2m, 4/mmm, 23, m$\bar{3}$, 432, and $\bar{4}$3m.
Zero, two, or four of 2+, 21+, 2−, and 21− simultaneously appear on the slab center plane. If 2+, 21+, 2−, or 21− appear on a slab centre plane with $\bar{1}$, the two-fold axes typically penetrate $\bar{1}$. The exceptions are space group type Fm$\bar{3}$c (number 226) where $\bar{1}$ coexists with either “2+ and 2−” or “21+ and 21−” that do not penetrate $\bar{1}$ in the slab centre plane or space group type Fd$\bar{3}$c (228) where $\bar{1}$ coexists with either “2+ and 21−” or “21+ and 2−” that do not penetrate $\bar{1}$, respectively. We attempt to use |aD| = |bD| and aD · bD = 0 if possible, but we need to compromise with 2|aD| = |bD| and aD · bD = 0 in a number of cases. The flow to determine the principal isometry is as follows:
[Flow 1: check whether T1 is a 45° or 135° rotation]
[Flow 2: primary isometry when T2 is a 45° or 135° rotation]
[Flow 3: T2 and primary isometry in other cases]
Mirror and glide can be an in-plane isometry in point groups 4/m, 4/mmm, m$\bar{3}$, and m$\bar{3}$m, but observation of symmetry operations of all relevant space groups show that inversion centres and/or two-fold screw axes are included in the mirror/glide plane (the former is absent only in space group type Pa$\bar{3}$, number 205) and therefore mirror and glide cannot be a principal isometry.
3.2.3 Potential centred rectangular (oc)The tentative basis vectors are taken such that aT > bT, γT = 90°, and det(T1) = 2. There are four possible T1 based on whether aP = bP or not, which arises from whether aT is smaller than $\sqrt{3} b_{\text{T}}$ or not, and whether γP is larger than 90° or not. The logic in section 3.2.4 is used, however, I and T90deg are used instead of T21 and T12, respectively, as T2.
3.2.4 OtherThe 2 × 1 supercell (T2 = T21) is used as the (hkl) n-double cell, if possible, to make the aspect ratio aD:bD closer to 1:1 under the constraint that aP < bP always hold. However, in some cases, the 1 × 2 supercell (T2 = T12) must be used instead. Two-fold rotation and screw axes do not necessarily penetrate inversion centres and vice versa. Principal isometries and are chosen in the following order.
The next step is to determine which atoms to remove from the (hkl) n-double cell. Atoms with z-coordinates outside ranges z− + m ≤ z ≤ z+ + m are removed. The centre of the largest interval between x-coordinates of atoms at the top layer is taken as xc where 0 ≤ xc < 0.5. However, two distinct choices could be possible if all gaps are 0.25, here the two smallest choices of xc are both considered. This interval is taken regardless of the atom species. For instance, if the topmost layer contains atoms at x = 0.3 + 0.5m and x = 0.4 + 0.5m, then xc = 0.1. In another example, if atoms exist at x = 0.1 + 0.5m and x = 0.35 + 0.5m, then both xc = 0.225 and xc = 0.475 are considered. Atoms with x-coordinates
| \begin{equation} x_{\text{c}} + m < x < x_{\text{c}} + 0.5 + m \end{equation} | (13) |
The final task is to remove half of the atoms in the bottom layer while retaining a principal isometry that keeps the slab nonpolar. Principal isometries m and b have the form
| \begin{equation} (\boldsymbol{{W}},\boldsymbol{{w}}) = \begin{pmatrix} 1 & 0 & w_{13} & w_{1}\\ 0 & 1 & * & *\\ 0 & 0 & \bar{1} & 2z_{c} \end{pmatrix} . \end{equation} | (14) |
| \begin{equation} x_{c} + w_{13}(z_{c} - z_{-}) + m < x < x_{c} + w_{13}(z_{c} - z_{-}) + 0.5 + m \end{equation} | (15) |
| \begin{equation} (\boldsymbol{{W}},\boldsymbol{{w}}) = \begin{pmatrix} \bar{1} & 0 & 0 & w_{1}\\ 0 & * & * & *\\ 0 & 0 & \bar{1} & 2z_{c} \end{pmatrix} , \end{equation} | (16) |
| \begin{equation} - (x_{c} + 0.5) + w_{1} + m < x < - x_{c} + w_{1} + m. \end{equation} | (17) |
Nonpolar and stoichiometric surfaces of normal spinel ZnAl2O4 are obtained as a demonstration of the proposed algorithm. The Wyckoff positions are Zn: 8b sites, Al: 16c sites, and O: 32e sites (u = 0.888 in origin choice 2). Figure 2 shows examples of obtained slabs. In a (111) slab, the top layer is always an Al layer but the second layer can be either an O layer or a Zn layer. T2 = I and the principal isometry is 2y penetrating $\bar{1}$ in all cases. There are two and one choices of xc in the former [Fig. 2(a)–(b)] and latter [Fig. 2(c)] cases, respectively. The (100) slab is a potential tp with T1 = I and T2 = T45deg. The top layer is either the Al–O layer or the Zn layer. One example is shown in Fig. 2(d). The principal isometry for the termination in Fig. 2(d) is 2y and $\bar{1}$ for even and odd layers containing Al, respectively. The (110) slab is not hp, tp, nor oc and T1 = I. The top layer is either the Al–Zn layer or the Al–Zn–O layer and the principal isometry is either 2y penetrating $\bar{1}$ or 21y penetrating $\bar{1}$ [an example of the latter is shown in Fig. 2(e)]. The (210) slab is a potential oc and T1 is determined accordingly. The top layer is either the Al–Zn layer or the Al–Zn–O layer. T2 = I in all cases and the principal isometry is $\bar{1}$ or 2y [an example of the latter is shown in Fig. 2(f)].
Examples of generated ZnAl2O4 surface slabs. Red, blue, and grey balls indicate Zn, Al, and O, respectively. The surface orientation is (a–c) (111), (d) (100), (e) (110), and (f) (210), respectively.
Here we consider making a nonpolar faceted model containing terraces and steps where the orientation of the terrace, which is one of the facets, is given. The flow of the algorithm is as follows:
A simple example of a face-centred cubic (fcc) crystal with the (111) terrace plane and (211) vicinal plane is shown in Fig. 3.24)
Generation of a face-centred cubic step model with a (111) terrace plane and (211) vicinal plane.
In a typical situation, the orientation of the terrace facet in a step model is given, which is denoted as (htOP, ktOP, ltOP). Based on a unit cell with basis vectors (a, b, c) where tOP = htOPa + ktOPb + ltOPc is an out-of-plane vector (terrace out-of-plane vector), e = hea + keb + lec, t = hta + ktb + ltc, v = hva + kvb + lvc, and the vicinal out-of-plane vector vOP = hvOPa + kvOPb + lvOPc must satisfy these relations:
| \begin{equation} (h_{\text{e}},k_{\text{e}},l_{\text{e}}) \cdot (h_{\text{tOP}},k_{\text{tOP}},l_{\text{tOP}}) = 0 \end{equation} | (18) |
| \begin{equation} (h_{\text{t}},k_{\text{t}},l_{\text{t}}) \cdot (h_{\text{tOP}},k_{\text{tOP}},l_{\text{tOP}}) = 0 \end{equation} | (19) |
| \begin{equation} (h_{\text{e}},k_{\text{e}},l_{\text{e}}) \cdot (h_{\text{vOP}},k_{\text{vOP}},l_{\text{vOP}}) = 0 \end{equation} | (20) |
| \begin{equation} (h_{\text{v}},k_{\text{v}},l_{\text{v}}) \cdot (h_{\text{vOP}},k_{\text{vOP}},l_{\text{vOP}}) = 0. \end{equation} | (21) |
An empty cell is generated where a virtual atom is placed at each lattice point. In other words, a virtual atom is at the position w if (I, w) is an isometry of the original cell.9) The empty cell is essentially the same as the actual cell in the fcc example because, by convention, atoms are positioned at lattice points that are also virtual atom sites of the empty cell. Edge vector candidates are simply positions of in-plane virtual atoms when basis vectors are chosen such that tOP is a terrace out-of-plane vector (eq. (18)).
The next step is to determine t that is the shortest among choices that results in the same step model. One vector in the terrace plane that is linearly independent from e, which is
| \begin{equation} \boldsymbol{{t}}' = h_{\text{t}'}\mathbf{a} + k_{\text{t}'}\mathbf{b} + l_{\text{t}'}\mathbf{c}, \end{equation} | (22) |
| \begin{equation} (h_{\text{t}'},k_{\text{t}'},l_{\text{t}'}) = (h_{\text{e}},k_{\text{e}},l_{\text{e}}) \times (h_{\text{tOP}},k_{\text{tOP}},l_{\text{tOP}}). \end{equation} | (23) |
| \begin{equation} \boldsymbol{{t}}'' = f\boldsymbol{{t}}' + g\boldsymbol{{e}} \end{equation} | (24) |
| \begin{equation} \boldsymbol{{t}} = \boldsymbol{{t}}'' + m\boldsymbol{{e}}. \end{equation} | (25) |
We then look for candidates of v in the adjacent plane above the plane spanned by t′ and e. A virtual atom s′ exists in this adjacent plane when (t × e) · s′ is the smallest allowed positive number. There are an infinite number of choices of s′, but any one can be used. The unique multiple of tOP that corresponds to a point on this plane, $\boldsymbol{{t}}'_{\text{OP}} = \alpha \boldsymbol{{t}}_{\text{OP}}$ where α > 0, satisfies $(\boldsymbol{{t}} \times \boldsymbol{{e}}) \cdot \boldsymbol{{s}}' = (\boldsymbol{{t}} \times \boldsymbol{{e}}) \cdot \boldsymbol{{t}}'_{\text{OP}}$. $\boldsymbol{{t}}'_{\text{OP}}$ is not necessarily a lattice vector. Candidates for v are obtained as
| \begin{equation} \boldsymbol{{v}} = \boldsymbol{{t}}'_{\text{OP}} + f'\boldsymbol{{t}}' + g'\boldsymbol{{e}} + m'\boldsymbol{{e}} \end{equation} | (26) |
| \begin{equation} (h_{\text{v${'}$OP}},k_{\text{v${'}$OP}},l_{\text{v${'}$OP}}) = (h_{\text{v}},k_{\text{v}},l_{\text{v}}) \times (h_{\text{e}},k_{\text{e}},l_{\text{e}}) \end{equation} | (27) |
| \begin{equation} (h_{\text{vOP}},k_{\text{vOP}},l_{\text{vOP}}) = \beta (h_{\text{v${'}$OP}},k_{\text{v${'}$OP}},l_{\text{v${'}$OP}}) \end{equation} | (28) |
Vectors projected to a plane normal to e are defined, namely
| \begin{equation} \boldsymbol{{t}}_{\bot} = \boldsymbol{{t}} - \frac{\boldsymbol{{t}} \cdot \boldsymbol{{e}}}{| \boldsymbol{{e}} |^{2}}\boldsymbol{{e}} \end{equation} | (29) |
| \begin{equation} \boldsymbol{{v}}_{\bot} = \boldsymbol{{v}} - \frac{\boldsymbol{{v}} \cdot \boldsymbol{{e}}}{| \boldsymbol{{e}} |^{2}}\boldsymbol{{e}} \end{equation} | (30) |
| \begin{equation} \boldsymbol{{s}}_{\bot} = \boldsymbol{{v}}_{\bot} - \boldsymbol{{t}}_{\bot} \end{equation} | (31) |
The angle between v⊥ and s⊥ is defined as σ, and that between and v⊥ and t⊥ as τ. σ + τ < 90° is imposed as a limitation to prevent an overhanging step. Therefore, t and v should be carefully chosen. σ depends only on v, thus a powerful strategy is to determine v first based on |v⊥| and σ, identify the maximum allowed |t⊥| that is
| \begin{equation} | \boldsymbol{{t}}_{\bot} |_{\text{max}} = | \boldsymbol{{v}}_{\bot} |\cos \sigma, \end{equation} | (32) |
A nonpolar slab without facets is obtained using the algorithm in Hinuma et al.9) However, (v, e, vOP) is used instead of (aP, bP, c1P), the in-plane isometry guaranteeing nonpolarity is limited to $\bar{1}$, 2y, 21y, m, and b, and stoichiometry is not checked at this point. The principal isometry is prioritized in the order $\bar{1}$, 2y or 21y, then m or b, and the isometry with smallest w1 is chosen if there are multiple isometries with the same priority.
The ridge point (r1, r2, r3)T is defined in fractional coordinates as the highest point of a facet, where 0 ≤ r1 < 1 and 0 ≤ r3 < 1 are imposed for uniqueness. It is possible to translate the origin by (r1, 0, r3)T and describe the slab in Cartesian coordinates where the x and y axes are parallel to v⊥ and e, respectively, and the positions of atoms in one supercell are limited to 0 ≤ x < |v⊥| and zbottom ≤ z ≤ ztop. ztop should be not negative and not excessively large. A list of atoms to be removed are made that includes
Removal of atoms in the wedges give the step model. However no atoms exist in wedges in the fcc example and thus there is no atom removal (Fig. 3).
4.3 Additional examplesAnatase TiO2 (space group type I41/amd, number 141, calculated lattice parameters a = b = 3.836 Å, c = 9.406 Å) with the (101) terrace plane is discussed here as a more complicated example. Figure 4(a) gives the crystal structure, where blue and red dots indicate Ti and O, respectively. The origin is taken at an inversion centre such that virtual atom sites lie on inversion centres (shown in green dots). As a consequence, Ti and O sites are not in the conventional positions. Figure 4(b) shows the (101) terrace plane (blue), the (314) vicinal plane (red), and various relevant vectors superimposed on Fig. 4(a). The choice of out-of-plane terrace and edge vectors are tOP = a + c (light blue arrow) and e = a/2 + b/2 − c/2 (green arrow), respectively. Figure 4(c) shows the crystal structure and vectors in Fig. 4(b) from a different direction such that the vicinal (314) plane appears like a line. This figure stresses that t and e lie on the terrace plane while v and e lie on the vicinal plane. Figure 4(d) describes the procedure to determine t. Equation (23) gives t′ = a/2 − b − c/2, and the green band shows the region of ft′ + ge (eq. (24)). One edge of the band is parallel to e and its length is |e|. The band extends along t′ while retaining its width, and protrudes toward the direction out of the sheet in Fig. 4(b); a lighter tone indicates points closer to the reader. t′ is shown with a gradated yellow arrow, and t′′ (gradated brown arrow) must be a lattice vector from the origin to some virtual atom (green dot) in the band. This example uses t′′ = 2t′ and therefore the arrows for t′′ and t′ coincidentally overlap; however, t′′ and t′ does not need to be parallel to each other. t = 2t′ − e = a/2 − 5b/2 − c/2 (gradated blue arrow) is chosen to have the minimum allowed norm with m = −1 in eq. (25) (this −e is shown as the green dashed arrow). The vicinal vector is chosen as v = tOP/2 + 7t′/3 − 4e/3 = a − 3b (gradated red arrow). tOP/2 + 7t′/3 − e/3 = 3a/2 − 5b/2 − c/2 is also a valid choice, but the former is marginally shorter. Both t and v point out of the sheet from the origin, therefore gradated arrows are used. The step vector is therefore s = tOP/2 + t′/3 − e/3 = a/2 − b/2 + c/2 (purple arrow). The vicinal out-of-plane vector is vOP = 3a + b + 4c (eq. (27)), indicating that the vicinal plane in the step model has the (314) orientation. Figure 4(e) is viewed such that e points vertically into the sheet. Both terrace and vicinal planes appear as lines in the figure. Angles σ and τ are inner angles of a triangle with edge lengths |v⊥| = 11.82, |t⊥| = 10.77, and |s⊥| = 3.59 Å, hence σ = 17.5° and τ = 64.3°, respectively. The ridge point in Fig. 4(f) is an O atom, and ztop = 0 is used. One wedge is defined by atoms that satisfy z > (−x + ε) tan σ and z > (−|v⊥| + x + ε) tan τ, where ε is a sufficiently small number, for example 10−4 Å. Atoms in the wedge and those related by translation and inversion symmetry are removed to obtain the step model in Fig. 4(f). However, further removal of atoms in the orange rectangles are necessary to obtain the step model reported in Ref. 20).
Generation of an anatase TiO2 step model.
As another example, step models for the CeO2 (111) surface in the literature can be obtained using the outlined algorithm. An example of the choice of vectors is e = −a/2 + c/2, t = −a/2 + 3b/2 − c, and v = −a/2 + 2b − c/2 for type I, e = a/2 − c/2, t = a/2 − 3b/2 + c, and v = a − 3b/2 + 3c/2 for type II, and e = a − b/2 − c/2, t = −b + c, and v = −b/2 + 3c/2 for type III in Ref. 19), respectively.
Algorithms to generate nonpolar slab-and-vacuum models after removal of surface atoms were developed. Outermost layer modification in a stripe pattern made possible generation of many nonpolar type C slabs including all spinel and perovskite aristotype surfaces. In addition, a procedure that assists generation of step models was proposed that facilitates investigation of step edges that are often active reaction sites.
This study was funded by the “Materials research by Information Integration” Initiative (MI2I) project of the Support Program for Starting Up Innovation Hub as well as a grant (No. JPMJCR17J3) from CREST of the Japan Science and Technology Agency (JST), and by a Kakenhi Grant-in-Aid (No. 18K04692) from the Japan Society for the Promotion of Science (JSPS). The VESTA code37) was used to draw Fig. 3 and 4.
