Modeling interfaces of fluorite-structure compounds using slab charge distribution

ABSTRACT Automated generation of reasonable atomic-level interface models, for example, at a grain boundary, is often a difficult task. The interface modeling algorithm for elementary substances based on charge densities of slab surfaces by Hinuma et al. [AIP Advances 11, 115020 (2021)] was applied to obtain Σ3 (111)/(11 ) and Σ5 (310)/(3 ) interface models of fluorite structure compounds reported in the ICSD database. The algorithm found only one type of in-plane rigid-body translation (RBT) in the former. In contrast, there were diverse interfaces with various RBTs in the latter; the RBT for each compound was identified by also testing a set of RBTs, given by the algorithm, from other compounds. The algorithm by Hinuma et al. can therefore be used, although with caveats, as a complementary tool to estimate the atom configuration at interfaces of compounds. Graphical abstract


Introduction
Information on the atomic configuration at an interface in materials or devices is important when investigating properties through computational methods to attain further improvement. Grain boundaries, in particular, have been studied extensively because these significantly affect various properties of materials, especially concerning diffusion [1][2][3], corrosion (intergranular corrosion) [4][5][6], fracture [7,8], and strength and deformation via affecting dislocation behaviors [9][10][11][12][13]. There are two large topics of interest, where one is understanding what properties, functions or phenomena grain boundaries have or induce in each material, and the other is investigating what types of grain boundaries are favorable, from geometric or crystallographic viewpoints, in various materials. Understanding the atomic configuration is crucial in both cases.
Atom-level resolution imaging has become possible through advanced microscopy methods, such as highangle annular dark-field scanning transmission electron microscopy (HAADF-STEM) [14]. Although information on atom positions at a specific grain boundary (interface) is helpful when understanding the atomic configuration interface, not all atoms are always visible in microscopy images, especially light atoms such as row 1 and 2 elements in the periodic table. Moreover, TEM images are in two dimensions (2D) even though atoms at an interface are arranged in 3D. Computational approaches complement experiments when obtaining the 3D configuration of all atoms.
The two sides of an interface are typically modeled computationally as slabs infinitely extending in 2D separated by vacuum in the remaining direction. Determination of the three degrees of translational freedom, ξ, η, and ζ (usually denoted as the rigidbody translation, or RBT) when forming an interface between two crystals is an important issue. The RBT must be specified to uniquely identify how to position two sides of an interface with respect to each other. Slabs 1 and 2 represents the two sides of an interface in Figures 1(b) and 2(b). Applying an arbitrary RBT of (ξ, η, ζ) = (0.1, 0.3, −0.235) to slab 2 results in Figures 1(c) and 2(c), respectively. Minor relaxations of atoms at the interface may be resolved through first-principles calculations, but the values of ξ, η, and ζ must be identified within a reasonable range to guide firstprinciples calculations to relax to a desirable interface. Although traditional grain boundary model generation typically requires a 3D commensurate site lattice (CSL), there is a procedure to obtain a 3D grain boundary model as long as there is a 2D commensurate lattice at the interface [15], but still the RBT must be obtained separately in Ref [15].
Originally introduced by Vítek [16], the generalized stacking fault energy, or γ-surface, is obtained as a 2D map of the minimum interface energy (typically denoted as γ) in ξη-space. Here, ζ is chosen to minimize γ for each combination of ξ and η to obtain a 2D surface. The γ-surface is typically derived computationally because experimental evaluation of the γsurface forces synthesis of unstable atomic configurations, but much computational resources are generally required if calculations are to be carried out for all combinations of ξ and η on a dense 2D grid. Gholizadeh et al. calculated a limited number of points to fit the coefficients of the Fourier expansion of the γsurface [17]. Kiyohara et al. used a Gaussian processbased kriging technique to obtain CSL grain boundaries of Fe with two orders of magnitude smaller number of first-principles calculations than using a brute-force grid search [18]. Gao et al. devised an interfacially confined swarm intelligence algorithm that is efficient in exploration of a potential energy surface [19]. The ARTEMIS code by Taylor et al. handled interfacial alignment, or translational shift of one side of the interface against another, by optimizing a bond ideality factor that considers the change in bonding environment [20]. Hinuma et al. used the overlap of charge density distribution at slab surfaces to obtain the RBT of an interface and demonstrated its effectiveness on interfaces of bcc, fcc, and hcp elementary metals [21]. The last method uses first principles calculations only when obtaining the slab charge densities, which consumes very little computational resources, and for final relaxation of interfaces once reasonable RBT candidates are found and applied.
This paper applied the procedure in Hinuma et al [21]. to binary compounds, namely fluorite structure compounds with stoichiometry MX 2 , because of its high symmetry (space group Fm � 3m, number 225) and diversity of constituent element combinations. The M-sites are eight-fold coordinated and form a bcc lattice, while the X-sites are four-fold coordinated, interstitial sites. The Σ3 (111)/(11 � 1) and Σ5 (310)/(3 � 10) interfaces, corresponding to tilt grain boundaries, were considered. Atom positions of (111)/(11 � 1) interfaces were quite limited. In contrast, a wide variety of atom positions were found in (310)/(3 � 10) interfaces, illustrating the necessity of computation as a tool to complement experiments in identifying atom positions. However, information on suggested RBTs from various compounds were often necessary to identify the (310)/(3 � 10) interface in many cases, thus the algorithm by Hinuma et al [21]. needs to be used with caution.

