Ab-initio calculations of carbon and boron nitride allotropes and their structural phase transitions using periodic coupled cluster theory

TL;DR

This study uses advanced quantum chemical methods to analyze the stability and phase transitions of carbon and boron nitride allotropes, providing insights into their thermodynamic properties and structural transformations.

Contribution

It demonstrates the effectiveness of coupled cluster theories in accurately predicting phase stability and transitions in BN and carbon allotropes, benchmarking against experimental data.

Findings

Coupled cluster theories predict degeneracy of BN phases at 0 K.

Quantum chemical methods improve accuracy over density functional theories.

Stacking order influences formation of meta-stable BN and carbon phases.

Abstract

We present an ab-initio study of boron nitride as well as carbon allotropes. Their relative thermodynamic stabilities and structural phase transitions from low- to high-density phases are investigated. Pressure-temperature phase diagrams are calculated and compared to experimental findings. The calculations are performed using quantum chemical wavefunction based as well as density functional theories. Our findings reveal that predicted energy differences often depend significantly on the choice of the employed method. Comparison between calculated and experimental results allows for benchmarking the accuracy of various levels of theory. The produced results show that quantum chemical wavefunction based theories allow for achieving systematically improvable estimates. We find that on the level of coupled cluster theories the low- and high-density phases of boron nitride become thermodynamically degenerate at 0 K. This is in agreement with recent experimental findings, indicating that cubic boron nitride is not the thermodynamically stable allotrope at ambient conditions. Furthermore we employ the calculated results to assess transition probabilities from graphitic low-density to diamond-like high-density phases in an approximate manner. We conclude that the stacking order of the parent graphitic material is crucial for the possible formation of meta-stable wurtzite boron nitride and hexagonal carbon diamond also known as lonsdaleite.

Tables

Table 1. Table 1: Structural parameters of the eight atomic monoclinic type unit cell ( β 𝛽 \beta  =  γ 𝛾 \gamma  =  90   ° times 90 degree 90\text{\,}\mathrm{\SIUnitSymbolDegree} ) including interlayer distance d 𝑑 d , the forming bond distance R 𝑅 R marked red in the transition state of Fig.  3  –  6 and its corresponding bond length in the other structures, atomic volume V 𝑉 V and density ρ 𝜌 \rho . See text for further details.

Carbon	G-ABC	G-AB	G-AA	pc-TS	pw-TS	bw-TS	l-pc-TS	$c$-D	$h$-D
$\left|\overrightarrow{a}\right|$ [Å]	$2.446$	$2.446$	$2.445$	$2.465$	$2.461$	$2.471$	$2.456$	$2.498$	$2.484$
$\left|\overrightarrow{b}\right|$ [Å]	$4.236$	$4.236$	$4.236$	$4.270$	$4.263$	$4.164$	$4.229$	$4.326$	$4.302$
$\left|\overrightarrow{c}\right|$ [Å]	$6.735$	$6.590$	$7.171$	$4.899$	$4.780$	$4.835$	$4.705$	$4.326$	$4.137$
$\alpha$	$77.9$	$90.0$	$90.0$	$73.1$	$90.0$	$90.0$	$80.3$	$70.5$	$90.0$
d [Å]	$3.292$	$3.295$	$3.585$	$2.344$	$2.390$	$2.418$	$2.319$	$2.039$	$2.069$
R [Å]	$3.292$	$3.295$	$3.585$	$2.058$	$2.108$	$2.100$	$2.472$	$1.530$	$1.549$
V [${\AA }^3{\mathrm{atom}}^{-1}$]	$8.528$	$8.535$	$9.285$	$6.167$	$6.268$	$6.220$	$6.021$	$5.510$	$5.526$
$\rho$ [$g{\mathrm{cm}}^{-3}$]	$2.417$	$2.415$	$2.220$	$3.343$	$3.289$	$3.315$	$3.424$	$3.742$	$3.731$
Boron nitride	$r$-BN	BN-AB	$h$-BN	pc-TS	pw-TS	bw-TS	l-pc-TS	$c$-BN	$w$-BN
$\left|\overrightarrow{a}\right|$ [Å]	$2.488$	$2.488$	$2.488$	$2.506$	$2.506$	$2.511$	$2.513$	$2.532$	$2.524$
$\left|\overrightarrow{b}\right|$ [Å]	$4.309$	$4.309$	$4.310$	$4.341$	$4.341$	$4.219$	$4.282$	$4.386$	$4.371$
$\left|\overrightarrow{c}\right|$ [Å]	$7.069$	$6.458$	$6.491$	$4.932$	$4.742$	$4.900$	$4.707$	$4.386$	$4.176$
$\alpha$	$77.5$	$90$	$90$	$72.9$	$90$	$90$	$79.5$	$70.5$	$90$
d [Å]	$3.229$	$3.229$	$3.246$	$2.357$	$2.371$	$2.450$	$2.315$	$2.068$	$2.088$
R [Å]	$3.229$	$3.229$	$3.246$	$2.088$	$2.102$	$2.128$	$2.447$	$1.551$	$1.564$
V [Å3]	$9.246$	$8.653$	$8.700$	$6.411$	$6.448$	$6.489$	$6.227$	$5.742$	$5.759$
$\rho$ [$g{\mathrm{cm}}^{-3}$]	$2.230$	$2.382$	$2.370$	$3.216$	$3.197$	$3.177$	$3.311$	$3.590$	$3.580$

Table 2. Table 2: Relative difference to experimental values Lynch and Drickamer ( 1966 ); Yoshiasa  et al. ( 2003 ); Slack and Bartram ( 1975 ); Furthmüller  et al. ( 1994 ); Nagakubo  et al. ( 2013 ) of lattice vectors | a → | → 𝑎 |\vec{a}| and | c → | → 𝑐 |\vec{c}| for stable structures of carbon an BN. In the case of w 𝑤 w -BN , h ℎ h -BN , h ℎ h -D  and G-AB we compare to a range of experimental values.

	$c$-BN	$w$-BN	$h$-BN	$c$-D	$h$-D	G-AB
