Crystal Structure and Mechanical Properties of ThBC2

Abstract

Thorium borocarbide compounds have fascinating physical properties and diverse structures, and hence have stimulated great interest. In this work, we determine the ground state structure of ThBC₂ by the unbiased structure prediction method based on first-principles calculations. The dynamical and elastic stabilities of our proposed ThBC₂ are verified by the calculations of phonon spectrum and elastic constants. To study the mechanical properties fundamentally, we estimated the elastic anisotropy of ThBC₂. The results show that the Young’s and shear moduli possess high degree of anisotropy. The ideal strength calculations reveal that ThBC₂ readily collapses upon applied stress due to small ideal strengths. The cleavage fracture probably occurs along the [111] direction while slip may easily appear along the [$$\overline{1}10$$] direction on the (111) plane for ThBC₂. In addition, we provide an atomic explanation for the different characteristics of the strain–stress curves under different strains.

1. Introduction

Due to the potential application in fission fuel, actinoid-metal borocarbides are of great interest, especially thorium borocarbide compounds [1,2,3,4,5,6,7]. Thorium borocarbide compounds have diversified stoichiometric proportions with different crystal structures due to the flexible B–C framework. Toth et al. [1] first investigated thorium borocarbide compounds, and observed four related compounds, including ThBC, ThB₂C, ThBC_2, and Th₂BC₂. The crystal structure of ThBC [2] shows a tetragonal symmetry with the P4₁22 space group, which contains tetrahedra and trigonal prisms of Th atoms, and two boron and two carbon atoms form isolated zigzag chains. ThB₂C [3] crystallizes in the rhombohedral space group ($R\overline{3}m$) with slightly puckered Th-metal layers, and boron atoms form B₆ hexagons connected by carbon atoms. Based on X-ray powder diffraction experiment, Rogl et al. reported the crystal structures of Th₃B₂C₃ [4] and Th₂B₂C₃ [5]. The monoclinic Th₃B₂C₃ belongs to the space group of P2/m, in which Th atoms form octahedra, tetrahedra, and trigonal prisms with B and C atoms occupying their centers, and non-metal atoms form isolated zigzag B–C chain fragments [4], which are similar to the B–C chains in ThBC [2]. The crystal structure of Th₂B₂C₃ consists of the infinite branched B–C chains and triangular prisms of Th atoms. Therefore, the flexible B–C framework leads to the complexity of the crystal structures for thorium borocarbide compounds. Consequently, the crystal structure of ThBC₂, observed by Toth et al. long ago [1], remains an open question in the past years.

Actually, the determinations of crystal structures of ternary borocarbide compounds would face many problems. First, the distributions of B and C atoms are so complex for a given crystal framework that it is hard to determine the lowest energy phases of ternary borocarbide compounds, which are called “coloring problems” in structural chemistry [8,9]. Second, from the point of view of the experiment, ternary borocarbide compounds contain heavy metal atoms, which can be readily distinguished by XRD or neutron-diffraction measurements. But for light B and C atoms, it is difficult to resolve their Wyckoff positions because the electronic and nuclear scattering cross-sections of B and C atoms are very similar [10]. Therefore, it is challenging to solve their crystal structures in experiments. To complement experiments, many crystal structure prediction techniques [11,12,13,14] based on the first-principles density functional method have been proposed recently. In particular, the particle swarm optimization (PSO) technique as implemented in the crystal structure analysis by particle swarm optimization (CALYPSO) [11] has been suggested to be an efficient tool to resolve the ground state of certain compound [15,16,17].

In the present work, we performed the PSO technique combined with first-principles calculations to determine the ground state structure of ThBC₂. The stability of ThBC₂ is verified by the phonon spectrum and elastic constants. To show the mechanical properties, we have systematically studied the elastic anisotropies of ThBC_2. In addition, the failure modes under tensile and shear strains are investigated by the ideal strength calculations.

2. Computational Methods