Methodology
The RBT determination algorithm by Hinuma et al [21]. was applied to fluorite structure compounds reported in the ICSD database [22,23] according to data downloaded in June 2021. Compounds containing lanthanides (Ce-Lu) or atoms heavier than Po inclusive were removed from consideration, resulting in a set of 91 compounds. Evaluation of compounds containing Ce is very attractive especially because some experimental results on CeO 2 are reported [24,25]. However, most of the Ce compounds encountered serious convergence problems. One possible reason is the difficulty in finding the stable state of electrons in under coordinated Ce and/or Ce with valence other than 4+. For instance, an electron added to Ce 4+ may enter either 4f or 5d states, but the choice is not always obvious. Adding Hubbard U to d or f states may control the outcome and could attain convergence, but choosing the appropriate U such that the resulting state is sufficiently reflecting reality is a non-trivial task with no straightforward answer. As a result, Ce containing compounds were excluded in this study.
First-principles calculations were conducted using the projector augmented-wave method [26] and the Perdew-Burke-Ernzerhof functional revised for solids (PBEsol) [27] as implemented in the VASP code [28,29]. There were three types of first-principles calculations: bulk calculations to determine the lattice parameters, slab calculations to obtain the charge density at the surface, and interface calculations to relax atom coordinates from an initial structure with a given RBT. The cutoff energy for the planewave basis (ENCUT tag) was explicitly fixed to 400 eV except for bulk relaxation calculations where 550 eV was used. The lattice parameters were fixed in slab and interface calculations, but not in bulk calculations. The k-point mesh in slab and interface calculations was sampled at 0.15 Å −1 in the in-plane directions (directions parallel to the surface) and one point in the remaining out-of-plane direction. The VASP pseudopotential (POTCAR) files and the number of valence electrons per atom used in this study are given in Table S1 in Supplemental material.
For charge density calculations using slab models, atom coordinates were fixed to those of cleaved bulk. The average atom size, defined as the cube root of the bulk volume per atom, ranged from 1.726 Å in CoH 2 to 3.685 Å in TeRb 2 . The supercells for {111} slabs were orthorhombic and contained three and nine repeat units of atoms and vacuum, respectively (Figure 1(a)). The supercell size ranged from a = 2.793, b = 4.838, and c = 23.370 Å in CoH 2 to a = 5.966, b = 10.334, and c = 58.456 Å in TeRb 2 . On the other hand, the supercells for {310} slabs were also orthorhombic and contained 10 and 30 repeat units of atoms and vacuum, respectively (Figure 2(a)). The supercell size ranged from a = 3.951, b = 6.246, and c = 24.985 Å in CoH 2 to a = 8.437, b = 13.341, and c = 53.363 Å in TeRb 2 . The charge density mesh was fixed to 60 × 108 × 576 and 84 × 140 × 540 for {111} and {310} models, respectively. The same charge density mesh was used regardless of the supercell size to ease comparison between compounds.
For interface relaxation calculations, the charge density mesh was automatically determined by the VASP code using the precision (PREC) tag set to Normal. The supercell size was the same as those in charge density calculations.
The algorithm to determine the RBT between two slabs [21] finds local maxima of the function The two slabs, 1 and 2, are centered at z = 0 and 0.5, respectively, and there is a vacuum region between the two slabs (Figures 1(b) and 2(b)). The RBT is defined as �; η; ζ ð Þ, which is the translation of slab 2 with respect to slab 1. Decreasing ζ from 0 has the effect of bringing the two slabs closer to each other (Figures 1(c) and 2(c)). The charge densities of the two slabs are obtained on a grid of X×Y×Z points. The filtered charge density, ρ 0 ¼ ρ exp À ρ=λ ð Þ, is obtained from the charge density ρ. The absolute value of the charge density is used for clarity; the function u is obtained as a product of two charge densities, thus the sign of charge densities is irrelevant. The value of ρ 0 is maximized at ρ ¼ λ and becomes smaller in vacuum (ρ ! 0) and in the slab (very large ρ everywhere). The value of u becomes larger where the filtered charge densities of slabs 1 and 2 overlap. From another viewpoint, the reaction front is defined as the region where the charge density is close to the parameter λ, and the function u is maximized when the reaction fronts of slabs 1 and 2 overlap. The RBT strongly depends on the value of λ, thus multiple choices need to be investigated; λ = 0.01, 0.016, and 0.025 q/Å 3 were considered in this study, where q is the elementary charge. The unit of the charge density in Ref [21]. is erroneously denoted as eV/Å 3 instead of q/Å 3 . A too large λ results in a large ρ 0 in bulk. Therefore, the skin of the electron charge cloud of each surface slab that acts as the reaction front for interface formation cannot be sampled. On the other hand, an excessively small λ results in reaction fronts that are too far apart, thus the charge density near the surface cannot be evaluated well. The value of ζ often results in too large separation of the two slabs, which could be remedied by relaxation using first-principles calculations, for example, to bring the slabs together. On the other hand, finding reasonable � and η (inplane RBT) is important. The value of u was explicitly obtained using Equation (1)  (310)/(3 � 10) interfaces (Σ3 and Σ5, respectively) predicted by the RBT algorithm were discussed for systems where the in-plane RBTs were correctly obtained for the (111)/(111) and (310)/(310) interfaces, respectively.