$\left|\overrightarrow{a}\right|$ [%]	0.94	$1.03\text{–}1.04$	$0.47\text{–}0.63$	0.97	$1.04\text{–}0.96$	$0.58\text{–}0.62$
$\left|\overrightarrow{c}\right|$ [%]		$0.63\text{–}1.21$	$2.53\text{–}2.55$		$-0.62\text{–}1.10$	$0.78\text{–}1.75$

Table 3. Table 3: Energy differences in meV per atom for carbon structures as obtained by various levels of theory.

System	LDA	PBE	PBE+MBD	SCAN	PBE0	B3LYP	HSE06	HF	MP2-TA-FS	CCSD-TA-FS	CCSD(T)-TA
$c$-D – G-ABC	$-12$	$121$	$52$	$94$	$45$	$206$	$49$	$110$	$-19\pm 29$	$-31\pm 5$	$-14$
$h$-D – $c$-D	$25$	$24$	$25$	$29$	$27$	$30$	$27$	$37$	$30\pm 14$	$35\pm 10$	$33$
pc-TS – G-ABC	$329$	$423$	$372$	$422$	$437$	$530$	$435$	$671$	$210\pm 46$	$426\pm 1$	$385$
bw-TS – G-ABC	$395$	$488$	$438$	$495$	$512$	$604$	$509$	$772$	$334\pm 49$	$529\pm 20$	$481$
pw-TS – G-ABC	$362$	$452$	$404$	$452$	$474$	$560$	$470$	$710$	$238\pm 61$	$469\pm 14$	$423$
l-pc-TS – G-ABC	$367$	$470$	$413$	$473$	$473$	$581$	$470$	$685$	$350\pm 41$	$469\pm 16$	$441$

Table 4. Table 4: Energy differences in meV per atom for boron nitride structures as obtained by various levels of theory.

System	LDA	PBE	PBE+MBD	SCAN	PBE0	B3LYP	HSE06	HF	MP2-TA-FS	CCSD-TA-FS	CCSD(T)-TA
$c$-BN – $r$-BN	$-56$	$64$	$-13$	$47$	$30$	$160$	$31$	$222$	$-58\pm 21$	$14\pm 11$	$2$
$w$-BN – $c$-BN	$18$	$17$	$18$	$20$	$19$	$20$	$19$	$23$	$23\pm 7$	$24\pm 4$	$22$
pc-TS – $h$-BN	$169$	$250$	$195$	$240$	$249$	$329$	$249$	$409$	$173\pm 17$	$243\pm 10$	$221$
bw-TS – $h$-BN	$214$	$292$	$240$	$286$	$294$	$372$	$294$	$454$	$233\pm 42$	$299\pm 30$	$281$
pw-TS – $h$-BN	$179$	$259$	$207$	$251$	$258$	$338$	$258$	$412$	$188\pm 18$	$256\pm 11$	$234$
l-pc-TS – $h$-BN	$202$	$293$	$230$	$288$	$286$	$381$	$287$	$461$	$213\pm 32$	$283\pm 19$	$259$

Equations

$$G$$

$$F_{\mathrm{vib}}$$

$$G^{\mathrm{CCSD}(\mathrm{T})}(V,T)\approx G^{\mathrm{LDA}}(V,T)-E_{\mathrm{tot}}^{\mathrm{LDA}}(V_0)+E_{\mathrm{tot}}^{\mathrm{CCSD}(\mathrm{T})}(V_0),$$

$$f=\exp (\frac{-\Delta G_{\mathrm{forward}}}{RT})-\exp (\frac{-\Delta G_{\mathrm{back}}}{RT})$$

Full text

I Introduction

The pressure-temperature phase diagrams of carbon and boron nitride reflect a delicate balance between weak and strong interatomic interactions. Although covalent bonds are the main source of their large cohesive energies, the accumulation of weak van der Waals interactions contributes significantly to the relative stability of their low- and high-density phases. Furthermore vibrational effects play a crucial role in the temperature dependence of the equilibrium phase boundary. Altogether this makes the prediction of phase diagrams and structural phase transition pathways a challenging task for modern electronic structure theories. In this work we seek to investigate boron nitride as well as carbon allotropes using various approximate electronic structure theories and compare theoretical with experimental findings. The aim is to benchmark their accuracy and help interpreting experimental results better if possible. To this end we employ a range of approximate density functional theories (DFT) and quantum chemical wavefunction based methods.

During the last decades approximate exchange and correlation (XC) density functionals have made significant progress in becoming more accurate and predictive for the description of interatomic interactions while keeping a high level of computational efficiency that allows for studying systems containing several hundreds of atoms routinely. The so-called Jacob’s ladder describes a ladder of approximations for the XC energy using increasingly complex as well as in general more accurate methods Perdew and Schmidt (2001). These rungs include functionals based on the local density approximation (LDA) Ceperley and Alder (1980); Perdew and Zunger (1981), the generalized gradient approximation (GGA) Perdew et al. (1996a), the meta generalized gradient approximation (mGGA) Sun et al. (2015) and hybrid functionals Perdew et al. (1996b); Becke (1993); Heyd et al. (2003); Krukau et al. (2006). The latter include a fraction of (screened) exact exchange energies and have a computational cost that is comparable to Hartree–Fock theory. However, all the functionals mentioned above suffer from shortcomings that are despite many efforts difficult to remedy Cohen et al. (2012). In the context of the present work a significant shortcoming is the inaccurate description of long range van der Waals interactions. In order to describe van der Waals and related interatomic interactions more accurately in the framework of approximate XC density functionals, a wide variety of dispersion corrections has been developed. As a consequence of the large number of available density functionals and corrections, there are numerous ground state energy functionals that could be considered in the present work, of which we have only chosen a small selection.

