Structural, mechanical, and thermodynamic properties of newly-designed superhard carbon materials in different crystal structures: A first-principles calculation

Computational Materials Science 171 (2020) 109229

Contents lists available at ScienceDirect

Computational Materials Science journal homepage: www.elsevier.com/locate/commatsci

Structural, mechanical, and thermodynamic properties of newly-designed superhard carbon materials in different crystal structures: A first-principles calculation

T

Wen-Duo Han, Ke Li , Jia Dai, Yao-Hui Li, Xun-Du Li ⁎

School of Mechanical and Electrical Engineering, Nanchang University, China

ARTICLE INFO

ABSTRACT

Keywords: Structural Mechanical Thermodynamic Carbon First principles

At present work, the ab-initio method was applied to systematically and deeply investigate the structural and mechanical properties, such as elastic, elastic anisotropy, thermal conductivity and thermodynamic properties of newly-designed superhard carbon in different crystal structures. The results of cohesive energies and formation enthalpies indicate that all carbon are energetically and structurally stable and P2221 carbon has the most stable structure among these carbon allotropes. On the basis of the single-crystal elastic constants, the elastic properties, such as bulk modulus B, shear modulus G and Young's modulus E, Poisson's ratio v and hardness HV were obtained, the values of G, E and HV is largest for P2221 carbon. The order of elastic anisotropy is C2/m28 > Imma > C2/m-16 > P21/m > C2/m-20 > Cm-32 > Cm-40 > Amm2 > P2221 > I-4. Clarke's and Cahill's models are employed to calculate the minimum thermal conductivities of these carbon, and the order is P2221 > P21/m > C2/m-20 > C2/m-16 > Cm-40 > Imma > Amm2 > Cm-32 > I-4 > C2/m-28. Furthermore, the minimum thermal conductivities of these carbon allotropes are anisotropic in different direction. The analysis of thermodynamic properties can be obtained that the thermal stability of these newlydesigned carbon allotropes is better than that of diamond and the order is Cm-40 > C2/m-28 > I4 > Imma > Amm2 > Cm-32 = C2/m-20 > P2221 > C2/m-16 > P21/m > Diamond, Cm-40 carbon has the strongest thermal stability. It can be concluded from phonon density that the C2/m-28 carbon is dynamically metastable structure. The study in this work give a comprehensive understanding of structure, elastic, elastic anisotropy, thermal conductivity and thermodynamic properties for the newly-designed superhard carbon, which can be an important guidance for the theoretical work and the fabrication in experiment in the future.

1. Introduction Superhard materials play a significant role in the application of modern industry on account of their superior mechanical properties, such as cutting, drilling, polishing, and surface protecting coating [1,2]. To satisfy the requirement of industrial manufacture, a great many investigations have been focused on developing and fabricating advanced superhard materials. Carbon has a variety of allotropes, such as graphite, graphene, C60, diamond, hexagonal diamond (lonsdaleite), and carbon nanotubes due to its unique bond character [3]. Meanwhile, carbon is also one of the most plentiful and cheapest element in the nature, which has been widely applied to synthesize different structural superhard materials by cold-compressing graphite process [4–8]. Hexagonal-diamond, M-carbon, bct4-carbon, Z-carbon, W-carbon, F-carbon, T12-carbon, chiral-framework carbon, and BC8-carbon have

⁎

been theoretically and experimentally confirmed the superhard carbon structures up to now. F. Tian, and X. Dong et al. [9] reported Hexagonal-diamond and predicted Hv close to 96 GPa. Oganov and Glass firstly theoretically proposed the monoclinic M-carbon in previous research [10], which has been verified to be the superior superhard carbon by Panzik and Wang et al. [11,12] later. Other superhard structures carbon mentioned above possess a Vickers hardness of over 80 GPa. Recently, Zhang et al. [13] designed various superhard carbon phases under the special conditions based on the CALYPSO algorithm, and the predicted structures can be synthesized by cold-compressing graphite [4–8], evolutionary metadynamics technique [14–17] due to their energies are close to that of M-carbon, which has been confirmed and fabricated by experiments at high pressure [11,12]. Carbon also posseses various hybridization form, You et al. [18] investigated the a novel topological carbon with sp, sp2, and sp3 hybridized bonds in one

Corresponding author at: School of Mechanical and Electrical Engineering, Nanchang University, 999 Xuefu Road, Honggutan District, Nanchang 330031, China. E-mail address: [email protected] (K. Li).

https://doi.org/10.1016/j.commatsci.2019.109229 Received 27 June 2019; Received in revised form 13 August 2019; Accepted 26 August 2019 0927-0256/ © 2019 Elsevier B.V. All rights reserved.

Computational Materials Science 171 (2020) 109229

W.-D. Han, et al.

Table 1 Experimental and theoretic lattice parameters, cohesive energies Ecoh and formation enthalpies ΔH of carbon. Space group

Lattice parameters (Å)

V (Å3)

Ecoh(eV/atom)

ΔH(eV/atom)

Refs.

P2221

3.559 5.544 3.509 3.555 8.537 3.506 3.549 8.521 3.499 7.600 4.287 4.366 7.596 4.283 4.363 2.522 5.556 8.642 2.520 5.550 8.636 5.576 5.576 5.521 5.571 5.571 5.516 4.728 5.057 7.165 4.726 5.053 7.161 4.721 5.046 7.154 5.810 7.574 9.158 5.805 7.570 9.150 6.056 2.533 4.025 6.041 2.527 4.018 14.263 2.527 6.617 14.265 2.524 6.608 9.030 2.525 8.723 9.019 2.524 8.716 12.652 2.559 5.741 12.246 2.532 5.741 3.577 3.577 3.577 3.567 3.567 3.567 3.566 3.566 3.566

106.696 106.404 105.813 142.247 142.912 121.079 120.780 171.657 171.192 95.612 95.424

−10.017

0.1028

Present [1319]

−9.987

0.1350

−9.960

0.1600

−9.984

0.1357

−9.941

0.1514

Present [13] Present [13] Present [13] Present [1319]

−9.976

0.1438

−9.935

0.1850

−9.958

0.1325

−9.995

0.1239

−9.374

0.7460

−10.145

−0.025

Imma Amm2 I-4 C2/m-16 Cm-32 P21/m Cm-40 C2/m-20 C2/m-28 Diamond

195.814 195.264 59.779 59.400 238.513 237.920 118.502 118.200 175.377 164.808 45.756 47.704 45.347

structure. Moreover, these predicted carbon structures are parallel to that of diamond [13]. Furthermore, the orthorhombic Imma and P2221 structures has lower transition energy and pressure than those in the Mcarbon at zero pressure [13], further indicating that these two structures are theoretically available. The 5 + 7 or 5 + 6 + 7 atom rings are contains in a large proportion of these designed monoclinic structure, such as the C2/m-16, P21/m, Cm-32, C2/m-20, Cm-40 and C2/m-28 structures. C2/m structure has been deeply investigated by Xing and Li et al. [19] in terms of structural, mechanical, and electronic properties under the pressure ranging from 0 to 100 GPa, and C2/m carbon was confirmed structurally and dynamically stable from 0 to 100 GPa, besides, they also observed the modulus of C2/m structure close to that of diamond, indicating the C2/m structure has a comparable mechanical properties with diamond. However, there are still insufficient systematic investigations on these newly-designed superhard carbon of the physical properties for the limitation of experiments, especially anisotropy of elastic modulus, which is significant in evaluating the mechanical durability. Besides, it is also worthwhile and essential to systematically investigate anisotropy of elastic, thermal conductivity and thermodynamic properties due to different degree of anisotropy of elastic can lead to microcracks in materials. On the meantime, the thermal conductivity and thermodynamic properties is also an important parameter to measure the performance of superhard materials under the high temperature conditions. All in all, a theoretical investigation benefits for breaking through the limitations of experiments. First-principle calculation, as a powerful and valid method, can systematically detail the mechanical and thermodynamic properties, solving the difficulties in experiments. Thus, we employ the first-principle calculation to comprehensively analyzed the structures and physic properties of carbon in C2/m, Cm, P21/m, I-4, Amm2 P2221 and Imma carbon. The P2221, Imma, Amm2 and I-4 carbon is the structure of orthorhombic and C2/m-16, Cm-32, P21/m, Cm-40, C2/m-20, C2/m-28 is monoclinic structure. The hardness and lattice parameters in our present work are consistent well with the previous literature, indicating the work present is feasible and reliable. The work mainly focuses on the elastic, elastic anisotropy, thermal conductivity and thermodynamic properties, aiming at providing a meaningful guidance to the experimental synthesis and application of newly-designed superhard carbon allotropes.