(111)/(111) interface
The correct RBTs of the (111)/(111) interface are (0,0,-0.25) and (0.5,0.5,-0.25). This is simply cleaving of a (111) slab, and there are two solutions because (0.5,0.5,0) is a lattice vector. Figure 3 shows whether the correct solution of the in-plane RBTs were obtained with the three considered λ choices in each compound. Symbols 'n', 'O', 'X', and 'M' means that no solution was obtained with the RBT Figure 3. In-plane RBT solutions found for the (111)/(111) interface with constituent elements M and X forming MX 2 . Symbols "n", "O", "X", and "M" means that no solution was obtained with the RBT algorithm over all considered λ, at least one of the correct solutions was included in all λ that gave solutions, the correct solutions were never found in λ that gave solutions, and one of the correct solutions was included in at least one λ but not in at least one λ, respectively. The red character for OAl 2 indicates that the fluorite structure was not retained in the (111)/(11 � 1) interface. algorithm over all λ (a result of choosing too large λ), at least one of the correct solutions was included in all λ that gave solutions, the correct solutions were never found in λ that gave solutions, and one of the correct solutions was included in at least one λ but not in at least one λ ('M'ixed), respectively.
A number of trends are found in Figure 3. In hydrides, H always occupies the X site, the correct solution was obtained when M is less electronegative (groups 1-5 and Cr), with the exception of M=Ba, while no solution was found for other M. The anti-fluorite structure forms with a group 16 X and a group 1 M. The correct solution was found when M is Li or Na, and not found (except for ORb 2 ) when M is K or Rb. This difference could come from the number of valence electrons in the group 1 element. There is one valence electron for Li and Na, while K and Rb has seven and nine, respectively. The charge density landscape is completely different between these two; Li and Na nominally has zero valence electrons at formal charge +1, but K and Rb has six and eight valence electrons, respectively. The antifluorite structure also forms with group 14 M and group 2 X, where no solution was obtained. There is a small cluster of fluorite structure compounds with group 15 M and group 9 X (PRh 2 , PIr 2 , and AsRh 2 ) with the correct in-plane RBT. Group 13 X compounds could not find the correct solution except AlO 2 and AuIn 2 , while the correct solution was obtained with late transition metal M (groups 8-10) and group 14 and 15 X (Si, Sn and N). The correct solution was always found for group 16 X, and for group 17 X combined with less electronegative M.