As a complement to the treatment of exchange and correlation on the level approximate density functionals, the computationally significantly more expensive quantum chemical wavefunction based theories are becoming more popular for the study of periodic systems Booth et al. (2013); Yang et al. (2014); Müller and Paulus (2012); Grüneis et al. (2011); Usvyat (2013); Nolan et al. (2009); Hirata et al. (2004); Rościszewski et al. (1999); Stoll and Doll (2012); Ren et al. (2012); Ben et al. (2013); Neuhauser et al. (2013); Neufeld and Thom (2017); Hermann and Schwerdtfeger (2008); Schwerdtfeger et al. (2010); McClain et al. (2017); Boese and Sauer (2016); Ochi and Tsuneyuki (2015); Usvyat et al. . This can partly be attributed to the increase in their computational efficiency, due to methodological developments, and to their ability to predict exchange and correlation energies in a systematically improvable manner. Quantum chemical methods constitute a hierarchy, which starting from the one-particle Hartree–Fock (HF) approximation, allows for a systematic treatment of the quantum many-body effects. The simplest form of such correlated methods is the second-order Møller–Plesset (MP2) perturbation theory Møller and Plesset (1934). The next level of theory that achieves a significantly improved trade-off between accuracy and computational cost is based on the coupled cluster ansatz for the many-electron wavefunction Cizek (1966). Coupled cluster singles and doubles theory provides a compelling framework of infinite-order approximations in the form of an exponential of cluster operators Bartlett and Musiał (2007). The coupled-cluster singles and doubles (CCSD) method where the triples are treated in a perturbative way, termed as CCSD(T), achieves chemical accuracy in the description of many molecular properties and is sometimes referred to as the gold standard method Raghavachari et al. (1989).

In this work we seek to investigate the accuracy of the electronic structure theories mentioned above for carbon and boron nitride allotropes. To this end we compare predicted ground state energy differences as well as calculated phase diagrams to experimental findings. Experimentally several phases have been synthesized as single crystals Kubota et al. (2008); Lu et al. (2015); Taniguchi and Yamaoka (2001); Austerman et al. (1967); Kanda (2000) or as powder Nagakubo et al. (2013); Matsui et al. (1981). Single crystals are usually synthesized in a closed chamber over a longer time period crystallizing from a solution. For synthesizing metastable structures the samples are put under static pressure and heated electrically or by laser Endo et al. (1994); Yagi et al. (1992); Britun et al. (1993); Kurdyumov et al. (1996); Taniguchi et al. (1997a). Another way of transforming samples into metastable phases is shock wave synthesis Wheeler and Lewis (1975); Sato et al. (1982). With these materials thermodynamic characterization can be performed Wagman (1945); Madelung et al. (2002); Solozhenko (1993, 1995); Jeong and Lee (2013); Wise et al. (1966) by determining relative enthalpy, entropy and heat capacity. These properties can be used to compare with the calculated energy differences and construct phase diagrams. Furthermore the phase diagrams can also be obtained by observing phase transitions directly Bundy et al. (1996); Clarke and Uher (1984); Corrigan and Bundy (1975); Eremets et al. (1998); Onodera et al. (1981); Sachdev et al. (1997); Wills (1985); Fukunaga (2000). In the present work we also investigate pressure-driven concerted phase transition pathways. In particular we study activation barrier heights for the transformation from low-density to high-density systems considering small unit cells that contain a few atoms at most. These models are far from realistic conditions under which phase transitions occur in experiment. Temperature-driven kinetic effects and catalysts are needed in practice to observe phase transitions close to the equilibrium phase boundary Berman and Simon (1955); Bundy et al. (1961); Bundy (1980); Bundy and Wentorf (1963). Hot liquid metals can be used to dissolve graphite and diamond will precipitate at the cooler region. However, the aim of the current work is to explore the accuracy of various electronic structure theories for transition states occurring in these phase transitions and provide a qualitative description of the various possible phase transition mechanisms. We believe that the small supercells considered are sufficient to describe these effects qualitatively correct.

This paper is organized as follows. Section II provides a description of the considered stable and metastable structures followed by an overview of the considered phase transition pathways. The employed structures can also be found in the Supplementary informations Sup . The results section IV summarizes the calculated ground state energy differences of the (meta-)stable structures and their activation energies at $0\text{\,}\mathrm{K}$ for the investigated phase transition pathways. By calculating the Gibbs energies, pressure-temperature phase diagrams are predicted and compared to experiment. Furthermore we assess the temperature and pressure dependence of the activation energies. Based on these results the experimentally observed phase transitions will be reviewed. Furthermore the existence of the wurtzite structure of carbon and boron nitride will be discussed. In the course of the discussion of these results we will assess the accuracy of the various approximate electronic structure theories.

II Crystal structures and phase transition pathways for carbon and boron nitride

II.1 (Meta-)stable structures

In this work we consider the most abundant crystal structures of carbon and boron nitride: the graphitic and diamond-like phases. Fig. 1 illustrates the corresponding structures. The low-density phases are graphitic with all atoms being $sp^{2}$ bonded and arranged in the planar honeycomb lattice with different stackings: AA, AB and ABC as depicted in Fig. 1. The G-AB and G-ABC have been observed experimentally for carbon and can be transformed into each other by translation of the layers. For boron nitride the stable low-density phase is (hexagonal) $h$-BN. $h$-BN exhibits an AA’ stacking order, indicating the other atom types for lattice sites on top of each other in the direction of stacking as shown in Fig. 1a. We note that the AA (G-AA) stacking is unstable for carbon.

The high-density phases of carbon and boron nitride are diamond-like. All atoms in the considered diamond-like phases can be assigned to chair or boat conformations of six-membered rings and two different stacking orders. For carbon and boron nitride the most stable high-pressure phases are cubic diamond ($c$-D) and zinc blende ($c$-BN), respectively. $c$-D and $c$-BN consist of six-membered rings of $sp^{3}$ bonded atoms in the chair conformation with an ABC stacking order as shown in Fig. 1e. There is a second high-density phase that is referred to as wurtzite for boron nitride and hexagonal diamond ($h$-D) for carbon with AA stacking order, exhibiting chair and boat conformations as illustrated in Fig. 1f. The mirror image of A is A, with the mirror parallel to the layer. The chair and boat conformations are parallel and perpendicular to the stacking, respectively. $h$-D is also known as lonsdaleite and serves as marker for shock impact events.

