research papers\(\def\hfill{\hskip 5em}\def\hfil{\hskip 3em}\def\eqno#1{\hfil {#1}}\)

IUCrJ
Volume 8| Part 5| September 2021| Pages 793-804
ISSN: 2052-2525

Geometrical prediction of cleavage planes in crystal structures

crossmark logo

aDepartment of Materials Science and Engineering, Tel Aviv University, Wolfson Building for Mechanical Engineering, Tel Aviv, 6997801, Israel, and bSchool of Mechanical Engineering, Tel Aviv University, Wolfson Building for Mechanical Engineering, Tel Aviv, 6997801, Israel
*Correspondence e-mail: gorfman@tauex.tau.ac.il

Edited by P. Lightfoot, University of St Andrews, United Kingdom (Received 26 May 2021; accepted 13 July 2021; online 20 August 2021)

Cleavage is the ability of single crystals to split easily along specifically oriented planes. This phenomenon is of great interest for materials' scientists. Acquiring the data regarding cleavage is essential for the understanding of brittle fracture, plasticity and strength, as well as for the prevention of catastrophic device failures. Unfortunately, theoretical calculations of cleavage energy are demanding and often unsuitable for high-throughput searches of cleavage planes in arbitrary crystal structures. A simplified geometrical approach (GALOCS = gaps locations in crystal structures) is suggested for predicting the most promising cleavage planes. GALOCS enumerates all the possible reticular lattice planes and calculates the plane-average electron density as a function of the position of the planes in the unit cell. The assessment of the cleavage ability of the planes is based on the width and depth of planar gaps in crystal structures, which appear when observing the planes lengthwise. The method is demonstrated on two-dimensional graphene and three-dimensional silicon, quartz and LiNbO3 structures. A summary of planar gaps in a few more inorganic crystal structures is also presented.

1. Introduction

Cleavage is the ability of single crystals to split easily along specifically oriented planes (Hurlbut & Klein, 1977[Hurlbut, C. S. Jr. & Klein, C. (1977). Manual of Mineralogy (after J. D. Dana), 19th ed. New York: John Wiley and Sons Inc..]). 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[Tutton, A. E. H. (1922). Crystallography and Practical Crystal Measurement, 2nd ed. London: Macmillan & Co. Ltd.]; Authier, 2013[Authier, A. (2013). Early Days of X-ray Crystallography. Oxford University Press.]). Cleavage is the most striking example of anisotropy of physical properties (Nye, 1985[Nye, J. F. (1985). Physical Properties of Crystals: Their Representation by Tensors and Matrices. Oxford: Clarendon Press.]). Finally, cleavage is the subject of many materials science oriented research (Gilman, 1960[Gilman, J. J. (1960). J. Appl. Phys. 31, 2208-2218.]; Lawn, 1974[Lawn, B. R. (1974). Mater. Sci. Eng. 13, 277-283.]): 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[Lawn, B., Lawn, T. & Wilshaw, R. (1993). Fracture of Brittle Solids. Cambridge University Press.]).

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[Spearing, S. M. (2000). Acta Mater. 48, 179-196.]). 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[Ramaseshan, S. (1946). Proc. Indian Acad. Sci. A, 24, 114.]). However, it becomes impractical for larger (e.g. ternary) structures. Accurate calculations of cleavage energies are demanding (Zhang et al., 2007[Zhang, S., Zhu, T. & Belytschko, T. (2007). Phys. Rev. B, 76, 094114.]; Bitzek et al., 2015[Bitzek, E., Kermode, J. R. & Gumbsch, P. (2015). Int. J. Fract, 191, 13-30.]) and should involve e.g. density functional theory (DFT) calculation of the surface energies (Ong et al., 2013[Ong, S. P., Richards, W. D., Jain, A., Hautier, G., Kocher, M., Cholia, S., Gunter, D., Chevrier, V. L., Persson, K. A. & Ceder, G. (2013). Comput. Mater. Sci. 68, 314-319.]; Tran et al., 2016[Tran, R., Xu, Z., Radhakrishnan, B., Winston, D., Sun, W., Persson, K. A. & Ong, S. P. (2016). Sci. Data, 3, 160080.]), which is half of the cleavage energy [for the case of slowly propagating cracks (Griffith, 1921[Griffith, A. A. (1921). Philos. Trans. R. Soc. A, 221, 163-198.])].

It is also possible to study cleavage experimentally (Cramer et al., 2000[Cramer, T., Wanner, A. & Gumbsch, P. (2000). Phys. Rev. Lett. 85, 788-791.]; Field, 1971[Field, J. E. (1971). Contemp. Phys. 12, 1-31.]; Jaccodine, 1963[Jaccodine, R. J. (1963). J. Electrochem. Soc. 110, 524.]; Lawn et al., 1993[Lawn, B., Lawn, T. & Wilshaw, R. (1993). Fracture of Brittle Solids. Cambridge University Press.]; Michot, 1987[Michot, G. (1987). Surf. Sci. 186, L561-L567.]; Sherman et al., 2008[Sherman, D., Markovitz, M. & Barkai, O. (2008). J. Mech. Phys. Solids, 56, 376-387.]). 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[Gilman, J. J. (1960). J. Appl. Phys. 31, 2208-2218.]; Gleizer & Sherman, 2014[Gleizer, A. & Sherman, D. (2014). Int. J. Fract, 187, 1-14.]; Sherman & Gleizer, 2014[Sherman, D. & Gleizer, A. (2014). Mater. Perform. Charact. 3, 189-213.]; Hirsh et al., 2020[Hirsh, Y., Gorfman, S. & Sherman, D. (2020). Acta Mater. 193, 338-349.]). 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[Momma, K. & Izumi, F. (2011). J. Appl. Cryst. 44, 1272-1276.])], 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.