Present [13] Present [13] Present [13] Present [13] Present [13] Present [2728]

2. Computational methods GGA of Perdew-Burke-Ernzerhof PBE is employed for all calculations [20]. The elastic tensor is further analyzed after the all carbon allotropes are fully relaxed to the equilibrium state. The stress-strain [21] relationships with PBE are employed to describe the elastic constants. The calculations of the total energy, elastic, elastic anisotropy properties, thermal conductivity and thermodynamic property are acquired by Cambridge sequential total energy package (CASTEP) based on density function theory (DFT) [22]. The interaction between ioncore and valence electrons was described by Ultrasoft pseudo-potentials [23]. The valence electron configurations of carbon was 2s22p2. The kpoint are 2 × 1 × 5, 2 × 3 × 3, 5 × 3 × 4, 3 × 6 × 3, 1 × 6 × 3, 5 × 2 × 3, 1 × 6 × 2, 2 × 6 × 4, 6 × 3 × 2, 3 × 3 × 3, 4 × 4 × 4 for P2221, Imma, C2/m-16, C/2m-20, C2/m-28, Cm-32, Cm-40, P21/m, Amm2, I-4 and diamond, respectively. A finite basis set correction and the Pulay scheme of density mixing was employed to describe the energy and stress [24,25]. The geometry optimization is performed based on Broyden-Fletcher-Goldfarb-Shannon (BFGS) method [26]. The total energy convergence was converged to 1.0 × 10−5 eV/atom, the force on each atom was converged to 0.03 eV/Å, the maximum strain amplitude was set as 0.05 GPa and the maximum displacement was converged to 0.0001 Å, the cut-off energy was set as 300 eV for plane wave expansions. 3. Result and discussion 3.1. Structural properties and phase stability The mechanical and physical properties, such as energy and elastic constants are calculated after the lattice constants and internal coordinates were optimized to equilibrium state. Table 1 lists the equilibrium parameters after optimization. Fig. 1 shows the structures of these novel carbon allotropes. Meanwhile, the dates of these carbon allotropes at previous literature [13] are also shown in Table 1 to compare with our present work. It is obvious from Table 1 that symmetries of these carbon alltropes unchanged after these structures are optimized on account of the unvaried space group. What’s more, the lattice parameters in this work is slightly overestimated due to the GGA 2

Computational Materials Science 171 (2020) 109229

W.-D. Han, et al.

3

(caption on next page)

Computational Materials Science 171 (2020) 109229

W.-D. Han, et al.

Fig. 1. Unit cell crystal structures of the novel carbon allotropes.

approximation applied in present wok. However, the structural parameters of carbon allotropes at present work is well consistent with the values at previous literature within an acceptable deviation of 1%, implying that the computational method we adopted and the setting in present calculation are feasible and credible. To a great extent, cohesive energy can be criterion to evaluate the structural stability of crystal, which is defined as the work that is needed when the crystal decomposes into single atom. Generally, a negative value of the cohesive energy corresponds to a stable structure. In addition, the value of formation enthalpy can be employed to predict the alloying ability. Generally, a negative formation enthalpy manifests crystal can be fabricated and exist stably, while the positive value of the formation enthalpy indicates an endothermic reaction, implying that the structure is not easy to be synthesized. As is known to us that the stability of structure play a significant role at the processing of designing and synthesizing the superhard material. Hence, the cohesive energies Ecoh and formation enthalpies ΔH should be paid attention to better design and fabricate excellent superhard material. Generally, lower values of ΔH and Ecoh correspond to a better structural stability. The calculated formation enthalpies ΔH and cohesive energies Ecoh of carbon allotropes have been listed in Table 1. It is obvious that the values of formation enthalpies and cohesive energies are slightly various. Meanwhile, the value of ΔH and Ecoh in this work would be a significant guidance for theoretical and experimental research due to lacking of relative literature for reference yet. Fig. 2 clearly reflect the the correlation between formation enthalpies and cohesive energies. It is obvious from Fig. 2 that a lower value of cohesive energies corresponds to a lower formation enthalpies. For orthorhombic structure, it can be seen that P2221 carbon has the lowest cohesive energies (−10.017 eV/atom) and formation enthalpy (0.1028 eV/atom), indicating P2221 carbon is the most energetically and structurally stable and easier to form under the same condition. While Amm2 carbon has the largest cohesive energies (−9.960 eV/ atom) and formation enthalpy (0.1600 eV/atom), implying the stability of Amm2 carbon is the weakest and it is difficult to form. As for monoclinic structure, the C2/m-28 carbon has the largest cohesive energies and formation enthalpy than that of C2/m-16, C/2m-20, Cm32, Cm-40, and P21/m, indicating that C2/m-28 has relatively weaker structural stability and is unlikely to form. It can also deduced from the analysis above that P2221 carbon is easier to be synthesized comparing with other carbon allotropes under the same condition.

3.2. The elastic constants of single crystal The elastic constants of single crystal, as a important parameter in evaluating the performance of the application in tribology and surface engineering, is essential to be further investigated for these superhard carbon structures. In this work, we adopted the the stress-strain approach to calculate the single-crystal elastic constants. The elastic constants of single crystal Cij, together with the relevant theoretical values available in previous work, is listed in Table 2 [13,19,28,29]. It is obvious that the elastic constants are consistent well with the values in previous literature [13,19,28,29], indicating the present work is feasible and reliable. On the basis of the lattice dynamical theory by Born and Huang [30], we obtained the mechanical stabilities of C2/m, Cm, P21/m, Amm2, I-4, P2221 and Imma carbon in orthorhombic and monoclinic structures. The data in Table 2 indicates that these carbon are structurally stable. For orthorhombic carbon, such as Amm2, I-4, P2221 and Imma, there is a relationship of C11≠ C22 ≠ C33, implying that these elastics constants are anisotropic, leading to the elastic modulus being dependent upon directions. Fig. 3 shows the elastic constants C11, C22, C33, C44 and C66. The resistances of x-axis, y-axis and z-axis to linear stress are correlated with the single-crystal elastic constants C11, C22, and C33, respectively. Generally, larger values of C11, C22, and C33 correspond to stronger resistance to compression of x-axis, y-axis and z-axis, respectively. For P2221, Amm2, C2/m, Cm and P21/m carbon except C2/m-20 structure, it is obvious the value of C22 is largest comparing to that of C11 and C33, indicating that these carbon structures are most incompressible under the y-axis uniaxial compression. When it to Imma and I-4 carbon, the incompressibility of z-axis is strongest due to their largest values of C33. What’s more, the values of C44 and C66, are related to the tolerance to shear deformation at (1 0 0) plane along the directions of [0 0 1] and [1 1 0], respectively. As is shown in Fig. 3, P2221 carbon possesses the stronger resistance to shear deformation at (1 0 0) plane along the directions of [0 0 1] due to its largest value of C44, while the tolerance to shear deformation at (1 0 0) plane along the [1 1 0] direction is the strongest for Cm-40 carbon. Meanwhile, it can be observed that the values of C44 and C66 are different for these carbon allotropes, implying that these structure is anisotropic for shear modulus. Furthermore, the deviation between C11 and C33 is not equal to that between C44 and C66 for P2221, Imma, I-4, C2/m-16, Cm-32 carbons, indicating the bulk modulus and shear modulus has different anisotropy for these structures. 3.3. Polycrystalline elastic modulus Elastic modulus, as a parameter evaluating mechanical properties of polycrystal, is also significant in theory and applications. On the basis of Voigt, Reuss and Hill approximations [33–35], we calculated the mechanical properties of polycrystal, such as Young's modulus E, bulk modulus B, shear modulus G, and Poisson's ratio v. Furthermore, theoretical polycrystalline elastic modulus are deduced in accordance with the average value of two bounds including of Voigt and Reuss bound. The formula can be listed as follows:

BH =

1 1 (BV + BR) GH = (GV + GR) 2 2

(1)