The crystal structures for C can easily be derived from BN, by substituting all B- and N-atoms with C-atoms. This makes (rhombohedral) $r$-BN (Fig. 1c), $c$-BN (Fig. 1e) and $w$-BN (Fig. 1f) equivalent to G-ABC, $c$-D and $h$-D, respectively.

Experimentally single crystals have been reported for $h$-BN Kubota et al. (2008); Lu et al. (2015), $c$-BN Taniguchi and Yamaoka (2001), G-AB Austerman et al. (1967) and for $c$-D Kanda (2000). There exist no single crystals for $w$-BN, but the XRD measurements show no sign of $c$-BN in the samples and a small amount ($<$2%) of the starting material $h$-BN Nagakubo et al. (2013). The synthesis of $r$-BN results in fibrous micro-crystals, but show no mixture with $h$-BN Matsui et al. (1981). The $h$-D phase has been investigated in a number of theoretical as well as experimental studies in the past Fahy et al. (1986, 1987); Tateyama et al. (1996); Bundy et al. (1996). However, recent experimental studies indicate that the previously believed samples of $h$-D are in fact $c$-D crystals that contain a large number of twins and stacking faults, creating x-ray diffraction patterns similar to the hypothetical $h$-D Németh et al. (2014).

II.2 Structural phase transition pathways

We now discuss the investigated structural phase transition pathways. For the present study we keep the computational cost of the coupled cluster theory calculations low by restricting ourselves to transition state geometries that contain at most four atoms in the unit cell. Similar transition states have already been investigated in Refs. Fahy et al. (1986, 1987); Tateyama et al. (1996); Dong et al. (2013); Wentzcovitch et al. (1988); Wang et al. (2011); Yu et al. (2003). To drive a transition from the low-density graphitic phases to the high-density diamond-like phases the application of pressure is needed. Under pressure the c-axis of the graphitic phase experiences a much larger compression than the other axes and will therefore be referred to as the compression axis. At high pressures the planar structure of graphite splits. Fig. 2 depicts two basic mechanisms by which the splitting of the planar six-membered rings present in the honeycomb lattice occurs. The mechanisms are referred to as buckling or puckering. Buckling and puckering creates the boat and chair conformation of six-membered rings, respectively. We employ the following naming convention for transition states. The first letter refers to the puckering (p) or buckling (b) mechanism and the second letter refers to the cubic (c) or wurtzite (w) structure corresponding to the final state of the considered transition.

pc-TS is the transition state in the G-ABC to $c$-D transition via the puckering mechanism. We note that the cubic phase contains only six-membered rings in the chair conformation. Therefore it is reasonable to consider this one transition state (Fig. 3b) only.

The wurtzite structure contains six-membered rings in the chair and boat conformation. In the pw-TS the chair conformation is perpendicular and the boat conformation is parallel to the compression axis (Fig. 4b), whereas in the bw-TS the orientations are switched (Fig. 5b). By comparing the resulting structures from the corresponding transition pathways one can see that the c-axis of $w$-BN is rotated by 90∘ (Fig. 4c and 5c).

We also consider the transition from $w$-BN to $c$-BN, which occurs in a stepwise layer-to-layer rearrangement through 4H intermediate structures and will be referred to as l-pc-TS in this work Britun et al. (1993). During this transformation the boat conformation along the c-axis of the $w$-BN structure transforms into the chair conformation, while the six-membered rings perpendicular to the c-axis break apart (e.g. atoms 3+8 bottom layers) and rebond differently (e.g. atom 4+7 bottom layers) (Fig. 6a \ce-¿[6b] 6c with plane (001)${}_{\mathrm{w}}\parallel$ (111)c and direction [10$\overline{1}$0]${}_{\mathrm{w}}\parallel$ [112]c).

In total we consider four different transition states including pc-TS, pw-TS, bw-TS and l-pc-TS for carbon and boron nitride allotropes. In the case of carbon all lattice sites are occupied by the same atomic species.

The geometries of the transition states have been determined as follows. For the four atomic unit cells of the pc-TS and pw-TS three degrees of freedom were considered: the intra and inter layer bond distance and the angle in between. The first-order saddle point on the corresponding potential energy surface defines the transition state geometry as well as its energy. Due to the small number of considered degrees of freedom a sweeping algorithm was sufficient to determine the pc-TS and pw-TS. Determining the bw-TS and l-pc-TS is slightly more complicated due to the larger number of degrees of freedom. For bw-TS the $z$-coordinate of the interlayer bond distance $R$, the horizontal lattice vectors in Fig. 5b, the out-of-plane displacement of the atoms and the in-plane coordinates of all atoms were considered as degrees of freedom. The bw-TS was obtained applying a sweeping algorithm to all degrees of freedom except for the in-plane coordinates of all atoms and the lattice vectors that were optimized by relaxing the structures accordingly for a given unit cell. The l-pc-TS was obtained in a similar manner. A sweeping algorithm was employed to calculate energies for all coordinates and cell parameters interpolated linearly between the initial ($w$-BN) and the final ($c$-BN) structure. For each interpolated structure all cell parameters and atomic positions were allowed to relax while keeping only the vertical-coordinate of the atoms in Fig. 6a frozen.

II.3 Crystal lattice parameters