(111)/(11 � 1) interface
(111)/(11 � 1) interface calculations were conducted, for all RBTs suggested by the RBT algorithm using three λ choices, for compounds labeled with symbol 'O' or 'M' in Figure 3 because the charge density of the {111} surface is reliable in the RBT algorithm. Relaxation using first-principles calculations showed that the most stable in-plane RBT is one of (0,0), (0.5, 0.5), (0.666, 0), and (0.166, 0.5) in all system. The relaxed structures for CaF 2 with RBTs (0,0) and (0.666, 0) are shown in Figure 4(a,b), respectively. These four inplane RBTs result in the same grain boundary, suggesting that this is the only reasonable solution for the (111)/(11 � 1) interface. This CaF 2 interface was what was found in most compounds using calculations and experimentally in yttrium-stabilized zirconia (YSZ) [30] and CeO 2 [24,25]. The exceptions are OAl 2 (fluorite structure was not retained), SiO 2 (gap between two slabs did not close enough), and FeSi 2 , RhN 2 , and NiS 2 (these interfaces are shown in Figure S1).

(310)/(310) interface
The correct RBT of the (310)/(310) interface is (0,0,-0.25). Figure 5 shows whether the correct solution was obtained using the RBT algorithm.  the correct solution, as with the (111)/(111) interface. Group 13 X compounds successfully derived the solution, while with the correct in-plane RBT. Group 13 and 14 X compounds found the correct solution while X = N and O did not, which is the opposite of the (111)/(111) interface. There were more 'no solution' results in the (310)/(310) interface than the (111)/(111) interface. The atoms on the (310) surface are more sparsely distributed than the (111) surface. When charge densities are concentrated near a certain element, a smaller charge density threshold is necessary to form a smooth reaction front covering the entire surface. This translates to using a smaller value of λ, which could be too small to obtain reasonable RBTs that bring the two surfaces sufficiently close together.