The crystal structure of ThBC₂ was predicted using the PSO technique as implemented in the CALYPSO package [11]. To improve the prediction efficiency and avoid possible breakdown of structure optimization, the initial minimum inter-atomic distance between different types of atoms are 0.6 times the length of the sum of atomic covalent radii for the random structures. The proportion of the structures generated by PSO technique is set to 0.6. Other parameters are set as the default values of the CALYPSO package. All the self-consistent field calculations and structural relaxations were carried out within the framework of density functional theory using the Vienna ab initio simulation package (VASP) [18]. The generalized gradient approximation (GGA) with the Perdew–Burke–Ernzerhof functional was employed to describe exchange–correlation interaction [19]. The interactions between electron and core were described by the frozen-core projector augmented wave (PAW) method [20], where the 6s²6p⁶6d²7s², 2s²2p¹, and 2s²2p² are valence electrons for Th, B, and C atoms, respectively. The PAW pseudo potential has been successfully used to study the compounds including Th atoms, such as thorium dicarbide [21], thorium dioxide [22], thorium dihydride [23], and thorium mononitrides [24]. The cutoff energy of 600 eV for plane wave and the 9 × 9 × 9 k-point meshes with Monkhorst–Pack method [25] were adopted, which can ensure that the total energy was converged to 1 meV per atom. All the structures were relaxed until the forces on each atom are below 0.001 eV/Å. To verify the dynamical stability, we performed the finite displacement method to calculate the phonon spectrum of ThBC₂ with 2 × 2 × 2 supercell of primitive cell and 4 × 4 × 4 k-point meshes as implemented in PHONOPY package [26]. The elastic constants were estimated by the strain–stress method described in our previous work [27], and a set of homogeneous strains with maximum strain value of 0.1% is applied. The elastic moduli and Pugh’s ratio were calculated using the Voigt–Reuss–Hill approximation [28]. To estimate the ideal tensile and shear strengths, the conventional cell of ThBC₂ is deformed incrementally and continuously along the direction of the applied strain. Then, the atomic basis vectors orthogonal to the applied strain and the positions of all atoms are relaxed simultaneously. As a result, the strain–stress curves for deducing the ideal strength are obtained [29,30]. Considering that the spin-orbit coupling (SOC) plays an important role on the compounds including heavy atoms [31], we tested the influence of SOC on the crystal structure and elastic constants of ThBC₂. The calculated results [see the Supplementary Material (SM)] imply that SOC has little impact on the structure and mechanical properties of ThBC₂. Therefore, we can ignore the SOC effect in the following. For orthorhombic structure, the Young’s modulus E along the [hkl] direction is described by the expression [32]

(1)$E^{-1}=s_{11}{\alpha }^4+s_{22}{\beta }^4+s_{33}{\gamma }^4+\left(2s_{12}+s_{66}\right){\alpha }^2{\beta }^2+\left(2s_{23}+s_{44}\right){\beta }^2{\gamma }^2+\left(2s_{13}+s_{55}\right){\alpha }^2{\gamma }^2,$

where s_ij are the elastic compliance constants of orthorhombic symmetry, α, β, and γ denote the direction cosines of the [hkl] direction. The shear modulus G along the arbitrary [uvw] directions on the (hkl) shear plane can be written as [32].

(2)$\begin{array}{ll}{G}^{-1}=& 4s_{11}{\alpha }_1^2{\alpha }_2^2+4s_{22}{\beta }_1^2{\beta }_2^2+4s_{33}{\gamma }_1^2{\gamma }_2^2+8s_{12}{\alpha }_1{\alpha }_2{\beta }_1{\beta }_2+8s_{23}{\beta }_1{\beta }_2{\gamma }_1{\gamma }_2+8s_{13}{\alpha }_1{\alpha }_2{\gamma }_1{\gamma }_2\\ & +s_{44}{\left({\beta }_1{\gamma }_2+{\beta }_2{\gamma }_1\right)}^2+s_{55}{\left({\alpha }_1{\gamma }_2+{\alpha }_2{\gamma }_1\right)}^2+s_{66}{\left({\alpha }_1{\beta }_2+{\alpha }_2{\beta }_1\right)}^2,\end{array}$