Tab. 1 summarizes the lattice parameters of the employed geometries for the low- and high-density phases and transition states for carbon and boron nitride allotropes. These parameters have been optimized using DFT in the LDA. This is necessary because forces are not yet implemented in the employed coupled cluster theory code. We believe that the LDA provides sufficiently accurate structures compared to experiment that allow for an unbiased comparison between the employed electronic structure theories and to experiment. We note that the LDA lattice parameters deviate by about $1\text{\,}\mathrm{\char 37\relax}$ only from experiment even for the lattice vector $|\vec{c}|$ parallel to the compression axis as summarized in Tab. 2. The only exception is $h$-BN, where the deviation is slightly larger. To allow for a direct comparison between the lattice parameters of the (meta-)stable structures as well as transition states we consider an eight atomic unit cell with monoclinic lattice vectors $\vec{a}$, $\vec{b}$ and $\vec{c}$ such that $|\vec{a}|$ $\neq$ $|\vec{b}|$ $\neq$ $|\vec{c}|$, $\beta$ = $\gamma$ = $90\text{\,}\mathrm{\SIUnitSymbolDegree}$ and $\alpha$ can also be $90\text{\,}\mathrm{\SIUnitSymbolDegree}$. The $\vec{a}\vec{b}$ plane can be seen in Fig. 1d. The length of $\vec{a}$ corresponds to the width of a honeycomb ring and the vector points out of plane in Fig. 3 – 6. $\vec{b}$ points from left to right and $\vec{c}$ from bottom to top and spans across two layers. The ratio $|\vec{a}|$:$|\vec{b}|$ is 1:$\sqrt{3}$ for all structures except for bw-TS and l-pc-TS. In these cases the ratio is larger by up to 3 %. We point out that in l-pc-TS $R$ is much larger compared to the other transition states. However, $R$ refers to all bonds between the layers in the other transition states, whereas this not the case for l-pc-TS. In l-pc-TS additional interlayer bonds exist with a bond length of $1.654\text{\,}\mathrm{\SIUnitSymbolAngstrom}$ (1+6 or 2+5 in Fig. 6b). For carbon (BN) the average of these two bond lengths is $2.063\text{\,}\mathrm{\SIUnitSymbolAngstrom}$ ($2.050\text{\,}\mathrm{\SIUnitSymbolAngstrom}$), respectively and comparable to $R$ of the other transitions states. The employed structures can be found in the Supplementary informations Sup .

III Methods

III.1 Density functional and Hartree–Fock theory

All electronic structure calculations have been performed using the projector augmented wave (PAW) method Blöchl (1994) as implemented in the Vienna $ab-initio$ simulation package (VASP) Kresse and Hafner (1994); Kresse and Furthmüller (1996). We present results obtained using some of the most widely-used $ab-initio$ methods to approximate the exchange and correlation energy in the framework of DFT. These methods include the LDA functional as parametrized by Perdew-Zunger Ceperley and Alder (1980); Perdew and Zunger (1981), the GGA functional as parametrized by Perdew-Burke-Ernzerhof (PBE) Perdew et al. (1996a), the dispersion corrected PBE functional using the many-body dispersion energy method (PBE+MBD) Tkatchenko et al. (2012); Bucko et al. (2016), the meta-GGA as parametrized for the SCAN functional Sun et al. (2015) and the hybrid density functionals PBE0 Perdew et al. (1996b), Becke-3-parameter-Lee-Yang-Parr (B3LYP) Becke (1993) and Heyd-Scuseria-Ernzerhof (HSE06) Heyd et al. (2003); Krukau et al. (2006). We note that we have chosen only a small selection of functionals that could be considered.

The B 2$s^{2}2p^{1}$, N 2$s^{2}2p^{3}$ and C 2$s^{2}2p^{2}$ states have been treated as valence states in all calculations.

The geometries have been relaxed until the forces on all atoms are smaller than ${10}^{-5}\text{\,}\mathrm{eV}\text{\,}{\mathrm{\SIUnitSymbolAngstrom}}^{-1}$. The total energies have been converged using the self-consistent field approach to within ${10}^{-8}\text{\,}\mathrm{eV}$. For all DFT and HF calculations we employed a 16 atom supercell and a $4\times 4\times 4$ Monkhorst-Pack $k$-point mesh. The corresponding supercell structures are summarized in the supplementary informations. The kinetic energy cutoff for the plane wave basis set was set to $1000\text{\,}\mathrm{eV}$. We note that smaller kinetic energy cutoffs would have sufficed but these calculations do not consitute a computational bottle neck compared to the more expensive coupled cluster theory calculations.

The phonon calculations have been performed using the Phonopy code Togo and Tanaka (2015), creating the displacements within a $2\times 2\times 2$ supercell of the 16 atom cell. The forces are calculated using VASP and the LDA. These calculations employed a kinetic energy cutoff of $800\text{\,}\mathrm{eV}$ and a $2\times 2\times 2$ $k$-mesh.

III.2 Quantum chemical wavefunction theories

Results obtained using post-Hartree–Fock methods have been converged with respect to several computational parameters including energy cutoffs defining the plane wave basis sets, the number of virtual orbitals and the $k$-mesh. We have employed kinetic energy cutoffs of 500 eV for definining the orbital plane wave basis set and 300 eV for definining an auxiliary basis set that is used in the calculation of electron repulsion integrals required in post-HF methods. For the twist averaging technique we have employed a 4$\times$4$\times$4 $k$-mesh. Furthermore 14 unoccupied orbitals have been used per atom. These parameters ensure a convergence of the energy difference between graphite and diamond to within a few meV per atom. The same parameters have been employed in Ref. Gruber et al. (2018). For the virtual orbital space we employ MP2 natural orbitals that are obtained using a procedure outlined in Ref. Grüneis et al. (2011). Our estimates of the remaining basis set incompleteness error indicate that the energy difference between carbon diamond and graphite should be converged to within approximately $4\text{\,}\mathrm{meV}\text{\,}{\mathrm{atom}}^{-1}$ as shown in Fig.2 of Ref. Gruber et al. (2018). A similar level of accuracy is expected for the other energy differences.

III.2.1 Finite size errors

We stress that the convergence of ground state energies obtained using post-HF methods such as MP2 theory with respect to the employed $k$-mesh or supercell size is slower than for their DFT counterparts. Wavefunction based methods account for non-local electronic correlation effects explicitly and therefore the observed interatomic interactions such as van der Waals forces lead to a slower rate of convergence with respect to the employed $k$-mesh. Recently we have introduced a finite size correction scheme that allows for accelerating the rate of convergence for periodic systems. We refer to results obtained using the finite size correction scheme by employing the following naming convention. Corrected MP2 and CCSD results are referred to as MP2-TA-FS and CCSD-TA-FS, respectively. TA and FS stand for twist averaging and an interpolation method, respectively. For the perturbative triples (T) correction on top of CCSD-TA-FS, we employ the twist averaging technique only. As such CCSD(T)-TA refers to CCSD-TA-FS plus the (T)-TA contribution. The improved $k$-mesh and supercell size convergence of CCSD-TA-FS was demonstrated and discussed in Ref. Gruber et al. (2018). If not stated otherwise, all MP2, CCSD and CCSD(T) results in this work include the finite size corrections.