Besides, the Young’s modulus (E) and poisson’s ratio (v) can also be deduced from the bulk modulus (B) and shear modulus (G) on the basis of the following formula [36]:

E=

Fig. 2. The calculated cohesive energies Ecoh and formation enthalpies ΔH of carbon in our work (a);

9BHGH 3BH + GH

=

3BH EH 6BH

(2)

Meanwhile, the calculated formulas of bulk modulus B, shear 4

Computational Materials Science 171 (2020) 109229

W.-D. Han, et al.

Table 2 Calculated elastic constants Cij (in GPa) of carbon. Space group

C11

C12

C13

C15

C25

C35

C22

C23

C33

C44

C55

C66

C46

Refs.

P2221

927.3 949.0 1047.2 1008.0 1126.0 1023.4 900.6 984.5 939.6 970.5 907.8 1118.4 768.4 1059.8 1055.0 1053.0 1076.0

75.5 112.0 47.1 61.0 86.0 6.4 87.2 60.3 53.6 3.5 69.2 28.7 −50.2 102.4 119.0 120.0 125.0

119.6 141.0 82.6 79.0 121.0 56.7 77.6 78.2 96.4 80.4 115.9 91.3 103.9 –

–

–

–

–

–

– – – 2.7 7.2 14.2 27.8 17.5 −60.2 –

922.9 949 0.0 1131.3 1080.0 1196.0 877.2 962.8 1020.1 964.1 1001.6 991.5 890.1 823.1 –

525.8 524.0 427.8 425.0 446.0 401.8 405.1 436.7 388.4 509.8 441.3 467.9 280.2 561.0 569.0 563.0 577.0

510.4 516.0 506.7 488.0 526.0 434.4 405.1 488.1 482.3 443.2 462.4 452.5 391.4 –

460.4 435.0 377.3 355.0 377.0 441.6 465.1 463.1 386.8 385.0 481.6 472.7 210.3 –

Present [19] Present [1913]

– – 54.2 −24.2 −29.5 −33.2 19.8 12.0 –

88.2 103.0 19.4 37.0 61.0 56.2 77.6 65.5 43.2 59.5 44.3 73.8 72.5 –

– -

– – 8.9 1.2 58.7 −22.0 0.5 −74.9 –

1020.4 1071.0 1024.4 992 0.0 1105.0 1111.2 910.1 1027.7 1029.6 1007.4 1079.2 1065.5 836.6 –

Imma Amm2 I-4 C2/m-16 Cm-32 P21/m Cm-40 C2/m-20 C2/m-28 Diamond

– – 73.8 −16.5 −11.16 −50.16 54.88 35.6 –

Present [192829]

modulus G in VRH theory are correlated to the orthorhombic and monoclinic symmetries. The calculated elastic modulus and Poisson's ratios of carbon are listed in Table 3 and shown in Fig. 4. The bulk modulus B can be employed to measure the performance of tolerance to volume change by external pressure [37]. Generally, a higher bulk modulus corresponds to stronger resistance to volume change under the applied stress. It can be observed that the bulk modulus of these superhard carbon ranging from 329.3 GPa to 383.4 GPa, indicating different structures of carbon correspond to different compressibility. It is obvious that C2/m-20 structure has highest bulk modulus close to that of diamond, which implies C2/m-20 structure has the strongest resistance to volume change by external pressure among these newlydesigned carbon. And C2/m-28 structure exhibits the weakest incompressible characteristic among these carbon structures due to its smallest bulk modulus (329.3 GPa). The shear modulus G indicates the tolerance to shape change upon shear stress. A larger value of G implies stronger resistance to external shear stress. It can be observed that P2221 has the largest shear modulus (487.7 GPa), indicating that P2221 has the strongest tolerance to shear process. Furthermore, it can be found that the least formation enthalpy carbon structure (P2221)

Fig. 3. The calculated elastic constants for carbon allotropes.

Table 3 Calculated bulk modulus B (in GPa), shear modulus G (in GPa), Young's modulus E (in GPa), Poisson's ratio v, Pugh's modulus ratio GH/BH, and Vickers hardness HV (in GPa) of carbon. Space group

B

G

E

ν

B/G

HV

ΔH(-1 × 10−3)

Refs.

P2221

383.5

487.7

1026.9

0.053

0.786

282

Imma

377.1

447.3

967.1

0.073

0.840

Amm2

373.6

440.2

948.2

0.077

0.849

I-4

360.9

422.0

910.9

0.080

0.855

C2/m-16

381.7

462.7

988.6

0.068

0.825

Cm-32

368.8

431.8

931.8

0.079

0.854

P21/m

382.3

474.5

1006.9

0.061

0.805

Cm-40

381.2

457.7

980.6

0.071

0.832

C2/m-20

383.4

467.8

997.6

0.066

0.820

C2/m-28

329.3

388.4

836.4

0.077

0.848

Diamond

421.5

526.5

1115.2

0.059

0.800

92.8 85.9 83.7 87.8 82.3 84.0 79.5 83.0 87.8 84.9 80.7 83.1 91.7 85.4 86.2 87.2 89.1 85.9 78.9 86.1 98.4 89.5 97.5

Present [13] Present [13] Present [13] Present [13] Present [13] Present [13] Present [13] Present [13] Present [13] Present [13] Present [1331]

5

256 216 252 206 227 201 213 264 187 291

Computational Materials Science 171 (2020) 109229

W.-D. Han, et al.

Fig. 4. Calculated elastic modulus B, G and E (a) Pugh's ratio B/G, and Poisson's ratio (b) for carbon allotropes.

corresponds to a largest value of shear modulus G, while the largest formation enthalpy carbon structure (C2/m-28) has the least value of shear modulus G, and the order of shear modulus G is P2221 > P21/ m > C/2m-20 > C2/m-16 > Cm-40 > Imma > Amm2 > Cm32 > I-4 > C2/m-28. Young’s modulus, defined as the ratio of stress to strain, is closely correlated to stiffness of materials. A larger Young’s modulus indicates a stiffer material. It can be observed from Table 3 that the order of stiffness of these carbon is P2221 > P21/m > C/2m20 > C2/m-16 > Cm-40 > Imma > Amm2 > Cm-32 > I4 > C2/m-28, and P2221 has the highest stiffness among these structures due to its largest value of Young’s modulus. The poisson’s ratio v, which is usually employed to evaluate the stability of crystal against to shear deformation, ranging from −1 to 0.5. Generally, the more the value of poisson’s ratio is, the better the plasticity. It can be seen from Table 3 that the value of poisson’s ratio varies from 0.053 to 0.079. Moreover, the I-4 carbon possesses the best plasticity while P2221 has the worst plasticity. Theoretically, poisson’s ratio, BH/GH [38–40] can be employed to evaluate the intrinsic ductility and brittleness performance of non-metallic materials carbon. Generally, the material with v > 0.26, BH/ GH > 1.75 is ductile. Otherwise, the material is brittle. The data from Table 3 shows that all these structure carbon have the values of v smaller than 0.26, BH/GH lower than 1.75, implying all these carbon exhibit brittle. Moreover, P2221 carbon exhibit the most brittle behavior due to its lowest value of v and BH/GH. The largest values BH/GH (0.080) and v (0.855) for I-4 carbon also indicates the relatively ductile characteristic among these carbon allotropes. On the basis of the macroscopic concept bulk and shear modulus, the Vickers hardness of polycrystalline materials can be deduced by the following formula [41]:

HV = 2(G3/B2)0.585

3

(3)

Theoretically, the macroscopic parameters is essential and significant to be further investigated on the account of its close correlation with the hardness of materials, and it also provide a powerful guidance to design and synthesize advanced superhard materials. The calculated results and the data of other literature at previous work are listed in Table 3 based on the theoretical model. It is known that the value of Vickers hardness with more than 40 GPa can be as superhard material. As is shown in Table 3, it can be observed that all these carbon has high value of Vickers hardness with more than 70 GPa, close to that of diamond, indicating these carbon can be excellent candidate for a superhard material. Moreover, the value the P2221 structure has largest value of Vickers hardness (92.8 GPa) except diamond. Generally, the hardness is correlated with shear modulus and bulk modulus. Therefore, we plot the correlation between the bulk modulus (shear modulus) and hardness in Fig. 5. Moreover, the bulk modulus can be a criterion to evaluate the hardness of material in theoretical work and practical application on account of the relationship between the bulk modulus and hardness that a larger bulk modulus indicates higher hardness [2,42]. However, it can be found from Fig. 5(a) that larger bulk modulus not entirely corresponds to higher hardness, which also have been confirmed by previous literature [43,44]. Shear modulus, as parameter closely related to dislocation formation and movement [44,45], is more suitable to predict the hardness. It is obvious from Fig. 5(b) that higher values of shear modulus indicate a larger hardness. Hence, C2/m-28 carbon has the smallest value of shear modulus (836.4 GPa) accordingly, it can be observed the value of C2/ m-28 carbon is the lowest (78.9 GPa). Similarly, P2221 carbon possesses the largest value of shear modulus (1026.9 GPa), implying the highest value (92.8 GPa) of hardness among these carbon allotropes. Indeed, 6

