Geometrical prediction of cleavage planes in crystal structures

A new geometrical approach is proposed for finding the most prominent planar gaps/cleavage planes in crystal structures. The approach is realized by using the GALOCS program package and is demonstrated on several inorganic structures.


Introduction
Cleavage is the ability of single crystals to split easily along specifically oriented planes (Hurlbut & Klein, 1977). Such planes are usually parallel to the reticular lattice planes with low Miller indices/large interplanar distances. The cleavage phenomenon was known long before the discovery of X-ray diffraction by crystals; the observation of cleavage in crystals contributed greatly to the ideas about periodicity/long-range order of their structures (Tutton, 1922;Authier, 2013). Cleavage is the most striking example of anisotropy of physical properties (Nye, 1985). Finally, cleavage is the subject of many materials science oriented research (Gilman, 1960;Lawn, 1974): acquiring the information concerning cleavage in a given crystalline material is essential for the understanding of failure, brittle fracture, toughness, plasticity and strength (Lawn et al., 1993).
The prediction and discovery of cleavage planes is important for microelectronics and electro-optics where durability of crystalline materials and prevention of their catastrophic mechanical failures is critical (Spearing, 2000). Whether performed theoretically or experimentally, such prediction presents a challenging task.
The simpler approach to the prediction of cleavage planes involves counting the number of broken chemical bonds per unit area of a candidate plane. Such counting can be carried out analytically for the simplest crystal structures (e.g. of rock salt, diamond or sphalerite type) (Ramaseshan, 1946).
However, it becomes impractical for larger (e.g. ternary) structures. Accurate calculations of cleavage energies are demanding (Zhang et al., 2007;Bitzek et al., 2015) and should involve e.g. density functional theory (DFT) calculation of the surface energies (Ong et al., 2013;Tran et al., 2016), which is half of the cleavage energy [for the case of slowly propagating cracks (Griffith, 1921)].
It is also possible to study cleavage experimentally (Cramer et al., 2000;Field, 1971;Jaccodine, 1963;Lawn et al., 1993;Michot, 1987;Sherman et al., 2008). While some cleavage planes may immediately appear because of a 'mechanical impact', such impact 'experiments' can hardly be helpful in the precise measurement of a cleavage energy. More complex arrangements are implemented for this purpose (Gilman, 1960;Hirsh et al., 2020). Specifically, a single crystal should be cut to a plate whose normal is parallel to a candidate cleavage plane while a uniaxial stress (so-called Mode I stress) must be applied in a controlled manner. In conclusion, existing experimental and theoretical methods are demanding and unsuitable for throughput studies of cleavage without preliminary suggestions of several likely possibilities.
Here, we suggest a simple geometrical algorithm and a computer program GALOCS (gaps' locations in crystal structures) which searches for promising cases of cleavage in single crystals of known structures. Apart from accepting the working hypothesis about the planar character of cleavage, we assume that such cleavage occurs along the planes, exposing the widest planar gaps (the intervals of empty space, appearing when observing the plane lengthwise). The manual variant of this approach is common in crystallography classes, where a lecturer exposes a three-dimensional structural model and demonstrates prominent clearances/gaps in the structure to the audience. While such inspection can be performed both physically and graphically [e.g. by using the VESTA program (Momma & Izumi, 2011)], there are no computational methods for doing it automatically. This methodological gap is filled by the GALOCS package. In addition to its great illustrative potential, GALOCS can serve as a rough predictor of cleavage planes in crystals. If necessary, such planes can be analysed rigorously using DFT. Alternatively, the results may suggest a particular experimental geometry that allows measuring the cleavage energy of the plane.

The algorithm behind GALOCS
The proposed approach realizes an automatic (as opposed to visual) inspection of the 'space-filling function' (SFF). Specifically, we suggest inspecting the sections of the SFF by the planes that are parallel to the lattice planes with reasonable Miller indices/interplanar distances. The dependence of the average SFF on the depth of the plane in the unit cell is produced in order to locate the most prominent planar gaps and measure their width and depth. The correlation between the width/depth of the gap and the cleavage ability of the candidate plane is the major hypothesis behind the algorithm.

Space-filling function
Calculation of the SFF requires knowledge of space group type, lattice parameters, coordinates, occupancies of atomic/ Wyckoff positions and types of atoms. The first definition of the SFF is the superposition of electron densities of individual pseudo-atoms: Here, R m = u mi a i (Einstein summation over the repeated index running from one to three is implemented everywhere throughout the article) are the positions of all the lattice nodes, u mi are arbitrary integers and a i are the lattice basis vectors i = 1, 2, 3. The vectors R ¼ x i a i (0 x i < 1) list the positions of all the atoms in a unit cell. ðrÞ are electron densities of a (pseudo)atom number , which can be approximated using e.g. electron density of an isolated atom (Su & Coppens, 1998) or atomic invarioms (Dittrich et al., 2004). All the atomic electron densities are normalized: R ðrÞdr ¼ Z [Z is the number of electrons, associated with the (pseudo)atom number ].
Another possible definition of the SFF involves an atomic probability density function (PDF). It describes the probability of finding an atom displaced by a vector u from its average position R . Such displacements may originate from either a thermal motion or a static disorder. A PDF corresponds to the average of all the atomic positions over time or over different unit cells. A three-dimensional Gaussian function is the simplest approximation of an atomic PDF (Coppens, 1997): Here, is the second-rank tensor of atomic displacement. The components of this tensor are known from an X-ray or neutron diffraction experiment and are represented as isotropic or anisotropic atomic displacement parameters (U ij or ij ) (Trueblood et al., 1996;Coppens, 1997;Tsirelson & Ozerov, 1996). If is represented by a 3 Â 3 matrix then u T and u are the rows or the column vectors, respectively. The SFF (r) can now be calculated as Here, W describe the 'weights' of every atom in the cumulative SFF. W ¼ 1 yields an SFF in the form of atomic density (it disregards the sort of atoms involved). Setting W to the atomic masses defines the SFF as the mass density. Finally, setting W to Z will produce an SFF in the form of nuclei charge density. Notably, all possible definitions of the density function sustain the periodicity so that (r) = (r + R m ). Therefore, it is sufficient to calculate the values of (r) inside the unit cell only. Because both p ðrÞ and ðrÞ are negligible when |r| > r max (r max ' 1Å ), only the limited number of atoms whose research papers 794 Uriel Vaknin et al. Geometrical prediction of cleavage planes in crystal structures centres are inside the unit cell or within r max from its borders are included in the sums (1) or (3).

Forming a section of a unit cell by an arbitrary candidate cleavage plane
GALOCS inspects the sections of the SFF by the planes that are perpendicular to a unitary vector n (|n| = 1). Let us define the depth D of the plane inside the crystal and introduce the orientation and depth-dependent average À(n, D) of the SFF (r) as Here the averaging of (r) is performed over all the positions r, satisfying the condition rn = D (the equation of a plane, normal to the vector n and standing at the distance D from the origin). The periodicity of the crystal structure means that D is non-negative and below the corresponding interplanar distance d(n) = d hkl . The latter is the inverse length of the primitive reciprocal lattice vector, It is sufficient to average the SFF (r) rn = D over translationally independent locations only, rather than over the entire plane. It is therefore worth transforming the coordinates of r to the coordinate system A 1 , A 2 , A 3 such that A i and a i are the bases of the same crystal lattice, but the vectors A 1 and A 2 are parallel to the plane of interest (hkl) and A 3 connects two adjacent lattice planes (A 3i h i = 1). The number-theoretical algorithms for such a transformation are described elsewhere (Gorfman, 2020) and the corresponding MATLAB-based program MULDIN is deposited there. This step presents the core of the algorithm because it uses the periodicity of the crystal structure. Using the coordinates X 1 , X 2 and X 3 such that r = X i A i = x i a i and expressing the SFF (x 1 , x 2 , x 3 ) as hkl (X 1 , X 2 , X 3 ) yields The cleavage planes are likely to appear where À(hkl, X 3 ) drops to the lowest possible values. Averaging of À(hkl, X 3 ) over X 3 would yield the average of the SFF hÀðhkl; X 3 Þi X 3 ¼ 0 over the entire unit cell. Let us introduce the constant threshold S = y 0 (typically y = 0.75) such that all the X 3 values where À(hkl, X 3 ) < S are considered as planar gaps. We define the maximum effective width of a planar gap as Here ÁX 3 is the length of the longest continuous range of X 3 values where À(hkl, X 3 ) < S and h S À Àðhkl; X 3 Þi ÁX 3 is the average gap depth in this range. The denominator 0 is introduced in order to reduce the dimension of C hkl to Å . Fig.  1 explains definition (7) graphically. It shows an arbitrary À(hkl, D) with the most prominent gap between two vertical dashed lines. Definition (7) suggests that C hkl must be centrosymmetric. Indeed, the sets of Miller indices (hkl) and ð " h h " k k " l lÞ define the same lattice planes. The only difference between them concerns the direction of the plane movement with the increasing value of X 3 . Specifically, Àð " h h " k k " l l; X 3 Þ ¼ Àðhkl; ÀX 3 Þ, meaning that Àð " h h " k k " l l; X 3 Þ has the same minima, the same average value and the same width of the structural gap as À(hkl, X 3 ). Accordingly, C hkl ¼ C" h h " k k " l l , so the symmetry of the C hkl is defined by one of the eleven Laue classes. A graphical illustration of the definition of C hkl according to equation (7). The solid line shows an arbitrary À(hkl, D) function extending over two unit cells (one unit cell on the positive and one unit cell on the negative side of the D axis). The dashed horizontal line shows the threshold value of S = 0.75 0 . The regions of À(hkl, D) < S are filled by colour. The numerator of equation (7) is equal to the area inside the coloured region.

The flow chart of the GALOCS algorithm
GALOCS realizes the following steps: (i) Calculating electron density (x 1 , x 2 , x 3 ) (or any other SFF in future releases) inside a unit cell (0 x i < 1) as a superposition of spherical atoms.
(ii) Generating the set of symmetry-independent (with respect to the relevant Laue class) lattice planes (hkl) such that d hkl > d min .
(iv) Calculating C hkl according to equation (7), and drawing the directional dependence of C hkl on the stereographic projections and polar plots.
The MATLAB-based software package includes several modules for: (a) Automatic reading of relevant structural information from CIFs.
(b) Calculating (x 1 , x 2 , x 3 ) in the unit cell on the predefined grid. The current version of the program uses pre-tabulated spherically symmetrical electron densities in isolated atoms (Su & Coppens, 1998).
(d) Calculation of C hkl (according to equation 7 with userdefined clearance threshold, typically y = 0.75) for the set of lattice planes with interplanar distance above some userdefined value d min .
The software package of GALOCS with corresponding manual is available in the supporting information and through the GitHUB platform. The functionality of the package will be extended (e.g. by introduction to the graphical user interface and by implementing various SFF models) subject to the interest of the user community. Fig. 2 illustrates the GALOCS output for a two-dimensional graphene structure. The conventional unit cell is based on the vectors a 1 and a 2 such that a 1 ¼ a 2 ¼ a ¼ 2:461 Å ; ¼ 120 . The SFF was defined by the superposition of electron densities of isolated carbon atoms. Fig. 2(a) shows the transformed unit cells [based on the vectors A 1 and A 2 with A 1 being parallel to (hk)] for the cases of (10), (11) and (12) planes. Fig. 2(b) illustrates the dependence of C hk on the reciprocal-space directions using the polar plot. The spikes extend in the directions of the reciprocal lattice vectors B hk , their lengths are proportional to the corresponding C hk values. It shows that the longest gap in graphene is present for {10} H and {11} H families of planes [the notation fhkg H stands for set of planes that are symmetry equivalent to (hk) with respect to the operations of the crystal class of two-dimensional hexagonal lattice]. {21} H planes have $15 times smaller C hk (see Table 1 for the numerical values).  A two-dimensional illustration of the GALOCS output for the case of graphene. (a) The SFF for the structure of graphene according to equation (1) and electron densities of isolated carbon atoms. The same figure shows the standard unit cell and its three transformations, suitable for the analysis of the corresponding planar gaps in this structure. In all cases, the transformed unit cell has the vector A 1 parallel to the (hk) planes. (b) The dependence of the C hk on the direction n k B hk . The Miller indices corresponding to each of the directions are marked explicitly. (c) The dependences of À(10, X 3 ), À(11, X 3 ) and À(12, X 3 ) on D = d hk X 3 (vertically shifted relatively to each other for clarity).

Examples
Here, we demonstrate the implementation of the algorithm for the case of three inorganic structures (Si, LiNbO 3 and SiO 2 ). directions. The structure has one symmetry-independent atom sitting at the standard origin 1 of the space group Fd3m. Each atom has four nearest neighbours at the corners of a regular tetrahedron. Table 2 summarizes the relevant structural information.

Silicon/diamond structural type
The calculation of C hkl over the lattice planes with the interplanar distance above d min ¼ 0:2Å involved 17 330 primitive reciprocal lattice vectors enclosed in a corresponding reciprocal-space sphere of radius d À1 min ¼ 5Å

À1
. This set included 446 symmetry-independent planes (with respect to the operations of the point-symmetry group m3m). The cutoff d min is justified by the fact that C hkl drops to zero for those planes that have small interplanar distances. Fig. 4 illustrates this statement by showing C hkl as a function of the length of the primitive reciprocal lattice vector (B hkl ¼ d À1 hkl ). Prominent gaps are present among low B hkl (high d hkl ) planes only. The same figure indicates the Miller indices of the planes with the highest cleavage ability. Note that the indices of the planes are given with respect to the conventional non-primitive (face centred) unit cell and therefore some of the indices are not coprime. The indices are coprime if expressed using a primitive unit-cell basis (Nespolo, 2015).
Figs. 5 and 6 illustrate the striking anisotropy of C hkl in silicon using stereographic projections and polar plots. The stereographic projections contain a false-colour map, which sets one-to-one correspondence between the colour and the C hkl values. The projection is viewed along the zone axis A uvw = [uvw] = ua 1 + va 2 + wa 3 {½1 " 1 10/[111] in Figs. 5(a)/6(a), respectively}. Figs. 5(b) and 6(b) show polar plots of C hkl for all the reciprocal lattice directions B hkl in the zone (i.e. such that B hkl Á A uvw = 0). The individual spikes extend in the directions that are normal to the anticipated cleavage planes. Such polar plots provide specific guidelines for the cleavage experiment, in which a crystal is prepared in the form of a wafer whose surface is normal to A uvw . Cleavage is initiated at a small pre-crack [see e.g. (Hirsh et al., 2020)] and at a uniaxial stress [also known as Mode I stress (Lawn et al., 1993)] along the direction n. Table 3 lists numerical values for the cleavage ability of the most prominent planes.
The calculations predict {111} C and {110} C as the most prominent cleavage planes in silicon [the notation fhklg C stands for the set of planes that are symmetry equivalent to ðhklÞ with respect to the operations of the cubic crystal class m " 3 3m]. These planes are well known from real cleavage experiments on single crystals of silicon .
The output is further illustrated in Figs. 7 and 8. Fig. 7 shows the corresponding sections of electron density by (111) planes and the À(n, D) function where the gap is seen clearly. Fig. 8 shows the view of the structure along the predicted cleavage plane. This figure was produced using the VESTA program, where (111) plane was added artificially and the structural model was rotated in a way that the normal to the plane is parallel to the horizontal axis of the screen. The easy ability to produce the images, exposing the structural gaps in the crystal clearly, is one of the goals of the suggested algorithm.

Figure 4
The B hkl dependence of C hkl values in silicon. Here B hkl is the length of the corresponding reciprocal lattice vector (inverse to d hkl ). The Miller indices are given with respect to the non-primitive (face centred) unit cell; accordingly, some Miller indices [e.g. (220)] have a common divider. When converted to the primitive unit cell, these indices are coprime. Table 2 The structure of silicon used for the calculations.

Space group
Fd " 3 3m ( Table 1 Magnitudes of C hk in graphene. (hk) (10) (11) (12) C hk ðÅ Þ 0.54 0.40 0.03 Table 3 Magnitudes of C hkl in silicon.  Table 4 summarizes the relevant structural information (Levien et al., 1980). Right-handed quartz crystallizes in the space group type P3 2 21 (No. 154) and corresponds to the Laue class 3m. The structure has two symmetryindependent atoms (Si and O). Enumeration of the lattice planes with the interplanar distance above d min ¼ 0:2Å results in 49 358 reciprocal lattice vectors, 4393 of them are symmetry independent (with respect to the point-symmetry operations of the Laue class). Figs. 10-12 and Table 5 show the results of the calculations (as in chapter 3.1). Fig. 10   Anisotropy of C hkl values in silicon. (a) The stereographic projection viewed along the ½1 " 1 10 crystallographic direction displays the false-colour map of the C hkl (dependence of C hkl on the direction). Each point on the stereographic projection corresponds to a specific reciprocal-space direction. The small white points indicate the stereographic projections of the normals to all the lattice planes involved in the calculation. The colour scheme expresses the values C hkl =maxðC hkl Þ according to the colour bar on the left. (b) The polar diagram of the cleavage ability for all the planes in the ½1 " 1 10 zone. The outer circle corresponds to the highest cleavage ability (C 111 = 1.025 Å for the case of silicon). The structure of silicon viewed along the (111) plane. The [111]* reciprocal lattice direction is horizontal. The figure exposes the corresponding structural gap (as also appears in Fig. 7).

Figure 9
The

Figure 10
Same as Fig. 4 but for the case of quartz.

Table 5
Magnitudes of C hkl in quartz.  Table 6 The structure of LiNbO 3 used for the calculations (Weigel et al., 2020).

Space group
R3c ( (organized in the same way as Figs. 7 and 8) illustrate the output for the case of (011) planes in quartz.
Enumeration of the lattice planes with the interplanar distance exceeding d min ¼ 0:2Å results in the generation of 46 340 reciprocal lattice vectors, where 4108 of them are symmetry independent (with respect to the point-symmetry group " 3 3m). The calculation results are presented in Figs. 15-19 and Table 7, organized as in the previous sections. The analysis predicts the most prominent cleavage plane with the Miller indices (012). This plane is also known from the measurements of cleavage energies in LiNbO 3 crystals (Hirsh et al., 2020).      Table 8 Magnitudes of C hkl in AlN [wurzite structural type, space group type P6 3 mc, the structure is reported by Xu & Ching (1993)].

Figure 14
The structure of quartz viewed along the (011) plane. The [011]* reciprocal lattice direction is horizontal.

Figure 15
Same as Fig. 4 but for the case of LiNbO 3 .

Table 9
Magnitudes of C hkl in CaF 2 [fluorite structural type, space group type Fm " 3 3m, the structure is reported by Speziale & Duffy, 2002)].  Table 10 Magnitudes of C hkl in diamond/C (space group type Fm " 3 3m).

Discussion
GALOCS is an easy and illustrative way to find planar gaps in arbitrary crystal structures. Although we do not claim a oneto-one correspondence between the size of these gaps and the cleavage energies, some correlation between them exists. Specifically, the calculations in Section 3.1 suggest that the most prominent gap is seen along (111) and the next most prominent is seen along (110) Bloss & Gibbs (1963), when crushed, quartz cleaves most readily along (101)/(011) planes and next most readily along (112) planes. All these planes appear at the top of the list calculated by our algorithm.
Cleavage of LiNbO 3 crystals was recently investigated by Hirsh et al. (2020), featuring (012), (010) and (116) cleavage planes [the cleavage energies of (012) and (010) were measured]. While all these planes appear on the most prominent planar-gap list in Table 7, there is some disagreement between the measured cleavage energy and the gap width. Specifically (see Table 14), while the measured cleavage energy of (010) is $40% smaller than that of (012), the gap of (010) is four times narrower than (012). It is important to reiterate that this type of disagreement is expected, while the presence of these planes on the list is the simplified GALOCS scheme main goal (since most of the planes do not exhibit any gaps at all).
Still, in order to investigate the matter deeper, we inspected the À(hkl, D) for both (012) and (010)

Table 14
Cleavage energies in LiNbO 3 (obtained from an experiment) against the length of the structural gaps calculated in this work. (012) Cleavage-energy experiment (Hirsh et al., 2020) ðJ=m 2 Þ 2.2 1.3 Gap length (this work) (Å ) 1.55 0.39 that the À(012, D) curve has a well, which separates the gap into two 'valleys'. This suggests that the chosen threshold value y = 0.75 results in overestimation of the actual gap size for (012) planes. Such close inspection of À(hkl, D) is therefore recommended for all planes that are selected as candidate cleavage planes, and can be carried out by using the GALOCS package. Additionally, it is possible to recalculate all the C hkl values using different y thresholds if necessary. The algorithm can be improved by adding dummy atoms into the structures (e.g. to mimic chemical bonds). Alternatively, it may implement advanced models of electron density, which, in turn, is obtained by DFT calculations or results from a multipole model refinement of X-ray diffraction intensities (Hansen & Coppens, 1978). Additionally, atomic electron densities may be convoluted with their PDFs; this can be particularly valuable for disordered materials with high atomic displacement parameters. In general, the SFF (x 1 , x 2 , x 3 ) can be customized, e.g. by adding 'sticks' (additional electron densities) along the bond lines and removing the atoms themselves. This way the algorithm and the program will be capable of counting the number of chemical bonds intersected by the planes.

Conclusions
We developed geometrical algorithm GALOCS to locate the planes that expose the most prominent planar gaps in crystal structures. Such planes are listed as candidate cleavage planes, to be explored experimentally or by using density-functionaltheory calculations. GALOCS implements some known generalized space-filling function (e.g. superposition of atomic electron densities or the probability density function). It calculates the average values of this function within specific planes and as a function of plane depth with respect to the unit-cell origin. We also provided detailed illustration of the algorithm for silicon, quartz and LiNbO 3 where a clear correlation between our calculations and existing experimental studies of cleavage energy is present.

Funding information
The following funding is acknowledged: Israel Science Foundation (grant No. 1561/18).