In the present work we employ box plots to depict finite size errors in Figs. 8, 9, 10 11, 12 and 13. The box plots show the distribution of the obtained results for a set of different $k$-meshes. We stress that finite size errors must not be confused with stochastic errors. However, small error bars indicate that results are not affected significantly by the size of the employed $k$-mesh and can therefore be considered converged with respect to the $k$-mesh density. The employed box plots mark the largest and lowest value, the second and third quartile as well as the mean value. The bars in Figs. 8, 9, 10 11, 12 and 13 mark the value with the largest $k$-point mesh. Some results for certain $k$-meshes differ significantly from the other $k$-meshes and were marked as outlier with a ”+“. We find that very anisotropic Brillouin zone sampling using $k$-meshes with one $k$-point along one direction only yields results that are considered outliers.

For the results depicted in Figs. 8 and 11 we employ two-atomic unit cells and the following $k$-meshes: $2\times 2\times 2$, $2\times 2\times 4$, $3\times 3\times 2$ and $3\times 3\times 3$.

For the calculations of the barrier heights shown in Figs. 9, 10 and Figs. 12, 13 we have employed four and eight atomic unit cells. The following $k$-meshes were employed to sample the Brillouin zone of the four atomic unit cells G-ABC/$h$-BN, pc-TS, pw-TS (bw-TS and l-pc-TS in parentheses): $2\times 2\times 2$ ($1\times 2\times 4$), $3\times 3\times 1$ ($1\times 3\times 3$) and $3\times 3\times 2$ ($2\times 3\times 3$). We note that two transition states (bw-TS and l-pc-TS) need four atoms per layer to be described correctly and have only one layer. Therefore the $c$-axis is just half as long and the $a$-axis is doubled compared to the other cells. Consequently the $k$-mesh was adjusted appropriately. For the eight atomic unit cell we have employed a $1\times 3\times 2$ $k$-mesh.

III.3 Thermodynamic properties

For the calculation of pressure-temperature phase diagrams we need to compute the Gibbs energies ($G$) of all phases. $G$ is defined as the sum of the ground state energy, as obtained from DFT or a similar approach ($E$tot), all entropic contributions and the $pV$ term. The vibrational ($F$vib) contribution is the largest entropy related contribution. Using the finite displacement method a phonon density of state ($D(\omega)$) has been obtained for the frequency range ($\omega$). This phonon density of state contains the vibrational information for a specific volume ($V$) and can be used to calculate $F_{\mathrm{vib}}$ at any temperature ($T$) with the Planck ($h$) and Boltzmann constant (kb) such that

[TABLE]

To account for volume expansion during temperature increase the quasi harmonic approximation is used. For the diamond-like phases isotropic expansion is assumed. The Gibbs energy is calculated for at least five different volumes. The universal equation of state (EOS) Vinet et al. (1989) has been used to find the minimum of $G$ with respect to the volume for a given temperature and pressure. For all graphite-like phases and transition states anisotropic expansion along the c-axis has been included. One parameter changes the unit cell isotropically and a second one changes only the c-axis. A fourth order polynomial fit is used to interpolate the EOS between the sampling points. This increases the number of data points and also the accuracy. Using the procedure described above we have calculated the Gibbs energies employing the LDA for a wide range of pressures and temperatures. However, due to the computational cost involved it is currently not possible to perform the same calculations on the level of CCSD(T). The CCSD(T) Gibbs energies have therefore been approximated using the following expression

[TABLE]

where $V_{0}$ corresponds to the LDA (equilibrium) volume of the (meta-)stable allotropes and the transition state geometries. As such the volume and temperature dependence of the CCSD(T) Gibbs energy is approximated using LDA, which achieves sufficiently accurate descriptions of the phonons and bulk moduli for the purpose of the present study.

III.4 Phase transition probability

In this work we will approximate the probability that a phase transition occurs using the activation energy only and disregard kinetic effects. The activation energy ($\Delta G$) is the difference between the Gibbs energies of the transition and initial states. The Gibbs energy can be calculated as described in the previous subsection. Close to the equilibrium phase boundary the back and forward reaction has to be taken into account and the probability $f$ depends exponentially on $\Delta G$

[TABLE]

with $R$ as the gas constant. A comparable ansatz was published in Refs. Yafei et al. (1994); Wang and Yang (1999).

Fig. 7a to 7f shows the calculated CCSD(T) phase transition probabilities and activation energies for carbon and a selection of transition states. We note that the behavior of the transition probabilities at low temperatures mostly arises from the explicit dependence of $f$ on the temperature in the exponent rather than the temperature dependence of $\Delta G_{\mathrm{forward}}$ or $\Delta G_{\mathrm{back}}$. Sec. IV.4 provides a more detailed discussion of the obtained results.

IV Results and discussion

We now turn to the discussion of the obtained DFT, HF and post-HF results. The following section is organized as follows. We first summarize the energy differences obtained using different methods for carbon (sec. IV.1) and boron nitride (sec. IV.2) allotropes, respectively. Subsequently a comparison between results obtained for carbon and boron nitride will be drawn in sec. IV.3. Sec. IV.4 employs the calculated ground state energies on the level of CCSD(T) theory and the DFT results to predict the pressure-temperature phase diagrams of carbon and boron nitride. In sec. IV.5 we review experimentally observed phase transitions and compare to the produced theoretical results. Section IV.6 focuses on the hexagonal form of diamond.

IV.1 Carbon allotropes