Computational Materials Science 171 (2020) 109229

W.-D. Han, et al.

Fig. 5. Correlations between bulk modulus B (shear modulus G) and hardness HV of carbon.

the more the value of shear modulus is, the larger the hardness.

deviates from 1implies the higher elastic anisotropy of bulk materials. Table 4 lists the results of anisotropic indexes (AU). Obviously, C2/ m-28 carbon has the highest elastic anisotropy for the largest values of anisotropic indexes (AU = 0.262). while I-4 (AU = 0.029) carbon possesses the most elastic isotropy among these structures. Besides, the value of Acomp and Ashear can also be employed to characterize the degree of elastic anisotropy. Obviously, C2/m-28 carbon possess the largest values of Acomp (32.8%) and Ashear (29.6%), implying the highest elastic anisotropy of C2/m-28 carbon, while the value of Acomp mostly close to zero for the other carbon allotropes, indicating bulk elastic isotropic characteristic of them. I-4 carbon has the least shear elastic anisotropy for its smallest value of Ashear (0.3%). Moreover, the degree of shear anisotropy in the (1 0 0), (0 1 0) and (0 0 1) plane can be reflected by the absolute value of (A1-1), (A2-1) and (A3-1), respectively. The values of (A1-1) and (A2-1) for P2221 carbon has the most deviation from 1, indicating the highest shear anisotropy of P2221 carbon in (1 0 0) and (0 1 0) planes, while C2/m-28 carbon possesses the largest value of (A3-1), suggesting the most highly shear anisotropy of C2/m-28 carbon in (0 0 1) plane. At present work, several anisotropic indexes are employed to describe the elastic anisotropies, and the order of universal elastic anisotropic index AU is C2/m-28 > Imma > C2/m-16 > P21/ m > C2/m-20 > Cm-32 > Cm-40 > Amm2 > P2221 > I-4 (See Fig. 6). To effectively describe the directional elastic anisotropy of crystal, the 3D surface construction of elastic modulus is employed to characterize the elastic anisotropy of carbon along the corresponding direction. The Young's modulus E and bulk modulus B are calculated based on the formula as follows [49,50]: For orthorhombic crystal class:

3.4. Anisotropy of elastic modulus Taking the correlation between the elastic anisotropy of material and mechanical durability into consideration, the anisotropy of elastic should be deeply investigated on account of which can induce microcracks in materials. In this work, several different approaches are employed to describe the elastic anisotropy such as universal elastic anisotropy index (AU) and the percent anisotropy (Acomp for the anisotropy of the compressibility and Ashear for the shear anisotropy) [46]. Generally, larger values of AU, Acomp and Ashear correspond higher elastic anisotropy. In addition, the shear anisotropic factors (A1 for the (1 0 0) plane, A2 for the (0 1 0) plane and A3 for the (0 0 1) plane) are calculated based on the formula as follows [46–48]:

AU = 5GV / GR + BV / BR Acomp = (BV = (GV

(4)

6

BR/ BV + BR) × 100% Ashear (5)

GR/GV + GR) × 100%

A1 = 4C 44/(C11 + C 33

2C13)

(6)

A2 = 4C55/(C 22 + C 33

2C 23)

(7)

A3 = 4C 66/(C11 + C 22

2C12)

(8)

where B and G correspond to the bulk and shear modulus, respectively, while the subscripts V and R represents the Voigt and Reuss limits, respectively. Cij are the single-crystal elastic constants. A1 = A2 = A3 = 1 indicates the elastic isotropic for bulk materials. While the more value of shear anisotropic factors (A1, A2 and A3)

1/ E = l14 S11 + l24 S22 + l34 S33 + l 22 l 32 (2S23 + S44) + l12 l32 (2S13 + S55) + l12

Table 4 Calculated elastic anisotropic indexes (AU, Acomp, Ashear, A1, A2 and A3) of carbon. Space group

AU

Acomp (%)

Ashear (%)

A1

A2

A3

P2221 Imma Amm2 I-4 C2/m-16 Cm-32 P21/m Cm-40 C2/m-20 C2/m-28 Diamond

0.046 0.083 0.055 0.029 0.077 0.067 0.074 0.062 0.068 0.262 0.030

0.1 0.1 0.3 0.1 0.1 0.1 0.1 0.1 0.3 32.8 0.0

0.4 0.7 0.5 0.3 0.7 0.6 0.7 0.6 0.6 29.6 0.0

1.313 0.850 0.899 0.949 0.945 0.908 1.126 1.059 1.021 0.810 –

1.155 0.957 0.926 0.943 1.019 1.011 0.938 0.933 1.001 1.034 –

1.024 0.764 0.867 1.013 0.979 0.831 0.782 1.041 0.888 0.493 –

l22 (2S12 + S66 )

(9)

1/ B = (S11 + S12 + S13) l12 + (S12 + S22 + S23) l22 + (S13 + S23 + S33) l32 (10) For monoclinic crystal class:

1/ E = l14 S11 + 2l12 l22 S12 + 2l12 l32 S13 + 2l13 l3 S15 + l24 S22 + 2l22 l32 S23 + 2l1 l22 l3 S25 + l34 S33 + 2l1 l33 S35 + l22 l32 S44 + 2l1 l22 S46 + l12 l32 S55 + l12 l22 S66 (11)

1/ B = (S11 + S12 + S13) l12 + (S12 + S22 + S23) l 22 + (S13 + S23 + S33) l32 + (S15 + S25 + S35) l1 l3 7

(12)

Computational Materials Science 171 (2020) 109229

W.-D. Han, et al.

characteristic than C2/m-28 in (0 0 1) plane while has no obvious anisotropy in (0 1 0) and (1 0 0) planes, moreover, (0 1 0) pane of these carbon allotropes has relatively weaker anisotropy than that of (1 0 0) and (0 0 1) planes overall. As for bulk modulus, C2/m-28 carbon exhibit obvious anisotropy in (0 1 0) plane while shows weak anisotropic characteristic in (1 0 0) and (0 0 1) planes, furthermore, (1 0 0) pane of these carbon allotropes exhibit relatively larger anisotropy than (0 1 0) and (0 0 1) planes overall. 3.5. Debye temperatures and anisotropy of sound velocities Debye temperature can be employed to characterize the phenomena of solid-state physics such as specific heat, melting point and thermal expansion. On the basis of the elastic modulus, Debye temperature can be obtained from the sound velocities [51,52]. Debye temperatures and sound velocities are plotted in Table 5 and shown in Fig. 10(a). Theoretically, Debye temperature can be employed to describe the strength of covalent bonds [53]. Generally, higher value of Debye temperature indicates stronger strength of covalent bonds. Obviously, P2221 carbon has the largest value of Debye temperature, implying the chemical bond of P2221 carbon is strongest among these carbon allotropes except that of diamond, while the value of Debye temperature for C2/m-28 carbon is smallest, indicating the weakest the chemical bond of C2/m-28 carbon. The order of Debye temperature is P2221 > P21/m > C2/m16 > C2m-20 > Imma > Cm-40 > Amm2 > Cm-32 > I-4 C2/m28. In addition, Debye temperature has a close correlation with the thermal conductivity, and a more value of Debye temperature indicate a higher thermal conductivity. Thus, the thermal conductivity of P2221 carbon should be largest, while C2/m-28 has the least thermal conductivity, which can be confirmed by the following results. The acoustic Grüneisen constant γɑ can be employed to characterize the anharmonic correlation between atoms in solids. Moreover, some significant parameters, such as the effect of temperature on elastic properties, acoustic wave absorption and thermal conduction, are dominated by acoustic Grüneisen constant γɑ [54]. Grüneisen constants are defined as follows [55]:

Fig. 6. Calculated elastic anisotropic indexes (AU, A1, A2 and A3).

where Sij is the usual elastic compliance constant in Table 4, l1, l2 and l3 are the direction cosines. Generally, a spherical shape of the 3D surface implies the crystal is isotropic, the degree of deviation from 3D shape sphere represents the degree of anisotropy. Bulk modulus and Young's modulus are plotted in Figs. 7 and 8, respectively. Obviously, C2/m-28 carbon exhibits the highest deviation from sphere in terms of bulk modulus, implying the largest anisotropic bulk modulus of C2/m-28 carbon among these carbon allotropes. Although 3D surface graphs of other carbon allotropes are much close to sphere, 3D constructions surface of C2/m-20 and Amm2 carbon show higher deviation than that of the other carbon except C2/m-28 carbon, indicating larger anisotropies of bulk modulus than the other carbon except C2/m-28 carbon. Moreover, the conclusion of the anisotropic bulk modulus based on 3D graphs evidently coincides with that of compressible anisotropy Acomp in Table 4 accordingly, C2/m-28 carbon exhibits the largest anisotropic bulk modulus due to its highest values of Acomp (32.8%), and then is C2/m-20 and Amm2 carbon. As for the 3D surface constructions of Young's modulus, it can be observed that these carbons show more deviation from sphere, indicating larger anisotropy than that of bulk modulus, meanwhile, Fig. 8. shows that C2/m-28 carbon exhibits anisotropic characterization, and the anisotropy of P2221 carbon is slightly than that of C2/m-28 carbon, while I-4 carbon exhibits the least anisotropic nature, consistent with the analysis of universal elastic anisotropic index AU. The elastic anisotropy can be characterized by the 3D surface of elastic modulus. To more intuitively show the elastic anisotropy, the projections of elastic modulus in different planes ((0 0 1), (0 1 0) and (1 0 0) crystal planes) are shown in Fig. 9 based on the 3D surface of elastic modulus. For Young's modulus, it can be found that C2/m-28 carbon exhibits the largest anisotropy nature in (0 0 1), (0 1 0) and (1 0 0) crystal planes, and Imma carbon shows slightly less anisotropic

a=

3 2

4vt2

3vl2

(13)

vl2 + 2vt2

where vl and vt are longitudinal sound velocity and transverse sound velocity, respectively. Table 5 and Fig. 10(b) show the value of γɑ. Generally, a larger value of γɑ indicates more temperature (or pressure) dependence of the crystal volume. It can be observed that C2/m-28 carbon has the largest value of γɑ, indicating the greatest temperature (or pressure) dependence of the crystal volumes for C2/m-28 carbon while the least effect on P21/m due its lowest value of γɑ. It is known to us that the transverse and longitudinal sound velocities are dependent on the elastic modulus. Therefore, it is essential to further investigate the anisotropy of velocity due to the elastic modulus of these carbon allotropes is anisotropic. In this work, the directions of [1 0 0], [0 1 0] and [0 0 1] for orthorhombic crystal are investigated. The calculation of sound velocities in different directions are based from the formula as follows: Orthorhombic crystal class:

[100] vl =

C11/ [010] vt 1 =

C 66/ [001] vt 2 =

C 55/

(14)

[010] vl =

C 22/ [100] vt1 =

C 66/ [001] vt 2 =

C 44/

(15)

[001] vl =

C 33/ [100] vt1 =

C55/ [010] vt 2 =

C 44/

(16)

Monoclinic crystal class: For the [1 0 0] direction

vt = 8

C66/

(17)

Computational Materials Science 171 (2020) 109229

W.-D. Han, et al.

Fig. 7. The surface construction of bulk modulus.

Fig. 7. (continued)

9

Computational Materials Science 171 (2020) 109229

W.-D. Han, et al.

Fig. 8. The surface construction of Young's modulus.

Fig. 8. (continued)

10

Computational Materials Science 171 (2020) 109229

W.-D. Han, et al.

11

(caption on next page)

Computational Materials Science 171 (2020) 109229

W.-D. Han, et al.

Fig. 9. The projections of bulk and Young's modulus on the (0 0 1), (0 1 0) and (1 0 0) crystal planes of carbon. Graphs (a), (c) and (e) are the projections of Young's modulus E, and graphs (b), (d) and (f) are the projections bulk modulus B.

[0 1 0] and [0 0 1] directions correspond to C11, C22 and C33, respectively. Generally, the more value of Cxx (x = 1,2,3) is, the more the longitudinal sound velocities. Obviously, Imma carbon has the largest longitudinal sound velocities in [1 0 0] and [0 0 1] directions due to the highest value of C11 and C33, similarly, the value of longitudinal sound velocities in [0 1 0] for Amm2 carbon is largest for the highest value of C22 for Amm2 carbon. The velocities of the faster mode (v+) in [1 0 0] and [0 0 1] directions is the largest among those three modes of velocity for monoclinic carbon. In the [0 1 0] direction, the velocities in the longitudinal mode (vl) are larger than those in two pure transverse modes (vt+ and vt-) for monoclinic carbon.

Table 5 The density ρ (Kg/m3), sound velocity (longitudinal vl, transverse vt and average vm) (m/s), acoustic Grüneisen parameterγɑ and Debye temperature ΘD(K) of carbon. Space group

ρ(Kg/m3)

vl(m/s)

vt(m/s)

vm(m/s)

ΘD(K)

γɑ

P2221 Imma Amm2 I-4 C2/m-16 Cm-32 P21/m Cm-40 C2m-20 C2/m-28 Diamond

3362.7 3357.7 3293.3 3440.6 3336.2 3257.8 3333.4 3343.2 3364.3 3182.1 3475.6

17533.6 17300.0 17029.2 16383.9 17301.2 17029.4 17449.5 17233.1 17302.8 16316.7 17979.3

12044.2 11733.5 11639.9 11074.9 11776.7 11514.2 11930.9 11708.9 11792.4 11048.1 12307.9

13123.3 12795.9 12682.4 12084.9 12914.9 12562.8 13003.9 12768.0 12856.8 12062.7 13407.4

2159.3 2104.3 2071.8 2003.6 2119.4 2045.5 2133.5 2096.9 2115.8 1948.7 2230.5

1.011 0.906 0.877 0.919 0.893 0.918 0.870 0.900 0.889 1.152 0.817

v+ =

v =

2 + 4C 2 (C11 + C55) + (C11 15

1 2 )2 2C11 C55 + C55

2

2 + 4C 2 (C11 + C55) + (C11 15

1 2 )2 2C11 C55 + C55

2

3.6. Heat conductivities The minimum heat conductivity play a significant role for material at high temperature since the heat conductivity will drop to an extreme as the rise of temperature [56]. The minimum heat conductivity kmin of carbon is obtained based on Clarke's model [57,58] and and Cahill’model [59]. The formula can be expressed as follows: Clarke’s model:

(18)

2

(19)

For the [0 1 0] direction

vl = v+ =

v =

2 + 4C 2 (C44 + C66) + (C44 46

1 2 )2 2C44 C66 + C66

2

(C44 + C66)

2 + 4C 2 (C44 46

1 2 )2 2C44 C66 + C66

2

(21) (22)

v+ =

v =

2 + 4C 2 (C33 + C55) + (C33 35

1 2 )2 2C33 C55 + C55

2

(C33 + C55)

2 + 4C 2 (C33 35

2

1 2 )2 2C33 C55 + C55

(27)

kB 2 n3 2.48

(vl + 2vt )

(28)

where vt and vl are the transverse and longitudinal sound velocities, respectively. n is density of the number of atoms in each unit cell. The models of Clarke and Cahill can be employed to obtain the minimum heat conductivities and the results are shown in Table 7 and Fig. 11. For Clarke's model, it can be observed that the minimum heat conductivity is dominated by Ma, E and ρ, implying a larger density, a larger elastic modulus and a lower average mass of atoms benefit a higher minimum heat conductivity. Moreover, the minimum heat conductivity is also governed by the Young's modulus E due to the values of Ma are 1.176 × 10−25 for all the carbon and the density of these carbon allotropes is close to each other. The order of E is P2221 > C2/m-20 > Imma > C2/m-16 > Cm-40 > P21/m >