2. 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.

2.1. 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:

[\rho \left({\bf r} \right) = \mathop \sum \limits_{m\mu } {\rho _\mu }({\bf r} - {{\bf R}_m} - {{\bf R}_\mu }). \eqno(1)]

Here, Rm = umiai (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, umi are arbitrary integers and ai are the lattice basis vectors i = 1, 2, 3. The vectors [{{\bf R}_\mu } = {x_{\mu i}}{{\bf a}_i}] (0 ≤ xμi < 1) list the positions of all the atoms in a unit cell. [{\rho _\mu }({\bf 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[Su, Z. & Coppens, P. (1998). Acta Cryst. A54, 357-357.]) or atomic invarioms (Dittrich et al., 2004[Dittrich, B., Koritsánszky, T. & Luger, P. (2004). Angew. Chem. Int. Ed. 43, 2718-2721.]). All the atomic electron densities are normalized: [\int\! {\rho _\mu }({\bf r} ){\rm d}{\bf r} = {Z_\mu }] [[{Z_\mu }] 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[Coppens, P. (1997). X-ray Charge Densities and Chemical Bonding. Oxford: International Union of Crystallography/Oxford University Press.]):

[{p_\mu }\left({\bf u} \right) = {{{{\left| {{\sigma ^{ - 1}_{\mu }}} \right|}^{1/2}}} \over {{{(2\pi)}^{3/2}}}}\exp \left({ - {1 \over 2}{{\bf u}^T}\hat{\sigma} _\mu ^{ - 1}{\bf u}} \right). \eqno(2)]

Here, [{\hat{\sigma}}_{\mu }] 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 (Uij or βij) (Trueblood et al., 1996[Trueblood, K. N., Bürgi, H.-B., Burzlaff, H., Dunitz, J. D., Gramaccioli, C. M., Schulz, H. H., Shmueli, U. & Abrahams, S. C. (1996). Acta Cryst. A52, 770-781.]; Coppens, 1997[Coppens, P. (1997). X-ray Charge Densities and Chemical Bonding. Oxford: International Union of Crystallography/Oxford University Press.]; Tsirelson & Ozerov, 1996[Tsirelson, V. G. & Ozerov, R. P. (1996). Electron Density and Bonding in Crystals: Principles, Theory and X-ray Diffraction Experiments in Solid State Physics and Chemistry. 1st ed. Boca Raton: CRC Press.]). If [{\hat \sigma _\mu }] is represented by a 3 × 3 matrix then uT and u are the rows or the column vectors, respectively. The SFF ρ(r) can now be calculated as

[\rho \left({\bf r}\right) = \sum _{m\mu }{W}_{\!\!\mu }\, {p}_{\!\mu }(r-{{\bf R}}_{m}-{{\bf R}}_{\mu }). \eqno(3)]

Here, [{W_\mu }] describe the `weights' of every atom in the cumulative SFF. [{W_\mu } = 1] yields an SFF in the form of atomic density (it disregards the sort of atoms involved). Setting [{W_\mu }] to the atomic masses defines the SFF as the mass density. Finally, setting [{W_\mu }] to [{Z_\mu }] 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 + Rm). Therefore, it is sufficient to calculate the values of ρ(r) inside the unit cell only. Because both [{p_\mu }({\bf r})] and [{\rho _\mu }({\bf r})] are negligible when |r| > rmax ([{r_{\rm max}}\simeq1\AA]), only the limited number of atoms whose centres are inside the unit cell or within rmax from its borders are included in the sums (1)[link] or (3)[link].

2.2. 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

[\Gamma \left({\bf n},D\right) = \langle \rho \left({\bf r}\right)\rangle {}_{{\bf r}{\bf n} = D}. \eqno(4)]

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) = dhkl. The latter is the inverse length of the primitive reciprocal lattice vector, [{{\bf B}_{hkl}} = {h_i}{\bf a}_i^*] {hi or h, k, l are coprime integers/Miller indices of the plane [the Miller indices are coprime if a primitive unit cell is used (Nespolo, 2015[Nespolo, M. (2015). J. Appl. Cryst. 48, 1290-1298.])] and [{\bf a}_i^*] are the basis vectors of the reciprocal lattice [{\bf a}_{i}{\bf a}_{j}^{{*}} = {\delta }_{ij}]} where Bhkln:

[{d_{hkl}} = {\left| {{B_{hkl}}} \right|^{ - 1}}. \eqno(5)]

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 A1, A2, A3 such that Ai and ai are the bases of the same crystal lattice, but the vectors A1 and A2 are parallel to the plane of interest (hkl) and A3 connects two adjacent lattice planes (A3ihi = 1). The number-theoretical algorithms for such a transformation are described elsewhere (Gorfman, 2020[Gorfman, S. (2020). Acta Cryst. A76, 713-718.]) 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 X1, X2 and X3 such that r = XiAi = xiai and expressing the SFF ρ(x1, x2, x3) as ρhkl(X1, X2, X3) yields

[\Gamma \left(hkl, {X}_{3}\right) = \langle {\rho }_{hkl}\left({X}_{1},{X}_{2},{X}_{3}\right)\rangle {}_{0\le {X}_{1},{X}_{2}\lt 1} \,\,\,\,{X}_{3}\in [0,1]. \eqno(6)]

The cleavage planes are likely to appear where Γ(hkl, X3) drops to the lowest possible values. Averaging of Γ(hkl, X3) over X3 would yield the average of the SFF [\langle\Gamma {({hkl,{X_3}})\rangle_{{X_3}}} = {\rho _0}] over the entire unit cell. Let us introduce the constant threshold ρS = yρ0 (typically y = 0.75) such that all the X3 values where Γ(hkl, X3) < ρS are considered as planar gaps. We define the maximum effective width of a planar gap as

[{C_{hkl}} = {1 \over {{\rho _0}}}{d_{hkl}}{{\Delta}}{X_3}\langle{\rho _S} - \Gamma {(hkl,{X_3})\rangle_{{{\Delta}}{X_3}}}. \eqno(7)]

Here ΔX3 is the length of the longest continuous range of X3 values where Γ(hkl, X3) < ρS and [\langle{\rho _S} - \Gamma {(hkl,{X_3})\rangle_{{{\Delta}}{X_3}}}] is the average gap depth in this range. The denominator ρ0 is introduced in order to reduce the dimension of Chkl to Å. Fig. 1[link] explains definition (7)[link] graphically. It shows an arbitrary Γ(hkl, D) with the most prominent gap between two vertical dashed lines.

[Figure 1]
Figure 1
A graphical illustration of the definition of Chkl according to equation (7)[link]. 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)[link] is equal to the area inside the coloured region.

Definition (7)[link] suggests that Chkl must be centrosymmetric. Indeed, the sets of Miller indices (hkl) and [(\bar h\bar k\bar l)] define the same lattice planes. The only difference between them concerns the direction of the plane movement with the increasing value of X3. Specifically, [\Gamma (\bar{h} \bar{k} \bar{l}, {X}_{3}) = \Gamma (hkl, -{X}_{3})], meaning that [\Gamma (\bar{h} \bar{k} \bar{l}, {X}_{3})] has the same minima, the same average value and the same width of the structural gap as Γ(hkl, X3). Accordingly, [{C}_{hkl} = {C}_{ \bar{h} \bar{k} \bar{l}}], so the symmetry of the Chkl is defined by one of the eleven Laue classes.

2.3. The flow chart of the GALOCS algorithm

GALOCS realizes the following steps:

(i) Calculating electron density ρ(x1, x2, x3) (or any other SFF in future releases) inside a unit cell (0 ≤ xi < 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 dhkl > dmin.

(iii) For each (hkl) transforming the basis vectors aiAi using the MULDIN algorithm (Gorfman, 2020[Gorfman, S. (2020). Acta Cryst. A76, 713-718.]), and expressing ρ(x1, x2, x3) as ρhkl(X1, X2, X3) (0 ≤ Xi < 1).

(iv) Calculating Chkl according to equation (7)[link], and drawing the directional dependence of Chkl 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 ρ(x1, x2, x3) 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[Su, Z. & Coppens, P. (1998). Acta Cryst. A54, 357-357.]).

(c) Drawing Γ(hkl, D) curves for any chosen hkl values.

(d) Calculation of Chkl (according to equation 7[link] with user-defined clearance threshold, typically y = 0.75) for the set of lattice planes with interplanar distance above some user-defined value dmin.

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.

2.4. Two-dimensional illustration

Fig. 2[link] illustrates the GALOCS output for a two-dimensional graphene structure. The conventional unit cell is based on the vectors a1 and a2 such that [{a}_{1} = {a}_{2} = a = 2.461\, \AA, \alpha =] [120^\circ]. The SFF was defined by the superposition of electron densities of isolated carbon atoms. Fig. 2[link](a) shows the transformed unit cells [based on the vectors A1 and A2 with A1 being parallel to (hk)] for the cases of (10), (11) and (12) planes. Fig. 2[link](b) illustrates the dependence of Chk on the reciprocal-space directions using the polar plot. The spikes extend in the directions of the reciprocal lattice vectors Bhk, their lengths are proportional to the corresponding Chk values. It shows that the longest gap in graphene is present for {10}H and {11}H families of planes [the notation [{\{{hk}\}_{\rm 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 Chk (see Table 1[link] for the numerical values).

Table 1
Magnitudes of Chk in graphene

(hk) (10) (11) (12)
[{C}_{hk} \,(\AA)] 0.54 0.40 0.03
[Figure 2]
Figure 2
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)[link] 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 A1 parallel to the (hk) planes. (b) The dependence of the Chk on the direction nBhk. The Miller indices corresponding to each of the directions are marked explicitly. (c) The dependences of Γ(10, X3), Γ(11, X3) and Γ(12, X3) on D = dhkX3 (vertically shifted relatively to each other for clarity).

3. Examples

Here, we demonstrate the implementation of the algorithm for the case of three inorganic structures (Si, LiNbO3 and SiO2).

3.1. Silicon/diamond structural type

Fig. 3[link] projects the structure of silicon along the crystallographic [001] and [110] directions. The structure has one symmetry-independent atom sitting at the standard origin 1 of the space group [Fd \overline{3}m]. Each atom has four nearest neighbours at the corners of a regular tetrahedron. Table 2[link] summarizes the relevant structural information.

Table 2
The structure of silicon used for the calculations

Space group [Fd\bar 3m] (No. 227)
Laue class [m\bar 3m]
Number of atoms per unit cell 8
Lattice parameter(s) (Å) a = 5.43

Asymmetric unit

Atom name Charge Position Wyckoff letter and multiplicity Local symmetry
Si 14 [000] 8a [\bar 43m]
[Figure 3]
Figure 3
The structure of silicon viewed along (a) [001] and (b) [110] crystallographic directions. The images were produced using the VESTA program (Momma & Izumi, 2011[Momma, K. & Izumi, F. (2011). J. Appl. Cryst. 44, 1272-1276.]).

The calculation of Chkl over the lattice planes with the interplanar distance above [{d}_{\rm min} = 0.2 \AA] involved 17 330 primitive reciprocal lattice vectors enclosed in a corresponding reciprocal-space sphere of radius [{d}_{\rm min}^{-1} = 5{\AA}^{-1}]. This set included 446 symmetry-independent planes (with respect to the operations of the point-symmetry group [m \overline{3}m]). The cut-off dmin is justified by the fact that Chkl drops to zero for those planes that have small interplanar distances. Fig. 4[link] illustrates this statement by showing Chkl as a function of the length of the primitive reciprocal lattice vector ([{B_{hkl}} = d_{hkl}^{ - 1}]). Prominent gaps are present among low Bhkl (high dhkl) 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[Nespolo, M. (2015). J. Appl. Cryst. 48, 1290-1298.]).

[Figure 4]
Figure 4
The Bhkl dependence of Chkl values in silicon. Here Bhkl is the length of the corresponding reciprocal lattice vector (inverse to dhkl). 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.

Figs. 5[link] and 6[link] illustrate the striking anisotropy of Chkl 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 Chkl values. The projection is viewed along the zone axis Auvw = [uvw] = ua1 + va2 + wa3 {[[1\bar 10]]/[111] in Figs. 5[link](a)/6[link](a), respectively}. Figs. 5[link](b) and 6[link](b) show polar plots of Chkl for all the reciprocal lattice directions Bhkl in the zone (i.e. such that Bhkl · Auvw = 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 Auvw. Cleavage is initiated at a small pre-crack [see e.g. (Hirsh et al., 2020[Hirsh, Y., Gorfman, S. & Sherman, D. (2020). Acta Mater. 193, 338-349.])] and at a uniaxial stress [also known as Mode I stress (Lawn et al., 1993[Lawn, B., Lawn, T. & Wilshaw, R. (1993). Fracture of Brittle Solids. Cambridge University Press.])] along the direction n. Table 3[link] lists numerical values for the cleavage ability of the most prominent planes.

Table 3
Magnitudes of Chkl in silicon

(hkl) (111) (110) (100) (311) (211) (331) (511)
[{C}_{hkl} \, (\AA)] 1.025 0.676 0.289 0.285 0.149 0.121 0.047
[Figure 5]
Figure 5
Anisotropy of Chkl values in silicon. (a) The stereographic projection viewed along the [[1\bar 10]] crystallographic direction displays the false-colour map of the Chkl (dependence of Chkl 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 Chkl/max(Chkl) according to the colour bar on the left. (b) The polar diagram of the cleavage ability for all the planes in the [[1\bar 10]] zone. The outer circle corresponds to the highest cleavage ability (C111 = 1.025 Å for the case of silicon).
[Figure 6]
Figure 6
Same as Fig. 5[link] but for the [111] direction in silicon. The outer circle in (b) corresponds to the highest cleavage ability C111 = 1.025 Å.

The calculations predict {111}C and {110}C as the most prominent cleavage planes in silicon [the notation [{\{{hkl}\}_{\rm C}}] stands for the set of planes that are symmetry equivalent to (hkl) with respect to the operations of the cubic crystal class [m\bar 3m]]. These planes are well known from real cleavage experiments on single crystals of silicon (Gleizer et al., 2014[Gleizer, A., Peralta, G., Kermode, J. R., De Vita, A. & Sherman, D. (2014). Phys. Rev. Lett. 112, 115501.]).

The output is further illustrated in Figs. 7[link] and 8[link]. Fig. 7[link] shows the corresponding sections of electron density by (111) planes and the Γ(n, D) function where the gap is seen clearly. Fig. 8[link] 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 7]
Figure 7
An illustration of the structural-gap location for the case of the (111) plane in silicon. The parts (a)–(c) show three electron-density sections by the (111) plane at three various positions of the plane in the unit cell, while (d) shows the Γ(111, D) function with the vertical lines at the positions corresponding to the electron-density sections in (a)–(c).
[Figure 8]
Figure 8
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[link]).

3.2. Quartz (SiO2)

Fig. 9[link] projects the structure of quartz along [100] and [001] directions, while Table 4[link] summarizes the relevant structural information (Levien et al., 1980[Levien, L., Prewitt, C. T. & Weidner, D. J. (1980). Am. Mineral. 65, 920-930.]). Right-handed quartz crystallizes in the space group type P3221 (No. 154) and corresponds to the Laue class [\overline{3}m]. The structure has two symmetry-independent atoms (Si and O). Enumeration of the lattice planes with the interplanar distance above [{d}_{\rm min} = 0.2 \AA] 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[link]–12[link][link] and Table 5[link] show the results of the calculations (as in chapter 3.1[link]). Fig. 10[link] illustrates the dependence of Chkl on the length of the reciprocal lattice vectors. Figs. 11[link] and 12[link] illustrate the anisotropy of Chkl using stereographic projections along [100] and [001] zone axes [defining the normal to X and Z cuts of quartz wafers (Brainerd et al., 1949[Brainerd, J. G., Jensen, A. G., Cumming, L. G., Batcher, R. R., Begun, S. G., Black, H. S., Grown, G. M., Burrows, C. R., Busignies, H., Cady, W. G. et al. (1949). Proc. IRE, 37, 1378-1395.])]. Table 5[link] lists the numerical values of Chkl of seven planes with the most prominent gaps. Figs. 13[link] and 14[link] (organized in the same way as Figs. 7[link] and 8[link]) illustrate the output for the case of (011) planes in quartz.

Table 4
The structure of quartz used for the calculations

Space group P3221 (No. 154)
Laue class [\bar 3m]
Number of atoms in the unit cell 9
Lattice parameters (Å) a = 4.91, c = 5.40

Asymmetric unit

Atom name Charge Position Wyckoff letter and multiplicity Local symmetry
Si 14 [0.47 00] 3a 2
O 8 [0.41 0.27 0.12] 6b 1

Table 5
Magnitudes of Chkl in quartz

(hkl) (011) (101) (203) (112) (031) (110) (010) (102)
[{C}_{hkl} \, (\AA)] 0.714 0.483 0.213 0.209 0.170 0.165 0.138 0.114
[Figure 9]
Figure 9
The structure of quartz viewed along [100] and [001] crystallographic directions. The images were produced using VESTA (Momma & Izumi, 2011[Momma, K. & Izumi, F. (2011). J. Appl. Cryst. 44, 1272-1276.]).
[Figure 10]
Figure 10
Same as Fig. 4[link] but for the case of quartz.
[Figure 11]
Figure 11
Same as Fig. 5[link] but for the [100] direction in quartz. The outer circle in (b) corresponds to the highest cleavage ability C011 = 0.714 Å.
[Figure 12]
Figure 12
Same as Fig. 5[link] but for the [001] direction in quartz. The outer circle in (b) corresponds to the highest cleavage ability C011 = 0.714 Å.
[Figure 13]
Figure 13
An illustration of the gap location for the case of the (011) plane in quartz. Parts (a)–(c) show three electron-density sections by the (011) plane at three various positions of the plane in the unit cell, while (d) shows the Γ(011, D) function with the vertical lines at the positions corresponding to the electron-density sections in (a)–(c).
[Figure 14]
Figure 14
The structure of quartz viewed along the (011) plane. The [011]* reciprocal lattice direction is horizontal.

3.3. Lithium niobate (LiNbO3)

LiNbO3 crystallizes in the space group type R3c (Weis & Gaylord, 1985[Weis, R. S. & Gaylord, T. K. (1985). Appl. Phys. A, 37, 191-203.]; Weigel et al., 2020[Weigel, T., Funke, C., Zschornak, M., Behm, T., Stöcker, H., Leisegang, T. & Meyer, D. C. (2020). J. Appl. Cryst. 53, 614-622.]), Laue class [\bar 3m]. Table 6[link] summarizes the relevant structural information. The structure has three symmetry-independent atoms (Li, Nb and O).

Table 6
The structure of LiNbO3 used for the calculations (Weigel et al., 2020[Weigel, T., Funke, C., Zschornak, M., Behm, T., Stöcker, H., Leisegang, T. & Meyer, D. C. (2020). J. Appl. Cryst. 53, 614-622.])

Space group R3c (No. 161)
Laue class [\bar 3m]
Number of atoms per unit cell 30
Lattice parameters (Å) a = 5.15160 (10) c = 13.8690 (6)

Asymmetric unit

Atom name Charge Position Wyckoff Local symmetry
Li 3 [00  0.3014] 6a 3
Nb 41 [00  0.0199] 6a 3
O 8 [\left[0.0488\,\,\, 0.3438\,\,\, {{1}\over{12}}\right]] 18b 1

Enumeration of the lattice planes with the interplanar distance exceeding [{d}_{\rm min} = 0.2 \AA] results in the generation of 46 340 reciprocal lattice vectors, where 4108 of them are symmetry independent (with respect to the point-symmetry group [\bar 3m]). The calculation results are presented in Figs. 15[link]–19[link][link][link][link] and Table 7[link], 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 LiNbO3 crystals (Hirsh et al., 2020[Hirsh, Y., Gorfman, S. & Sherman, D. (2020). Acta Mater. 193, 338-349.]).

Table 7
Magnitudes of Chkl in LiNbO3

(hkl) (012) (104) (110) (001) (101) (116) (122) (010)
[{C}_{hkl}\, (\AA)] 1.556 0.945 0.776 0.572 0.562 0.482 0.392 0.390
[Figure 15]
Figure 15
Same as Fig. 4[link] but for the case of LiNbO3.
[Figure 16]
Figure 16
Same as Fig. 5[link] but for the [100] direction in LiNbO3. The outer circle in (b) corresponds to the highest cleavage ability C012 = 1.556 Å.
[Figure 17]
Figure 17
Same as Fig. 5[link] but for the [001] direction in LiNbO3. The outer circle in (b) corresponds to the highest cleavage ability C012 = 1.556 Å.
[Figure 18]
Figure 18
An illustration of the gap location for the case of the (012) plane in LiNbO3. Parts (a)–(c) show three electron-density sections by the (012) plane at three various positions of the plane in the unit cell, while (d) shows the Γ(012, D) function with the vertical lines at the positions corresponding to the electron-density sections in (a)–(c).
[Figure 19]
Figure 19
The structure of LiNbO3 viewed along the (012) plane. The [012]* reciprocal lattice direction is horizontal.

4. Further examples

Tables 8[link]–12[link][link][link][link] present a brief summary of the Chkl calculations for the structures of wurzite (AlN), fluorite (CaF2), diamond (C), pyrite (FeS2) and corundum (Al2O3). They are organized in the same way as Tables 3[link], 5[link] and 7[link]. All the examples are included in the GALOCS user manual in the supporting information.

Table 8
Magnitudes of Chkl in AlN [wurzite structural type, space group type P63mc, the structure is reported by Xu & Ching (1993[Xu, Y. N. & Ching, W. Y. (1993). Phys. Rev. B, 48, 4335-4351.])]

(hkl) (001) (010) ([1\bar 20)] (011) (013) ([\bar 122)] (012) (021)
[{C}_{hkl} \,(\AA)] 0.774 0.654 0.462 0.385 0.257 0.227 0.189 0.100

Table 9
Magnitudes of Chkl in CaF2 [fluorite structural type, space group type [Fm\bar 3m], the structure is reported by Speziale & Duffy, 2002[Speziale, S. & Duffy, T. S. (2002). Phys. Chem. Miner. 29, 465-472.])]

(hkl) (110) (111) (100) (311) (211) (331) (511) (310)
[{C}_{hkl}\, (\AA)] 0.718 0.648 0.338 0.285 0.182 0.141 0.067 0.040

Table 10
Magnitudes of Chkl in diamond/C (space group type [Fm\bar 3m])

(hkl) (111) (110) (311) (200) All the rest
[{C}_{hkl}\, (\AA)] 0.465 0.222 0.045 0.016 0

Table 11
Magnitudes of Chkl in pyrite [FeS2, space group type [Fm\bar 3m], the structure is reported by Ramsdell (1925[Ramsdell, L. S. (1925). Am. Mineral. 10, 281-304.])]

(hkl) (011) (113) (111) (102) (001) (012) (112) (115)
[{C}_{hkl}\, (\AA)] 0.378 0.299 0.249 0.154 0.151 0.128 0.057 0.040

Table 12
Magnitudes of Chkl in corundum [space group type [R\bar 3c], the structure is reported by Lewis et al. (1982[Lewis, J., Schwarzenbach, D. & Flack, H. D. (1982). Acta Cryst. A38, 733-739.])]

(hkl) (012) [(\bar 114)] (010) [(1\bar 20)] [(1\bar 2\bar 6)] [(2\bar 3\bar 4)] [(0\, \overline{1} \,10)] [(\bar 1\bar 13)]
[{C}_{hkl}\, (\AA)] 0.520 0.406 0.314 0.260 0.259 0.131 0.130 0.072

5. Discussion

GALOCS is an easy and illustrative way to find planar gaps in arbitrary crystal structures. Although we do not claim a one-to-one correspondence between the size of these gaps and the cleavage energies, some correlation between them exists. Specifically, the calculations in Section 3.1[link] suggest that the most prominent gap is seen along (111) and the next most prominent is seen along (110) planes. This result is proven by both experiments (Gleizer et al., 2014[Gleizer, A., Peralta, G., Kermode, J. R., De Vita, A. & Sherman, D. (2014). Phys. Rev. Lett. 112, 115501.]) and calculations (Pérez & Gumbsch, 2000[Pérez, R. & Gumbsch, P. (2000). Phys. Rev. Lett. 84, 5347-5350.]), suggesting that single crystals of silicon will indeed break most readily along (111) and next most readily along (110) planes. Specifically, these literature results imply that cleavage along (110) planes requires ∼20% more energy input than cleavage along (111) planes. Table 13[link] shows the numerical values of cleavage energies against the calculated parameters of the gaps.

Table 13
Cleavage energies of silicon (obtained from DFT calculations/experiment) against the length of the structural gaps calculated in this work

  (111) (110)
Cleavage-energy DFT (Pérez & Gumbsch, 2000[Pérez, R. & Gumbsch, P. (2000). Phys. Rev. Lett. 84, 5347-5350.]) (J/m2) 2.88 3.46
Cleavage-energy experiment (Gleizer et al., 2014[Gleizer, A., Peralta, G., Kermode, J. R., De Vita, A. & Sherman, D. (2014). Phys. Rev. Lett. 112, 115501.]) (J /m2) 2.2 ± 0.2 2.7 ± 0.3
Gap length (this work) (Å) 1.025 0.672

Cleavage of quartz single crystals is debated in the literature [see e.g. (White, 2006[White, J. S. (2006). Rocks & Minerals, 81, 389-391.]; Bloss & Gibbs, 1963[Bloss, F. D. & Gibbs, G. V. (1963). Am. Mineral. 48, 821-838.])]. We are unaware of any accurate calculations or precise measurements of the cleavage energies in quartz. Nonetheless, according to Bloss & Gibbs (1963[Bloss, F. D. & Gibbs, G. V. (1963). Am. Mineral. 48, 821-838.]), 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 LiNbO3 crystals was recently investigated by Hirsh et al. (2020[Hirsh, Y., Gorfman, S. & Sherman, D. (2020). Acta Mater. 193, 338-349.]), 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[link], there is some disagreement between the measured cleavage energy and the gap width. Specifically (see Table 14[link]), 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).

Table 14
Cleavage energies in LiNbO3 (obtained from an experiment) against the length of the structural gaps calculated in this work

  (012) (010)
Cleavage-energy experiment (Hirsh et al., 2020[Hirsh, Y., Gorfman, S. & Sherman, D. (2020). Acta Mater. 193, 338-349.]) (J/m2) 2.2 1.3
Gap length (this work) (Å) 1.55 0.39

Still, in order to investigate the matter deeper, we inspected the Γ(hkl, D) for both (012) and (010) planes. Fig. 20[link] shows 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 Chkl values using different y thresholds if necessary.

[Figure 20]
Figure 20
Γ(012, D) and Γ(010, D) curves for LiNbO3. The horizontal dashed line is drawn at the chosen threshold level of 0.75ρ0.

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[Hansen, N. K. & Coppens, P. (1978). Acta Cryst. A34, 909-921.]). 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 ρ(x1, x2, x3) 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.

6. 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-functional-theory 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 LiNbO3 where a clear correlation between our calculations and existing experimental studies of cleavage energy is present.

Supporting information


Funding information

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

References

First citationAuthier, A. (2013). Early Days of X-ray Crystallography. Oxford University Press.  Google Scholar
First citationBitzek, E., Kermode, J. R. & Gumbsch, P. (2015). Int. J. Fract, 191, 13–30.  CrossRef Google Scholar
First citationBloss, F. D. & Gibbs, G. V. (1963). Am. Mineral. 48, 821–838.  CAS Google Scholar
First citationBrainerd, J. G., Jensen, A. G., Cumming, L. G., Batcher, R. R., Begun, S. G., Black, H. S., Grown, G. M., Burrows, C. R., Busignies, H., Cady, W. G. et al. (1949). Proc. IRE, 37, 1378–1395.  Google Scholar
First citationCoppens, P. (1997). X-ray Charge Densities and Chemical Bonding. Oxford: International Union of Crystallography/Oxford University Press.  Google Scholar
First citationCramer, T., Wanner, A. & Gumbsch, P. (2000). Phys. Rev. Lett. 85, 788–791.  CrossRef PubMed CAS Google Scholar
First citationDittrich, B., Koritsánszky, T. & Luger, P. (2004). Angew. Chem. Int. Ed. 43, 2718–2721.  Web of Science CSD CrossRef CAS Google Scholar
First citationField, J. E. (1971). Contemp. Phys. 12, 1–31.  CrossRef Google Scholar
First citationGilman, J. J. (1960). J. Appl. Phys. 31, 2208–2218.  CrossRef CAS Google Scholar
First citationGleizer, A., Peralta, G., Kermode, J. R., De Vita, A. & Sherman, D. (2014). Phys. Rev. Lett. 112, 115501.  CrossRef PubMed Google Scholar
First citationGleizer, A. & Sherman, D. (2014). Int. J. Fract, 187, 1–14.  CrossRef CAS Google Scholar
First citationGorfman, S. (2020). Acta Cryst. A76, 713–718.  Web of Science CrossRef IUCr Journals Google Scholar
First citationGriffith, A. A. (1921). Philos. Trans. R. Soc. A, 221, 163–198.  Google Scholar
First citationHansen, N. K. & Coppens, P. (1978). Acta Cryst. A34, 909–921.  CrossRef CAS IUCr Journals Web of Science Google Scholar
First citationHirsh, Y., Gorfman, S. & Sherman, D. (2020). Acta Mater. 193, 338–349.  CrossRef CAS Google Scholar
First citationHurlbut, C. S. Jr. & Klein, C. (1977). Manual of Mineralogy (after J. D. Dana), 19th ed. New York: John Wiley and Sons Inc..  Google Scholar
First citationJaccodine, R. J. (1963). J. Electrochem. Soc. 110, 524.  CrossRef Google Scholar
First citationLawn, B., Lawn, T. & Wilshaw, R. (1993). Fracture of Brittle Solids. Cambridge University Press.  Google Scholar
First citationLawn, B. R. (1974). Mater. Sci. Eng. 13, 277–283.  CrossRef CAS Google Scholar
First citationLevien, L., Prewitt, C. T. & Weidner, D. J. (1980). Am. Mineral. 65, 920–930.  CAS Google Scholar
First citationLewis, J., Schwarzenbach, D. & Flack, H. D. (1982). Acta Cryst. A38, 733–739.  CrossRef ICSD CAS Web of Science IUCr Journals Google Scholar
First citationMichot, G. (1987). Surf. Sci. 186, L561–L567.  CrossRef CAS Google Scholar
First citationMomma, K. & Izumi, F. (2011). J. Appl. Cryst. 44, 1272–1276.  Web of Science CrossRef CAS IUCr Journals Google Scholar
First citationNespolo, M. (2015). J. Appl. Cryst. 48, 1290–1298.  Web of Science CrossRef CAS IUCr Journals Google Scholar
First citationNye, J. F. (1985). Physical Properties of Crystals: Their Representation by Tensors and Matrices. Oxford: Clarendon Press.  Google Scholar
First citationOng, S. P., Richards, W. D., Jain, A., Hautier, G., Kocher, M., Cholia, S., Gunter, D., Chevrier, V. L., Persson, K. A. & Ceder, G. (2013). Comput. Mater. Sci. 68, 314–319.  Web of Science CrossRef CAS Google Scholar
First citationPérez, R. & Gumbsch, P. (2000). Phys. Rev. Lett. 84, 5347–5350.  PubMed Google Scholar
First citationRamaseshan, S. (1946). Proc. Indian Acad. Sci. A, 24, 114.  CrossRef Google Scholar
First citationRamsdell, L. S. (1925). Am. Mineral. 10, 281–304.  CAS Google Scholar
First citationSherman, D. & Gleizer, A. (2014). Mater. Perform. Charact. 3, 189–213.  Google Scholar
First citationSherman, D., Markovitz, M. & Barkai, O. (2008). J. Mech. Phys. Solids, 56, 376–387.  CrossRef CAS Google Scholar
First citationSpearing, S. M. (2000). Acta Mater. 48, 179–196.  CrossRef CAS Google Scholar
First citationSpeziale, S. & Duffy, T. S. (2002). Phys. Chem. Miner. 29, 465–472.  Web of Science CrossRef CAS Google Scholar
First citationSu, Z. & Coppens, P. (1998). Acta Cryst. A54, 357–357.  CrossRef CAS IUCr Journals Google Scholar
First citationTran, R., Xu, Z., Radhakrishnan, B., Winston, D., Sun, W., Persson, K. A. & Ong, S. P. (2016). Sci. Data, 3, 160080.  CrossRef PubMed Google Scholar
First citationTrueblood, K. N., Bürgi, H.-B., Burzlaff, H., Dunitz, J. D., Gramaccioli, C. M., Schulz, H. H., Shmueli, U. & Abrahams, S. C. (1996). Acta Cryst. A52, 770–781.  CrossRef CAS Web of Science IUCr Journals Google Scholar
First citationTsirelson, V. G. & Ozerov, R. P. (1996). Electron Density and Bonding in Crystals: Principles, Theory and X-ray Diffraction Experiments in Solid State Physics and Chemistry. 1st ed. Boca Raton: CRC Press.  Google Scholar
First citationTutton, A. E. H. (1922). Crystallography and Practical Crystal Measurement, 2nd ed. London: Macmillan & Co. Ltd.  Google Scholar
First citationWeigel, T., Funke, C., Zschornak, M., Behm, T., Stöcker, H., Leisegang, T. & Meyer, D. C. (2020). J. Appl. Cryst. 53, 614–622.  CrossRef CAS IUCr Journals Google Scholar
First citationWeis, R. S. & Gaylord, T. K. (1985). Appl. Phys. A, 37, 191–203.  CrossRef Web of Science Google Scholar
First citationWhite, J. S. (2006). Rocks & Minerals, 81, 389–391.  CrossRef Google Scholar
First citationXu, Y. N. & Ching, W. Y. (1993). Phys. Rev. B, 48, 4335–4351.  CrossRef ICSD CAS Google Scholar
First citationZhang, S., Zhu, T. & Belytschko, T. (2007). Phys. Rev. B, 76, 094114.  CrossRef Google Scholar

This is an open-access article distributed under the terms of the Creative Commons Attribution (CC-BY) Licence, which permits unrestricted use, distribution, and reproduction in any medium, provided the original authors and source are cited.

IUCrJ
Volume 8| Part 5| September 2021| Pages 793-804
ISSN: 2052-2525