where α₁, β₁, and γ₁ are the direction cosines of the [uvw] direction, while α₂, β₂, and γ₂ are the direction cosines of the normal direction of the (hkl) shear plane. The anisotropy of Young’s modulus and shear modulus are discussed based on above two expressions.

3. Results and Discussion

To determine the ground state structure of ThBC₂, we performed the structural search by the unbiased PSO technique [11] with the cell up to 8 formula units (f.u.). As a result, we obtained the ground state structure of ThBC₂ according to comparing their total energies. As shown in Figure 1a,b, the predicted crystal structure of ThBC₂ belong to the space group of Immm (no. 71) within orthorhombic symmetry (16 atoms/cell), in which Th atoms occupy the 4e (0.26226, 0, 0) Wyckoff position, B atoms occupy 4j (1/2, 0, 0.36170) Wyckoff position, and C atoms occupy 8l (0, 0.17557, 0.26910) Wyckoff position, respectively. The optimized lattice constants of the conventional cell are a = 7.2403 Å, b = 3.9978 Å, and c = 6.7686 Å. Clearly, we can see that the B and C atoms form the armchair-like B–C chains along (010) direction (see Figure 1c). These armchair-like B–C chains are separated by the adjacent Th atomic layers. Note that the armchair-like B–C chains had been uncovered in other thorium borocarbide compounds, such as ThBC [2,7]. One Th atom and four C atoms can constitute a square pyramid (see Figure 1d), which serves as a basic building-block of ThBC₂. The average B–C bond length of ThBC₂ is 1.57 Å, which is shorter than those of t-BC₄ (1.66 Å) [27,33] and d-BC₃ (1.61 Å) [34], but is slightly longer than those of ThBC (1.54 Å) [2] and ThB₂C (1.49 Å) [3]. The interatomic distance between adjacent B atoms is 1.87 Å, which is also larger than those of ThBC (1.77 Å) [2] and ThB₂C (1.85 Å) [3]. In addition, to facilitate the possible experimental synthesis, we discuss the metastable ThBC₂ in the SM, and the crystal structures of the most stable form and other four predicted forms are available in the SM.

To verify the stability, we first calculated the phonon spectrum of ThBC₂. As shown in Figure 2a, it is found that there is no imaginary frequency in the Brillouin zone along the high-symmetry direction, which indicates the dynamical stability of ThBC₂. Moreover, due to the large differences of atomic masses among Th, B, and C atoms, the optical branches and acoustic branches of ThBC₂ are separated distinctly, leading to large acoustic-optical branch gap. Meanwhile, the high-frequency optical phonon bands are also discrete. The formation energy can reflect the relative stability of the compound and whether the compound can be readily synthesized experimentally. Here, we choose the face-centered cubic thorium, rhombohedral α-boron and graphite as the reference phases. Therefore, the formation energy $\Delta E_f$ of ThBC₂ can be defined by

(3)$\Delta E_f=E{(\mathrm{ThBC}}_2)-E\left(\mathrm{Th}\right)-E\left(\mathrm{B}\right)-2\times E\left(\mathrm{C}\right),$

where $E{(\mathrm{ThBC}}_2)$ represents the total energy of ThBC₂, and $E\left(\mathrm{Th}\right)$, $E\left(\mathrm{B}\right)$, and $E\left(\mathrm{C}\right)$ denote the total energies of bulk thorium, α-boron, and graphite, respectively. The calculated formation energy of ThBC₂ is −1.454 eV/f.u., suggesting the thermodynamic stability and possible synthesis of ThBC₂ in experiment. In addition, the electronic band structure calculation reveals that ThBC₂ is metallic, as shown in Figure 2b,c. The density of states (DOS) for ThBC₂ implies that the Th and C atoms mainly contribute to the total density of states around the Fermi level.