Fig. 8 depicts the electronic ground state energy differences (\Delta E=E_{\text{c-D}}-E_{G-ABC}) between carbon diamond ($c$-D) and graphite (G-ABC) obtained using a range of DFT and quantum chemical wavefunction based theories. Experimentally G-AB is the most stable form of graphite. However, G-ABC and G-AB are degenerate to within a few meV per atom. Positive and negative energy differences in Fig. 8 indicate the thermodynamic stability of graphite and diamond, respectively. We stress that these calculations employ the DFT-LDA relaxed structures and that further relaxation effects of the respective functionals are not taken into account. Further relaxation effects can be significant for functionals that fail to describe the interlayer binding in graphite. However, this section will focus on benchmarking the accuracy of the employed functionals for a fixed geometry only. In passing we note, however, that for comparison we have repeated the DFT calculations summarized in Fig. 8 using geometries relaxed on the level of the PBE+MBD functional. The corresponding energy differences did not change by more than five percent as a result of the small changes in the employed geometries.

The grey bar in Fig. 8 corresponds to the experimental estimate of the ground state energy difference corrected for zero-point vibrational energies (ZPVEs). The experimental value of the difference in the Gibbs energy between graphite and diamond is $25\text{\,}\mathrm{meV}\text{\,}{\mathrm{atom}}^{-1}$ and has been obtained from the heat of combustion and extrapolation to $0\text{\,}\mathrm{K}$ using the heat capacity Wagman (1945). Therefore the latter value includes ZPVEs that stabilize graphite compared to diamond. To allow for a direct comparison between experiment and theory we have removed $9\text{\,}\mathrm{meV}\text{\,}{\mathrm{atom}}^{-1}$ ZPVE contributions (estimated using DFT-LDA) from the experimental energy difference.

We now turn to the discussion of the energy differences in Fig. 8 obtained using XC functionals in the framework of DFT. LDA underestimates the energy difference, predicting diamond to be more stable than graphite by $12\text{\,}\mathrm{meV}\text{\,}{\mathrm{atom}}^{-1}$. On the level of the GGA using the PBE functional we find that the stability of graphite is significantly overestimated by almost $100\text{\,}\mathrm{meV}\text{\,}{\mathrm{atom}}^{-1}$ compared to experiment. This overestimation is partly reduced by including dispersion effects on the level of MBD or by switching to the SCAN functional, yielding energy differences of $52\text{\,}\mathrm{meV}\text{\,}{\mathrm{atom}}^{-1}$ (PBE+MBD) and $94\text{\,}\mathrm{meV}\text{\,}{\mathrm{atom}}^{-1}$ (SCAN), respectively. Furthermore hybrid functionals such as PBE0 or HSE06 constitute a further improvement compared to SCAN, overestimating the stability of graphite compared to experiment by a few $10\text{\,}\mathrm{meV}\text{\,}{\mathrm{atom}}^{-1}$ only. However, we note that the B3LYP hybrid functional does not follow this trend and gives the worst agreement with experiment out of all theories considered in the present study. Therefore a systematic improvability of the employed XC functionals with respect to their rung and computational cost can not be achieved in the present case. Furthermore one conclusion of the above findings is that non-van der Waals corrected higher-level functionals (PBE, SCAN, PBE0 and HSE06) predict the graphitic phase to be more stable than diamond, whereas the inclusion of van der Waals corrections (MBD) can reverse their ordering. Indeed we note in passing that in contrast to HSE06, HSE06+MBD predicts $c$-D being more stable than G-ABC by $30\text{\,}\mathrm{meV}\text{\,}{\mathrm{atom}}^{-1}$.

We now turn to the discussion of the results obtained using wavefunction based theories as depicted in Fig. 8. HF predicts graphite to be more stable than diamond by approximately $100\text{\,}\mathrm{meV}\text{\,}{\mathrm{atom}}^{-1}$, albeit neglecting dispersion effects that play an important role in the interlayer binding of graphite. We note that due to the neglect of these contributions, HF would predict the isolated graphene sheets to be more stable than graphite. Second-order Møller-Plesset (MP2) perturbation theory corresponds to the next level of wavefunction based method and predicts diamond to be slightly more stable than graphite by approximately $19\text{\,}\mathrm{meV}\text{\,}{\mathrm{atom}}^{-1}$. However, the finite size errors of the obtained MP2 results are significant as indicated by the box plot, which is described in Sec. III.2.1. The $k$-point mesh convergence using CCSD-TA-FS theory is much faster compared to MP2-TA-FS theory as indicated by the smaller error bar. We find that CCSD-TA-FS theory predicts diamond to be more stable than graphite by $31\text{\,}\mathrm{meV}\text{\,}{\mathrm{atom}}^{-1}$. Including the perturbative triples contribution to CCSD-TA-FS theory yields an even better agreement with experiment albeit predicting diamond to be slightly more stable than graphite by $14\text{\,}\mathrm{meV}\text{\,}{\mathrm{atom}}^{-1}$. In passing we note that DMC has been used in Ref. Shin et al. (2014) to predict an energy difference in almost perfect agreement with experiment, whereas the random-phase approximation (RPA) predicts both allotropes to be exactly degenerate Lebègue et al. (2010). The good agreement between DMC and experiment is partly fortuitous due to remaining errors from the stochastic sampling, the fixed-node approximation and the employed pseudo-potentials. However, the agreement between the high-level methods such as wavefunction based theories DMC and CCSD(T), the RPA results and experiment to within a few ten $\mathrm{meV}\text{\,}{\mathrm{atom}}^{-1}$ is encouraging. The remaining finite size and basis set errors in CCSD(T) theory calculations do not allow for predicting which carbon allotrope is more stable, although we can conclude that they are expected to be degenerate to about 10–20 $\mathrm{meV}\text{\,}{\mathrm{atom}}^{-1}$ including ZPVE. Our findings indicate that quantum chemical wavefunction based theories allow for a systematic improvability of the predicted energy differences as one increases the level of theory ranging from HF, MP2, CCSD to CCSD(T).