(310)/(3 � 10) interface
Similar to as in the (111)/(11 � 1) interface, (310)/(3 � 10) interface calculations were conducted, for all RBTs suggested by the RBT algorithm using the three λ choices, in compounds labeled with symbol 'O' or 'M' in Figure 5. Relaxation during first-principles calculations of the (310)/(3 � 10) interface resulted in significant relaxation of atom coordinates, to the extent that the fluorite structure was not retained, in OAl 2 , IrSn 2 , RuSi 2 , NiS 2 , ORb 2 , SH 2 , SeH 2 , FeSi 2 , MoB 2 , CaH 2 , and WSi 2 , and no RBT was suggested in SrF 2 and BaH 2 , thus these are not discussed hereon. These compounds are denoted using red characters in Figure 5. There was only one symmetrically distinct RBT within 0.005 eV/unit cell from the lowest energy interface, and the corresponding RBTs are Figure 5. In-Plane RBT solutions found for the (310)/(310) interface with constituent elements M and X forming MX 2 . Symbols "n", "O", "X", and "M" means that no solution was obtained with the RBT algorithm over all considered λ, the correct solution was included in all λ that gave solutions, the correct solution was never found in λ that gave solutions, and one of the correct solutions was included in at least one λ but not in at least one λ, respectively. Red characters indicate that the fluorite structure was not retained in the (310)/(3 � 10) interface.
shown in Table 1 In-plane RBTs of (x,y), (-x,y), (x,0.4-y), and (x,0.4-y) are symmetrically equivalent based on the internal coordinates of slab models. The in-plane RBTs in Table 1  2) (inplane RBT type, labeled as A-I, respectively). The inplane RBT type derived using the charge density algorithm is labeled as 'Algorithm' in Table 1.
Interface models were obtained with one in-plane RBT from A-I and ζ=-0.235, and relaxation calculations were performed for the nine initial interface models. This procedure is based on the concept of transfer learning in machine learning, where knowledge acquired from a particular situation is used to refine understanding of a related situation. The inplane RBT(s) that resulted in 0.005 eV/unit cell from the lowest interface is given in the column 'Transfer' in Table 1. Multiple initial RBTs resulted in the same interface in some situations, which means that some systems can accommodate more relaxation compared to others, and multiple characters are given in the column 'Transfer' in Table 1. The RBT in 'Algorithm' was found in 'Transfer' in only one-third of the compounds in Table 1. Therefore, the RBT predicting power of the algorithm may appear relatively bad. However, the algorithm has provided some interesting insights. Figure 6 shows the most stable interface model for selected compounds, while Figures S2-4 shows results for all compounds. The characters A-I in these figures correspond to the nine RBTs that resulted in the stable interface model; multiple initial RBTs resulted in the same model in some cases. A kite-shaped pattern as viewed along the [001] direction, with M (blue) atoms at the vertices and ratio of edge lengths approximately 1:3:3:1 (shown in orange in the left panel in each of Figure 6(a-c)), was frequently encountered, but was not found in some systems (example in Figure 6(d)). This kite motif is experimentally observed in CeO 2 [31] but not in YSZ [14,32]. Appearance of the same kite motif of M atoms does not mean that the positions of the X atoms are similar, nor that there is only one type of pattern when viewed perpendicular to the [001] direction (right panel in each of Figure 6(a-d)). Two small kite patterns, both with edge lengths roughly 2:1:1:2, where one has M atoms and the other has X atoms at the three vertices with small angles, alternate along the interface in the direction perpendicular to [001] in PtSn 2 and PIr 2 ( Figure 6(a,  b), respectively). However, the planes of M and X atoms perpendicular to [001] at the same position on both sides of the interface in PtSn 2 (all atoms on the yellow dashed lines are M) but M atoms on one side of the interface and X atoms on the other share the same plane in PIr 2 (atoms on one side on the interface is M and the other side is O, as indicated by the yellow dashed line in Figure 6(b)). In CaF 2 , the large kites are found but the small kites are not, and X atoms are arranged in a complicated pattern (Figure 6(c)). In PtAl 2 , M and X sites overlap at the interface when viewed along the [001] direction, which is the result of a very narrow interface region and out-of-phase M and X layers perpendicular to [001] (Figure 6(d)).
The diverse interfaces in Figure 6 and Figures S2-S4, which were obtained from the same set of RBTs, are very difficult to derive with intuition only. The X atoms are typically lighter than M atoms that can form complicated patterns, (exceptions in Table 1 are PRh 2 , AsRh 2 , and PIr 2 ), and computational efforts are crucial to the understanding of atom positions the interface. The algorithm in Hinuma et al [21] was found to be a powerful tool, although simulations of interfaces with various compounds of the same crystal structure to find a diverse collection of initial RBTs is strongly recommended for more precise estimation.

Summary
This study applied the interface modeling algorithm based on charge densities of slab surfaces in Hinuma et al [21], which was originally applied to elementary substances with fcc, bcc, and hcp crystal structures, to known compounds with the fluorite structure. One type of RBT was consistently found for (111)/ (11 � 1) interface, and the 3D model was often that in Figure 4. In contrast, diverse interfaces shown in Figure 6 and Figures S2-S4 were suggested based on the algorithm. The predicting power of the algorithm for a single compound is relatively weak, but can be improved by using of a collection of RBTs suggested from diverse collection of compounds with the same structure. Further verification of the algorithm in Ref [21]. would improve the breadth of tools for uncovering atom ordering in complex interfaces.