To further confirm the stability and assess the mechanical properties, we have calculated the elastic constants and elastic moduli of ThBC₂. The orthorhombic ThBC₂ has nine independent elastic constants, which are C₁₁ = 315 GPa, C₁₂ = 146 GPa, C₁₃ = 40 GPa, C₂₃ = 67 GPa, C₂₂ = 286 GPa, C₃₃ = 576 GPa, C₄₄ = 40 GPa, C₅₅ = 49 GPa, and C₆₆ = 117 GPa. Clearly, these elastic constants satisfy the elastic stability criteria for a stable orthorhombic structure, that is, C₁₁ > 0, C₄₄ > 0, C₅₅ > 0, C₆₆ > 0, C₁₁C₂₂ > C₁₂², C₁₁C₂₂C₃₃ + 2C₁₂C₁₃C₂₃ − C₁₁C₂₃² − C₂₂C₁₃² − C₃₃C₁₂² > 0 [35]. The estimated bulk modulus (B) and shear modulus (G) of ThBC₂ are 184 and 86 GPa, which are larger than the bulk modulus (123 GPa) and shear modulus (61 GPa) of tetragonal ThBC [7]. That is to say, our predicted orthorhombic ThBC₂ has better ability to resist the tensile and shear deformations than that of ThBC. Accordingly, we can obtain the Pugh’s modulus ratio B/G, which is an indicator for determining whether material is ductile or brittle; if Pugh’s ratio is larger than 1.75, a material will tend to be ductile, otherwise brittle [36,37]. The calculated Pugh’s ratio of 2.13 implies that our predicted ThBC₂ tends to be ductile.

The elastic anisotropy refers to the directional dependence of mechanical properties, which is closely related with plastic deformation and crack behavior of materials. To get a comprehensive understanding of elastic anisotropy, we calculated the Young’s modulus E as a function of the directions based on their analytical expressions discussed above, in which the elastic compliance constants can be deduced by the elastic constants. As shown in Figure 3a,b, the Young’s moduli of ThBC₂ along different orientations present a high degree of anisotropy. In particular, we can clearly see that the maximum Young’s modulus is along the [001] direction, which is significantly larger than the Young’s modulus along other directions. This can be ascribed to the fact that there exist a number of strong covalent B–C bonds along the [001] direction, which facilitates ThBC₂ to resist tensile deformation and results in largest Young’s modulus along this direction.

To further study the mechanical anisotropy, we estimated the orientation dependence of Young’s modulus and shear modulus as the tensile or shear strain is within the (001), (100), and (010) planes. For instance, when the tensile strain is within the (001) plane, the tensile direction can be denoted as [hk0], and we refer to the angle between the [hk0] and [100] direction as θ. As a result, the direction cosines of the [hk0] can be obtained as $\alpha =\cos \theta$, $\beta =\sin \theta$, and $\gamma =0$. Thus, Equation (1) translates into

(4)$E^{-1}=s_{11}{\cos }^4\theta +s_{22}{\sin }^4\theta +\left(2s_{12}+s_{66}\right){\cos }^2\theta si{n}^2\theta ,$

in which the Young’s modulus is dependent with the defined orientation θ. The other orientation dependence of Young’s modulus and shear modulus can be obtained using similar method, and the detailed results can be found in Ref. [32]. Based on the above approach, the calculated results are shown in Figure 3c,d. The maximum Young’s modulus is found along the [001] direction (E_[001] = 560 GPa), which shows excellent consistency with the three-dimensional (3D) plots of Young’s modulus. It is found that Young’s modulus along principal crystal directions possess the sequence of E_[101] < E_[011] < E_[010] < E_[100] < E_[110] < E_[001]. To measure the ability of resisting shear deformation for ThBC₂, the dependences of shear moduli on the directions are illustrated in Figure 3d. We can see that the largest shear modulus (G = 117 GPa) appears along the [001] direction on the (100) shear plane. Meanwhile, the obtained smallest shear modulus is along the [100] direction on the (001) shear plane, i.e., G_(001)[100] = 40 GPa, which is equal to the elastic constant of C₄₄.