Having demonstrated that CCSD(T) theory is expected to yield accurate energy differences for the thermodynamically most stable carbon allotropes we now seek to investigate the pressure-driven transition pathways introduced in Sec. II.2. To this end we focus on the activation barrier height that is defined as the difference in the electronic ground state energy between graphite and the corresponding transition state $E_{A}=E_{\rm TS}-E_{\rm G-ABC}$. The considered transition states are referred to as pc-TS, pw-TS, bw-TS and l-pc-TS. The activation barriers can not be compared to experimental observations directly but serve as theoretical benchmark systems and qualitative models for realistic phase transitions. Fig. 9 depicts the difference of calculated activation barrier heights between various methods and CCSD(T) (including finite size corrections); for example, $\Delta E_{A}^{\rm LDA}=E_{A}^{\rm LDA}-E_{A}^{\rm CCSD(T)-TA}$. The depicted results confirm well-known trends for the accuracy of DFT methods, assuming that CCSD(T) theory can be considered an accurate benchmarking reference for the activation barrier height. LDA underestimates the activation barrier heights for all investigated transition states by $50100\text{\,}\mathrm{meV}\text{\,}{\mathrm{atom}}^{-1}$, showing that this level of theory suffers from larger errors in the description of XC energies for transition states compared to initial and final states (G-ABC and $c$-D). Including gradient corrections on the level of the PBE functional improves the agreement with CCSD(T) theory noticeably, yielding overestimated activation energies with errors smaller than $50\text{\,}\mathrm{meV}\text{\,}{\mathrm{atom}}^{-1}$ for all four transition states. This is in contrast to PBE results for molecular activation barrier heights in the gas phase that are in general underestimated Cohen et al. (2012). However, we believe that the overestimation of the PBE barriers for the studied solids is caused by the neglect of interatomic van der Waals forces, which play an important role for the present systems. We stress that the inclusion of dispersion effects to PBE on the level of PBE+MBD theory yields again underestimated activation energies that agree with CCSD(T) to within $50\text{\,}\mathrm{meV}\text{\,}{\mathrm{atom}}^{-1}$. The SCAN functional yields barrier heights that are almost identical to our PBE findings. Furthermore the inclusion of non-local exchange in the PBE0 and HSE06 hybrid functionals yields on average slightly larger barrier heights. We note that adding the MBD effect (from the difference between PBE and PBE+MBD calculations) to these hybrid functionals would yield barrier heights in almost perfect agreement with CCSD(T) theory. On the other hand, we find that the B3LYP hybrid functional yields overestimated barrier heights, exhibiting errors on a scale of more than $100\text{\,}\mathrm{meV}\text{\,}{\mathrm{atom}}^{-1}$.

We now turn to the discussion of activation barrier heights calculated using wavefunction based theories starting with the HF method. Our findings are depicted in Fig. 9 and show that HF yields strongly overestimated barrier heights with errors on the scale of almost $300\text{\,}\mathrm{meV}\text{\,}{\mathrm{atom}}^{-1}$ compared to CCSD(T). This trend is known from molecular quantum chemistry and can be explained by the fact that HF neglects electronic correlation effects, which are in general larger in the transition state compared to the initial and final state of most chemical reactions. Accounting for electronic correlation effects on the level of MP2 theory corrects for this tendency although it yields underestimated barriers by about $100200\text{\,}\mathrm{meV}\text{\,}{\mathrm{atom}}^{-1}$ for all transition states. We note that the box plots of the MP2 results in Fig. 9 also indicate that the remaining finite size errors for these estimates are significant. We attribute this to the observation that some transition states exhibit a metallic character in DFT calculations (l-pc-TS and bw-TS) and that MP2 theory suffers from severe shortcomings in metals such as $k$-point mesh divergence Shepherd and Grüneis (2013). CCSD results for the barrier heights constitute a substantial improvement over MP2 findings, exhibiting errors compared to CCSD(T) that are smaller than $50\text{\,}\mathrm{meV}\text{\,}{\mathrm{atom}}^{-1}$. Furthermore we note that the box plot for CCSD-TA-FS results is significantly smaller, indicating that the remaining finite size errors are below a few $10\text{\,}\mathrm{meV}\text{\,}{\mathrm{atom}}^{-1}$. From these findings we conclude that quantum chemical wavefunction based theories including MP2 and CC theories have the potential of achieving results for activation barrier heights in solid-solid phase transitions with systematically improvable accuracy.

Having discussed the accuracy of DFT and wavefunction based methods for predicting the activation barrier heights, we now seek to address the question: which transition states are energetically the most favorable? This is an important question because it affects through which transition state a pressure-driven phase transition proceeds and which (meta-)stable carbon allotrope will be the outcome. Fig. 10 depicts the energy differences of the activation barrier heights with respect to the pc-TS for the respective electronic structure theories. Unequivocally all theories predict the pc-TS to be the energetically most favorable transition state, implying that the puckering mechanism is expected to play the most important role in pressure-driven graphite to diamond transitions. In passing we note that our LDA results are comparable to previous work Fahy et al. (1987); Tateyama et al. (1996) and that the energy difference between the boat and chair conformation of graphane is $55\text{\,}\mathrm{meV}\text{\,}{\mathrm{atom}}^{-1}$, favoring the chair conformation Sofo et al. (2007). As regards the ordering of the remaining transition states (bw-TS, pw-TS, and l-pc-TS), we find that the DFT methods shown in Fig. 10 predict all very similar orderings. bw-TS is energetically the least favorable transition state, whereas pw-TS and l-pc-TS agree to within a few meV per atom, except for the PBE, SCAN and B3LYP functionals that predict the pw-TS to be slightly more favorable in energy than the l-pc-TS. In the case of results obtained using wavefunction based methods depicted in Fig. 10 we find that the bw-TS corresponds to the largest activation barrier height and that pw-TS and l-pc-TS agree to within the remaining finite size errors. However, MP2 theory deviates from this trend by predicting equally large activation barrier heights for the bw-TS and l-pc-TS, making the pw-TS the second most favorable transition state. However, we stress that MP2 results are perhaps not meaningful due to the metallic character of some transition states. From the above results we conclude that interatomic van der Waals forces play a minor role in the ordering of the respective transition states. Furthermore the ordering is already correctly described on the level of the LDA to the XC functional.

Tab. IV.1 summarizes all energy differences discussed above for the (meta-)stable carbon allotropes and the transition states. Furthermore the table also lists the energy difference between $c$-D and $h$-D that is predicted by all methods to be about $30\text{\,}\mathrm{meV}\text{\,}{\mathrm{atom}}^{-1}$.