(23)

C44/

Ma = [M /(m . NA)]

k min =

For the [0 0 1] direction

vt =

(26)

where Ma is the mean mass of atoms in the unit cell, kB is Boltzmann’s constant, NA is the Avogadro’s number, E is Young's modulus, M is the molar mass, ρ is the density, and m is the number of atoms in a molecule. Cahill’s model:

(20)

C22/

1 1 6

k min = 0.87kB Ma 3 E 2

(24) (25)

Table 6 lists the directional sound velocities of carbon. For the orthorhombic structure, the longitudinal sound velocities in [1 0 0],

Fig. 10. Calculated Debye temperature ΘD(K) and the Grüneisen parameterγɑ (b) for carbon allotropes. 12

Computational Materials Science 171 (2020) 109229

W.-D. Han, et al.

Table 6 The anisotropic sound velocity (in m/s) of carbon. Space group P2221 Imma Amm2 I-4 C2/m-16 Cm-32 P21/m Cm-40 C2m-20 C2/m-28

Diamond

[1 0 0] [1 0 0]vl 16607.8 17661.9 17628.9 16180.3 [1 0 0] vt 12228.3 11782.1 10747.2 12002.6 11854.0 10023.6 [1 0 0] [1 0 0]vl 17463.6

[0 1 0] [0 1 0]vt1 11702.2 10601.5 11580.3 11627.7

[0 0 1]vt2 12321.3 12285.7 11485.5 10851.8

v+ 18822.8 17692.5 17120.4 16532.1 18234.3 15638.3

v13352.9 11310.8 11445.7 11753.2 11597.3 11104.9

[0 1 0]vt1 12705.9

[0 1 0]vt2 12705.9

[0 0 1]

[0 1 0]vl 17421.5 17468.6 18369.6 16265.4 [0 1 0] vl 18635.3 16443.3 17385.0 17967.3 17797.1 156548.2 [1 1 0] [1 1 0]vl 18129.0

Amm2 > Cm-32 > I-4 > C2/m-28, while the order of kmin is P2221 > P21/m > C2m-20 > C2/m-16 > Cm-40 > Imma > Amm2 > Cm-32 > I-4 > C2/m-28, the slight different order between the E and kmin can be attributed to the variety of their density and Young's modulus. However, a larger elastic modulus corresponds to a higher value of minimum thermal conductivity overall. Additionally, it can be found from Fig. 11 that the minimum thermal conductivities evaluated by the model Clarke are slightly lower than those by that of Cahill, which can be attributed to the consideration of phonon spectrum in Cahill’s model while the lacking contribution of optical phonon in in Clarke's model. Theoretically, thermal conductivity are contributed by lattice and electronic heat conductivities. Generally, the heat conductivity is dominated by electron-phonon scattering at low temperature while governed by the lattice thermal conductivity at the temperature of ground state. In addition, Debye temperature determines the lattice thermal conductivity at low temperature, indicating a higher the Debye temperature corresponds to a larger lattice thermal conductivity [60]. Therefore, P2221 and P21/m carbon have relative higher minimum thermal conductivity due to their relative larger Debye temperature. While C2/m-28 carbon exhibits the lowest thermal conductivity, which coincides well with the conclusion of Debye temperature. To show more details about heat conductivity, the model of Cahill is also employed to describe the anisotropy of heat conductivity, the expression are listed as follow:

k min =

kB 2 n3 2.48

(vl + vt1 + vt 2)

[1 0 0]vt1 11702.2 10601.5 11580.3 11627.7

[0 0 1]vt2 12505.8 11288.7 11046.1 10851.8

vt+ 13870.1 10301.2 12379.4 12426.5 12494.4 10126.5

vt11311.9 9871.2 10733.5 11040.0 11103.2 9542.6

[11–0]vt1 11736.9

[0 0 1]vt2 12705.9

[0 0 1]vl 16568.3 18357.5 16321.3 16729.7 [0 0 1] vt 12185.6 10099.6 12366.2 11489.4 11793.7 10265.4 [1 1 1] [1 1 1]vl 18345.5

[1 0 0]vt1 12321.3 12285.7 11485.5 10851.8

[0 1 0]vt2 12505.8 11288.7 11046.1 10851.8

v+ 18693.4 15929.1 17337.8 17226.9 16481.2 15463.5

v12886.7 11252.9 11526.4 11742.4 11589.0 10644.9

[11¯2]vt1 12068.5

[11¯2] vt2 12068.5

Fig. 11. The calculated minimum thermal conductivity kmin based on Clarke and Cahill model.

For orthorhombic structure, Table 7 lists the minimum heat conductivities in [1 0 0], [0 1 0] and [0 0 1] directions. Obviously, the minimum thermal conductivities is ansotropic in three directions for these carbon, which can be attributed to the difference of vl, vt1 and vt2. It can be found the value of kmin[0 1 0] is the largest except Imma carbon, while that of kmin[0 0 1] is the smallest for Amm2 and I-4 carbon, indicating that the conduction velocity of heat in [0 1 0] direction is the fastest except Imma carbon, while conduction velocity of heat in [0 0 1] direction is the slowest for Amm2 and I-4 carbon. As for Imma carbon,

(29)

where vl is longitudinal sound velocities, vt1 and vt2 are two transverse branches of sound velocities. Table 7 Calculated thermal conductivities kmin (W·m−1·K−1) of carbon. Space group

Clarke model

P2221 Imma Amm2 I-4 C2/m-16 Cm-32 P21/m Cm-40 C2m-20 C2/m-28 Diamond

Ma(10−25) 1.176 1.176 1.176 1.176 1.176 1.176 1.176 1.176 1.176 1.176 1.176

Cahill model Kmin 6.416 6.225 6.141 6.066 6.287 6.080 6.344 6.264 6.325 4.942 6.720

n(1028) 16.67 16.87 16.52 16.31 16.73 16.34 16.73 16.78 16.88 15.96 17.48

Kmin 7.210 6.933 6.677 6.404 7.170 6.669 6.938 6.888 6.976 6.003 7.417

13

Kmin[1 0 0] 6.789 6.775 6.799 6.460 – – – – – – –

Kmin[0 1 0] 6.956 6.576 6.850 6.474 7.321 6.118 6.767 6.923 6.917 6.403 –

Kmin[0 0 1] 6.917 7.006 6.492 6.422 – – – – – – –

Computational Materials Science 171 (2020) 109229

W.-D. Han, et al.

Fig. 12. Calculated total phonon density of state.

Fig. 13. (a) Calculated enthalpy E(T), (b) the partial enlargement.

the performance of thermal conductivity in [0 0 1] is the most favorable while the worst in [0 1 0]. The different performance of thermal conductivity can be attributed to the the variety of phonon mean free path in different direction. Accordingly, the phonon collision probability is

highest for P2221, Amm2 and I-4 carbons in [0 1 0] direction, while Imma carbon has the most phonon collision probability in [0 0 1] direction.

14

Computational Materials Science 171 (2020) 109229

W.-D. Han, et al.

Fig. 14. (a) Calculated Gibbs free energy G(T),(b) the partial enlargement.

Fig. 15. (a) Calculated entropy S(T),(b) the partial enlargement.

Fig. 16. (a) Calculated heat capacity at constant volume CV(T), (b) the partial enlargement.

3.7. Thermodynamic properties.

(T), Gibbs free energy G(T), entropy S(T) and heat capacity CV(T), respectively. Theoretically, the phonon can be employed to describe the dynamical stability for material. It can be observed from Fig. 12 that the phonon frequency is positive for these carbon allotropes except negative frequency of C2/m-28 carbon, indicating C2/m-28 carbon is

As the superhard materials, it is crucial to further investigate the thermodynamic properties due to high temperature produced in the process of work. Figs. 12–16 show the phonon density state, enthalpy E 15

Computational Materials Science 171 (2020) 109229

W.-D. Han, et al.