The ideal strength of a material is the upper limit of mechanical strength that it can sustain. If the stress upon material is larger than ideal strength, it will collapse. In particular, the ideal tensile and shear strengths are intimately related with the failure mode of material, and thus can be employed to uncover when the cleavage fracture and slip will appear for certain compounds. Therefore, we next study the ideal tensile and shear strengths of ThBC₂ to further explore its mechanical property. The tensile and shear strains are continuously applied upon ThBC₂ to obtain the strain–stress relationships, and consequently the ideal strength can be deduced. The calculated strain–stress relationships are shown in Figure 4. It is found that the largest ideal tensile strength is σ = 59 GPa along the [010] direction. The main reason is that the armchair-like B–C chains are along [010] direction, which can form strong covalent B–C against the tensile deformation. The lowest ideal tensile strength appears in the [111] direction, i.e., σ_[001] = 7 GPa, suggesting that the cleavage fracture can easily occur along this direction. An inspection of the Th-C bond length suggests that the distance between the Th4 and C6 increases from 2.66 Å to 3.01 Å when the [111] tensile strain increases from 0 to 0.07. This fact implies these Th-C bonds are the main load bearing component along the [111] direction. Therefore, the weak Th-C bonds along the [111] direction is responsible for this low tensile strength. Besides the [111] direction, the [011] and [110] tensile deformation also easily stretch the weak Th-C bonds and result in small ideal tensile strength along these directions. It is worth noting that the ideal shear strengths of ThBC₂ on different shear planes are all very small, which is the inevitable consequence of the layered structure of ThBC₂. We can see that the highest ideal shear strength is along the [100] direction on the (001) plane, i.e., τ_(001)[100] = 18 GPa, and the obtained lowest ideal shear strength occurs within the [$\overline{1}10$] direction on the (111) plane, indicating that ThBC₂ can readily slip under shear along this direction.

Interestingly, the strain–stress curves under tensile and shear deformation exhibit different characteristic along various directions. For instance, the stress along the [010] direction drops abruptly under tensile deformation, while along the [111] direction the stress decreases gradually rather than a sharp decrease past the elastic limit. To get deeper insights into this interesting characteristic, we next analyze atomic bondings and electron localization function (ELF) to clarify this point. Figure 5a shows the bond lengths of our chosen typical Th-C, B–C, C–C, Th-B, and B–B bonds as a function of the [010] tensile strain. Clearly, we can see that the distance between the B3 and C8 atoms increases abruptly when the tensile strain exceeds the critical value of 0.17, while other bonds show small changes compared with the B–C bonds. Therefore, the alignment of strong covalent B–C bonds along the [010] direction is the main load bearing component, which leads to the largest ideal tensile strength 59 GPa. Nevertheless, once the tensile strain is larger than the critical value of 0.17, as show in Figure 4, a sudden ‘‘hard’’ breaking of the strong covalent B–C bonds occurs, which accounts for the sharp decrease of the stress along this direction. To further confirm this bond-breaking pattern, the three-dimensional isosurfaces of ELF with an isovalue of 0.75 before and after the critical tensile strains are plotted in Figure 5b,c. Generally, the large ELF value corresponds to high electron localization and is the characterization of covalent bondings. As shown in Figure 5b,c, the charge distribution localized in between B3–C8 bond exhibits a sharp decrease when the tensile strain changes from 0.15 to 0.18, further suggesting the bond-breaking of strong covalent B–C bond under the [010] tensile deformation. However, for the [111] tensile deformation, as shown in Figure 5d, the bond lengths of B–C, C–C, Th-B, and B–B bonds nearly keep invariant with the increase of tensile strain, while the Th-C bond increases gradually as the tensile loading increases, even though beyond the critical strain 0.06. Therefore, the stress under this tensile mode change gradually beyond the critical strain, corresponding to the softer bond-deformation behavior. In a word, the ‘hard’’ breaking of B–C bonds correspond to the sharp decrease of the stress in the strain–stress curve, while the weak Th-C bonds result in softer bond-deformation behavior as characterized by a gradual decrease of the stress.

4. Conclusions

In summary, we predicted the crystal structure of ThBC₂ and calculated the mechanical properties of ThBC₂. Our predicted ThBC₂ consists of the Th atomic layers and the armchair-like B–C chains in parallel with Th atomic layers. Based on the phonon spectra and elastic constants, we verified the stability of ThBC₂. We also calculated the anisotropy of Young’s modulus and shear modulus for ThBC₂. The obtained results show that ThBC₂ exhibits significant anisotropy. The calculated strain–stress curves suggest ThBC₂ has low ideal strengths, implying that the cleavage fracture and slip can easily occur in ThBC₂. We expect that our findings can favor the identification of structure for ThBC₂, and provide insight understanding into the mechanical properties of ThBC₂.

Supplementary Materials

The following are available online at https://www.mdpi.com/2073-4352/9/8/389/s1, Figure S1: Crystal structure and phonon spectrum of Pmmm ThBC₂, Figure S2: Crystal structure and phonon spectrum of C2/m ThBC₂, Figure S3: Crystal structure and phonon spectrum of $\mathrm{I}\overline{4}m2$ ThBC₂, Figure S4: Crystal structure and phonon spectrum of Amm2 ThBC₂, Table S1: Lattice parameters, atomic positions and the total energies per atom of four metastable ThBC₂, Table S2: Calculated elastic constants C_ij of the most stable ThBC₂ without and with SOC, most stable structure.cif (The most stable structure of ThBC₂), no.47.cif (Crystal structure of Pmmm ThBC₂), no.12.cif (Crystal structure of C2/m ThBC₂), no.119.cif (Crystal structure of $\mathrm{I}\overline{4}m2$ ThBC₂), no.38.cif (Crystal structure of Amm2 ThBC₂).

Author Contributions

Formal analysis, X.Z. and B.Z.; Investigation, X.Z.; Writing—original draft, X.Z. and B.Z.

Funding

This work is supported by the Natural Science Foundation of China (Grant No. 11647007), the Natural Science Basic Research plan in Shaanxi Province of China (Grant Nos. 2018JQ1083), the Baoji University of Arts and Sciences Key Research (Grant No. ZK2017009), and the National Undergraduate Innovative Training Program (Grant No. 201610721017).

Conflicts of Interest

The authors declare no conflict of interest.

Figures

Figure 1. Side (a) and top (b) views of the crystal structure of ThBC2, in which the purple (large), red (middle), and blue (small) spheres denote Th, C, and B atoms, respectively. (c) The infinite armchair-like B–C chain along y-axis direction. (d) The square pyramid building block consisting of one Th and four C atoms.

Figure 2. (a) Phonon spectrum of ThBC2 along the high-symmetry lines. (b) Band structure of ThBC2 along the high-symmetry lines. (c) DOS of ThBC2. Here the vertical blue dashed lines represent the positions of high-symmetry point in the Brillouin zone, and the horizontal blue dashed line denotes the Fermi level.

Figure 3. 3D plots of Young’s modulus (a) and its plane projections (b) as the function of the directions, in which the distances between the origin and the surface are the values of Young’s moduli along corresponding directions. Orientation dependence of Young’s modulus (c) and shear modulus (d), in which the principal crystal orientations [110], [101], and [011] are marked by arrows with different color.

Figure 4. The stress as the function of tensile (a) and shear (b) strains along different tensile and shear strain directions.

Figure 5. Bond lengths as a function of strain under the [010] (a) and [111] (d) tensile modes. The ELF for ThBC2 before (b) and after (c) the critical [010] tensile strain.