dynamically metastable structure and other carbon allotropes are dynamically stable structures. What should be pointed is that the maximum phonon frequency of diamond (38.15 THz) is the lowest among these carbon. allotropes. Moreover, the maximum phonon frequency of P21/m carbon is the largest (43.8 THz) and then is Cm-40 carbon. Generally, the shorter bond and lower value of lattice parameter correspond to the larger binding force between atoms, leading to the larger phonon frequency [61]. Thus, the binding force between the atoms of P21/m carbon is the largest and then is Cm-40 carbon. The order of maximum phonon frequency is P21/m > Cm-40 > Cm-32 = C2/m20 > C2/m-16 > P2221 > Imma > Amm2 > I-4 > Diamond. The maximum phonon frequency of Cm-32 and C2/m-20 is equal, indicating the binding force between the atoms of both the carbon allotropes is same. As is shown in Figs. 13–16, the values of enthapy E(T), entropy S(T) and volume heat capacity CV increase with the rising temperature while the value of Gibbs free energy G(T) decreases with the increase temperature. In addition, the values of enthapy E(T) and Gibbs free energy G(T) is close to zero at low temperature (0–200 K) for all these carbon allotropes. For entropy S(T) and volume heat capacity CV, the values is zero at the temperature ranging from 0 to 100 K. Normally, the volume heat capacity CV is proportional to T3 at low temperature and converges to the Dulong-Petit limit of 3NR when the temperature is higher than Debye temperature. The the volume heat capacity CV has not converged due to the Debye temperature of these carbon allotropes higher than 1000 K. Moreover, the thermodynamic properties of Cm-32 and C2/m20 is almost identical as the overlapped curve of enthapy E(T), entopyS (T), Gibbs free energy G(T) and heat capacity, which is well consistent with the analysis of phonon density. At the same temperature, the order of enthapy E(T), entropy S(T) and heat capacity is Cm-40 > C2/m28 > I-4 > Imma > Amm2 > Cm-32 = C2/m-20 > P2221 > C2/ m-16 > P21/m > Diamond, while the order is opposite in terms of Gibbs free energy G(T) and the order is Cm-40 < C2/m-28 < I4 < Imma < Amm2 < Cm-32 < C2/m-20 < P2221 < C2/m16 < P21/m < Diamond. Generally, the lower value of Gibbs free energy corresponds to the better thermal stability [62,63]. Hence, Cm40 carbon has the strongest thermal stability while P21/m possesses the weakest thermal stability. Furthermore, the thermal stability of these carbon allotropes is better than that of diamond, indicating the newlydesigned carbon crystals can be favorable superhard materials, overcoming the defect of weak thermal stability for diamond.

(5) The order of the elastic anisotropy for these carbon allotropes is C2/ m-28 > Imma > C2/m-16 > P21/m > C2/m-20 > Cm32 > Cm-40 > Amm2 > P2221 > I-4. (6) The sound velocities of these carbon allotropes are anisotropic, the Debye temperature and thermal conductivity is the largest for P2221 carbon while the least for C2/m-28. (7) The thermal stability of these newly-designed carbon allotropes is better than that of diamond and the order is Cm-40 > C2/m28 > I-4 > Imma > Amm2 > Cm-32 = C2/m20 > P2221 > C2/m-16 > P21/m > Diamond, Cm-40 carbon has the strongest thermal stability. In addition, the C2/m-28 carbon is dynamically metastable structure. Acknowledgements The authors gratefully acknowledge the financial supports from the National Natural Science Foundation of China (No. 51665036), the Natural Science Foundation of Jiangxi Province (No. 20152ACB20014). References [1] C.M. Sung, M. Sung, Carbon nitride and other speculative superhard materials, Mater. Chem. Phys. 43 (1996) 1. [2] J. Haines, J. Leger, G. Bocquillon, Synthesis and design of superhard materials, Annu. Rev. Mater. Res. 31 (2001) 1. [3] M. Hu, Z.S. Zhao, F. Tian, A.R. Oganov, Q.Q. Wang, M. Xiong, Compressed carbon nanotubes: a family of new multifunctional carbon allotropes, Sci. Rep. 3 (2013) 1331. [4] M. Amsler, J.A. Flores-Livas, L. Lehtovaara, F. Balima, Crystal structure of cold compressed graphite, Phys. Rev. Lett. 108 (2012) 065. [5] J.T. Wang, C. Chen, Y. Kawazoe, Orthorhombic carbon allotrope of compressed graphite: Ab initio calculations, Phys. Rev. B 85 (2012) 033. [6] C.Y. He, L.Z. Sun, C.X. Zhang, J.X. Zhong, Two viable three dimensional carbon semiconductors with an entirely sp 2 configuration, Phys. Chem. Chem. Phys. 15 (2013) 680. [7] C.Y. He, L.Z. Sun, C.X. Zhang, X.Y. Peng, K.W. Zhang, J.X. Zhong, Four superhard carbon allotropes: a first-principles study, Phys. Chem. Chem. Phys. 14 (2012) 8410. [8] Y.A. Kvashnina, A.G. Kvashnin, P.B. Sorokin, Investigation of new superhard carbon allotropes with promising electronic properties, J. Appl. Phys. 114 (2013) 183708. [9] F. Tian, X. Dong, Z. Zhao, J. He, H.T. Wang, J. Phys.Condens. Matter 24 (2012) 165504. [10] A.R. Oganov, C.W. Glass, Crystal structure prediction using ab initio evolutionary techniques: principles and applications, J. Chem. Phys. 124 (2006) 244704. [11] Y. Wang, J.E. Panzik, B. Kiefer, K.M. Lee, Crystal structure of graphite under roomtemperature compression and decompression, Sci. Rep. 2 (2012) 520. [12] Y. Wang, K.M. Lee, From soft to superhard: fifty years of experiments on coldcompressed Graphite, J. Superhard Mater. 34 (2012) 360. [13] X.X. Zhang, Y.C. Wang, J. Lv, C.Y. Zhu, First-principles structural design of superhard materials, J. Chem. Phys. 138 (2013) 114101. [14] R. Martoňák, A.R. Oganov, C.W. Glass, Crystal structure prediction and simulations of structural transformations: metadynamics and evolutionary algorithms, Phase Transit. 80 (2007) 277. [15] Q. Zhu, A.R. Oganov, A.O. Lyakhov, Evolutionary metadynamics: a novel method to predict crystal structures, CrystEngComm 14 (2012) 3596. [16] S.E. Boulfelfel, Q. Zhu, A.R. Oganov, Novel sp 3 forms of carbon predicted by evolutionary metadynamics and analysis of their synthesizability using transition path sampling, J. Superhard Mater. 34 (2012) 350. [17] G.J. Finkelstein, P.K. Dera, S. Jahn, A.R. Oganov, C.M. Holl, Y. Meng, T.S. Duffy, Phase transitions and equation of state of for sterite to 90 GPa from single-crystal Xray diffraction and molecular modeling, Am. Miner. 99 (2014) 35. [18] J.Y. You, X.Y. Ma, et al., Carboneyane: a nodal line topological carbon with sp-sp2sp3 chemical bonds, Carbon 152 (2019) 909–914. [19] M.J. Xing, B.H. Li, Z.T. Yu, Q. Chen, C2/m-carbon: structural, mechanical, and electronic properties, J. Mater. Sci. 50 (2015) 7104. [20] J.P. Perdew, A. Ruzsinszky, G.I. Csonka, O.A. Vydrov, G.E. Scuseria, L.A. Constantin, X. Zhou, K. Burke Perdew, et al., Reply. Phys. Rev. Lett. 100 (2008) 136406. [21] Y.L. Page, P. Saxe, Symmetry-general least-squares extraction of elastic data for strained materials from ab initio calculations of stress, Phys. Rev. B 65 (2002) 104104. [22] Y. Li, Y. Gao, B. Xiao, T. Min, Z. Fan, First-principles study on the structural, elastic, and thermodynamics properties of Ni3X (X: Al, Mo, Ti, Pt, Si, Nb, V, and Zr) intermetallic compounds, J. Alloys Comp. 89 (2016) 1. [23] D. Vanderbilt, Ab initio study of ferroelectric domain walls in PbTiO3, Phys. Rev. B 65 (2002) 263. [24] B. Hammer, L.B. Hansen, J.K. Nørskov, Nitrogen adsorption and dissociation on Fe (111), Phys. Rev. B 59 (1999) 7413. [25] G.P. Francis, M.C. Payne, Structural complexity in grain boundaries with covalent

4. Conclusion Based on the density functional theory, the ab-initio method is employed to systematically and deeply investigated the structural property, elastic, elastic anisotropy, sound velocities, thermal conductivities and thermodynamic properties of carbon allotropes. The equilibrium lattice parameters after optimization are consistent well with the results at previous literature. (1) The obtained cohesive energies and formation enthalpies indicate that these carbon allotropes in different structures are energetically and structurally stable and the order of structural stability is P2221 > C2/m-20 > Imma > I-4 > Cm-32 > Amm2 > Cm40 > C2/m-16 > P21/m > C2/m-28. (2) These carbon allotropes are mechanically stable due to the calculated elastic constants of single crystal satisfy the criterions of mechanical stability. (3) A larger bulk modulus not entirely correspond to a higher hardness, while there is a relationship between shear modulus and hardness that a larger shear modulus indicates a higher hardness. (4) All these carbon phases can be favorable candidate for superhard materials due to the hardness is more than 70GPa, and the order of hardness is P2221 > P21/m > C2/m-20 > C2/m-16 > Cm40 > Imma > Amm2 > Cm-32 > I-4 > C2/m-28. 16

Computational Materials Science 171 (2020) 109229

W.-D. Han, et al. bonding, J. Phys. Condens. Matter 2 (1990) 4395. [26] B.G. Pfrommer, M. Côte, S.G. Louie, M.L. Cohen, Ab initio study of silicon in the R8 phase, Phys. Rev. B 56 (1997) 6662. [27] M.L. Petrescu, Boron nitride theoretical hardness compared to carbon polymorphs, Diamond Relat. Mater. 13 (2004) 1848. [28] Q.Y. Fan, Q. Wei, C.C. Chai, H.Y. Yan, M.G. Zhang, Z.Z. Lin, Structural, mechanical, and electronic properties of P3m1-BCN, J. Phys. Chem. Solids 79 (2015) 89. [29] M. Grimsditch, E.S. Zouboulis, A. Polian, Elastic constants of boron nitride, J. Appl. Phys. 76 (1994) 832. [30] M. Born, K. Huang, Dynamical Theory of crystal Lattices, Oxford University Press, Oxford, 1998. [31] L. Ma, Y. Duan, R. Li, Structural, elastic and electronic properties of C14-type Al2M (M=Mg, Ca, Sr and Ba) Laves phases, Phys. B 507 (2017) 147. [33] W. Voigt, Handbook of Crystal Physics, Taubner, Leipzig, 1928. [34] A. Reuss, Z. Angew, Math. Mech 9 (1929) 49. [35] R. Hill, Proc. Phys. Soc. A 65 (1952) 349. [36] Z.J. Wu, E.J. Zhao, H.P. Xiang, X.F. Hao, X.J. Liu, J. Meng, Publisher's Note: crystal structures and elastic properties of superhard IrN[sub 2] and IrN[sub 3] from first principles, Phys. Rev. B 76 (2007) 054115. [37] M.M. Wu, L. Wen, B.Y. Tang, L.M. Peng, W.J. Ding, First-principles study of elastic and electronic properties of MgZn2 and ScZn2 phases in Mg–Sc–Zn alloy, J. Alloys Comp. 506 (2010) 412. [38] S.F. Pugh, Philos. Mag. 45 (1954) 823. [39] D.G. Pettifor, Theoretical predictions of structure and related properties of intermetallics, Mater. Sci. Technol. 8 (1992) 345349. [40] B. Huang, Y.H. Duan, Y. Sun, M.J. Peng, S. Chen, Electronic structures, mechanical and thermodynamic properties of cubic alkaline-earth hexaborides from first principles calculations, J. Alloys Compd. 635 (2015) 213. [41] Y. Tian, B. Xu, Z. Zhao, Microscopic theory of hardness and design of novel superhard crystalsInt, J. Refract. Metals Hard Mater. 33 (2012) 93. [42] J.J. Gilman, Chemistry and Physics of Mechanical Hardness, John Wiley & Sons, New Jersey, 2009. [43] J.B. Levine, S.H. Tolbert, R.B. Kaner, Advancements in the search for superhard ultra-incompressible metal borides, Adv. Funct. Mater. 19 (2009) 3519. [44] D.M. Teter, Computational alchemy: the search for new superhard materials, MRS Bull. 23 (1998) 22. [45] A.L. Ivanovskii, Mechanical and electronic properties of diborides of transition 3d–5d metals from first principles: toward search of novel ultra-incompressible and superhard materials, Prog. Mater. Sci. 57 (2012) 184. [46] D.H. Chung, W.R. Buessem, F.W. Vahldiek, S.A. Mersol (Eds.), Anisotropy in Single Crystal Refractory Compound, Plenum, New York, 1968. [47] S.I. Ranganathan, M. Ostoja-Starzewski, Phys. Rev. Lett. 101 (2008) 055504.

[48] P. Ravindran, L. Fast, P.A. Korzhavyi, B. Johansson, Ab initio prediction of the structural, electronic, elastic and thermodynamic properties of the tetragonal ternary intermetallics. XCu2Si2 (X = Ca, Sr), J. Appl. Phys. 84 (1998) 4891. [49] J.F. Nye, Physical Properties of Crystals, Oxford University Press, Oxford, 1985. [50] R.F.S. Hearmon, A.A. Maradudin, Phys. Today 14 (1961) 48. [51] D. Music, A. Houben, R. Dronskowski, J.M. Schneider, Ab initio study of ductility in M_ 2 AlC (M= Ti, V, Cr), Phys. Rev. B 75 (2007) 174102. [52] S. Chen, Y. Sun, Y.H. Duan, B. Huang, M.J. Peng, Phase stability, structural and elastic properties of C15-type Laves transition-metal compounds MCo2 from firstprinciples calculations, J. Alloys Compd. 630 (2015) 202. [53] H.W. Shou, Y.H. Duan, Anisotropic elasticity and thermal conductivities of (ɑ, β, γ)LiAlSi2O6 from the first-principles calculation, J. Alloys Compd. 756 (2018) 40. [54] C.M.I. Okoye, Structural, elastic and electronic structure of LiCu2Si, LiCu2Ge and LiAg2Sn intermetallic compounds, Comput. Mater. Sci. 92 (2014) 141. [55] F. Arab, F.A. Sahraoui, K. Haddadi, A. Bouhemadou, L. Louail, Phase stability, mechanical and thermodynamic properties of orthorhombic and trigonal MgSiN2: an ab initio study, Phase Transit. 89 (2016) 480. [56] Y. Shen, D.R. Clarke, P.A. Fuierer, Anisotropic thermal conductivity of the Aurivillus phase, bismuth titanate (Bi4Ti3O12): a natural nanostructured superlattice, Appl. Phys. Lett. 93 (2008) 102907. [57] D.R. Clarke, Materials selection guidelines for low thermal conductivity thermal barrier coatings, Surf. Coating Technol. 163–164 (2003) 67. [58] D.R. Clarke, C.G. Levi, Annu. Rev. Mater. Res. 33 (2003) 383. [59] D.G. Cahill, S.K. Watson, R.O. Pohl, Lower limit to the thermal conductivity of disordered crystals, Phys. Rev. B 46 (1992) 6131. [60] P. Li, L.S. Ma, M.M. Peng, B.P. Shu, Y.H. Duan, Elastic anisotropies and thermal conductivities of WB2 diborides in different crystal structures: a first-principles calculation, J. Alloys Compd. 747 (2017) 905. [61] E. Deligoz, K. Colakoglu, The first principles study on boron bismuth compounds, Comput. Mater. Sci. 39 (2007) 533. [62] N. Sun, X. Zhang, J. Qin, J. Ning, S. Zhang, S. Liang, R. Liu, First principle study of elastic and thermodynamic properties of ZrZn2 and HfZn2 under high pressure, J. Appl. Phys. 115 (2014) 083514. [63] D.W. Zhou, J.S. Liu, S.H. Xu, Thermal stability and elastic properties of Mg3Sb2 and Mg3Bi2 phases from first-principles calculations, Phys. B 405 (2010) 2863.

Further reading [32] Z. Zhao, B. Xu, X.F. Zhou, L.M. Wang, B. Wen, J. He, Z. Liu, H.T. Wang, Y. Tian, Novel superhard carbon: C-centered orthorhombic C8, Phys. Rev. Lett. 107 (2011) 215502.

17