High-Temperature Thermodynamics Modeling of Graphite

: We present high-temperature thermodynamic properties for graphite from ﬁrst-principles anharmonic theory. The ab initio electronic structure is obtained from density-functional theory coupled to a lattice dynamics method that is used to model anharmonic lattice vibrations. This combined approach produces free energies and speciﬁc heats for graphite that compare well with available experiments and results from models that empirically represent experimental data, such as CALPHAD. We show that anharmonic theory for the phonons is essential for accurate thermodynamic quantities above about 1000 K.


Introduction
Carbon atoms can form materials in three varieties, namely, graphite, diamond, and fullerenes. Chemically, the carbon bonding is different in these varieties involving sp 2 (graphite) and sp 3 (diamond) orbitals bonding, while fullerenes feature mostly sp 2 , but with a small admixture of sp 3 . Consequently, carbon polymorphs exhibit contrasting material properties that make them useful for diverse applications in areas ranging from casting and molds to electronics and bioengineering [1]. Graphite, the primary focus of this study, consists of layers of carbon sheets (graphene sheets) that are only weakly bonded with each other causing great anisotropy in the material. Thus, thermal and electronic conductivity, strength, and chemical reactivity, are properties that depend very strongly on the direction in the material. Significant anisotropic properties are rather anomalous in elemental solids, and for graphite they can be leveraged in electronic components, lubricants, pencils, and casting mold materials [2]. Thermodynamically, graphite has a high melting point, modest thermal expansion, and high specific heat, making it relevant for composites and structural materials and as heat shields in certain applications, including nuclear reactors [3].
Academically, carbon is of course a very important element and one of the most abundant, with fundamental properties for life and biochemistry. In planetary science, high-pressure forms such as diamond are relevant; however, we target the graphite phase that can exist at high temperatures but only at moderate pressures [4]; see Figure 1. We are particularly interested in thermodynamics at high temperatures but are limiting ourselves to ambient pressure in this study.
Our investigation relies on first-principles theory methods that are combined with an anharmonic lattice dynamics approach. With this methodology, we provide an ab initio model that accounts for strong anharmonic phonon-phonon interactions that are often important at high temperatures. Density-functional theory (DFT) has certainly been applied for carbon and graphite electronic structures in the past, e.g., [5][6][7][8][9][10][11][12], but we find no studies that target high-temperature thermodynamics in a way that incorporates substantial lattice anharmonicity. Much has been written on carbon systems, and it is not possible to provide a complete overview here. Instead, we recommend a comprehensive review of carbon and graphene [13] as a good starting point. Our investigation relies on first-principles theory methods that are combined with an anharmonic lattice dynamics approach. With this methodology, we provide an ab initio model that accounts for strong anharmonic phonon-phonon interactions that are often important at high temperatures. Density-functional theory (DFT) has certainly been applied for carbon and graphite electronic structures in the past, e.g., [5][6][7][8][9][10][11][12], but we find no studies that target high-temperature thermodynamics in a way that incorporates substantial lattice anharmonicity. Much has been written on carbon systems, and it is not possible to provide a complete overview here. Instead, we recommend a comprehensive review of carbon and graphene [13] as a good starting point.
We address high-temperature anharmonic phonons in graphite by coupling DFT first-principles electronic structure with a self-consistent lattice dynamics method able to capture anharmonic phonon-phonon interactions. Our results are compared to more approximate quasi-harmonic theory and thermodynamic models that have been constrained by experimental data in a thermodynamically consistent fashion. As one would expect, we will show that anharmonicity cannot be ignored when the temperature rises to thousands of degrees in graphite.
In the next section, we introduce our computational approaches followed by a presentation of our thermodynamic results in Section 3, and then a conclusion and discussions in Section 4.

Computational Approaches
Because we attempt to model graphite at high temperatures, we must consider the potential importance of anharmonic atomic vibrations. Although very convenient, quasiharmonic (QH) modeling tends to break down at high temperatures above approximately 1000 K depending on the material in question. Our goal here is to calculate Gibbs free energies and heat capacities up to 4000 K using parameter-free ab initio theory and explore the necessity of including anharmonic contributions to the lattice dynamics.
The conventional procedure to study temperature effects in materials from first principles is to perform quantum molecular dynamics (QMD) simulations. In principle, QMD answers all thermodynamic questions, but the computational simulations tend to be very costly, and consequently, it is difficult to address greater domains of a materials phase space. Alternatively, we chose an efficient lattice dynamics approach, namely the self-consistent ab initio lattice dynamics (SCAILD) method [14,15], that we recently employed for some actinide carbides and nitrides [16] with very encouraging results. Initially, the approach focused on dealing with dynamically unstable phonons at low temperatures [14] but it has more recently been applied for free-energy calculations [16,17]. SCAILD is efficient, because it constrains phonons along normal-mode directions of commensurate vibrations, which have amplitudes that are mean-square displacements determined by the We address high-temperature anharmonic phonons in graphite by coupling DFT firstprinciples electronic structure with a self-consistent lattice dynamics method able to capture anharmonic phonon-phonon interactions. Our results are compared to more approximate quasi-harmonic theory and thermodynamic models that have been constrained by experimental data in a thermodynamically consistent fashion. As one would expect, we will show that anharmonicity cannot be ignored when the temperature rises to thousands of degrees in graphite.
In the next section, we introduce our computational approaches followed by a presentation of our thermodynamic results in Section 3, and then a conclusion and discussions in Section 4.

Computational Approaches
Because we attempt to model graphite at high temperatures, we must consider the potential importance of anharmonic atomic vibrations. Although very convenient, quasiharmonic (QH) modeling tends to break down at high temperatures above approximately 1000 K depending on the material in question. Our goal here is to calculate Gibbs free energies and heat capacities up to 4000 K using parameter-free ab initio theory and explore the necessity of including anharmonic contributions to the lattice dynamics.
The conventional procedure to study temperature effects in materials from first principles is to perform quantum molecular dynamics (QMD) simulations. In principle, QMD answers all thermodynamic questions, but the computational simulations tend to be very costly, and consequently, it is difficult to address greater domains of a materials phase space. Alternatively, we chose an efficient lattice dynamics approach, namely the self-consistent ab initio lattice dynamics (SCAILD) method [14,15], that we recently employed for some actinide carbides and nitrides [16] with very encouraging results. Initially, the approach focused on dealing with dynamically unstable phonons at low temperatures [14] but it has more recently been applied for free-energy calculations [16,17]. SCAILD is efficient, because it constrains phonons along normal-mode directions of commensurate vibrations, which have amplitudes that are mean-square displacements determined by the temperature. The QMD calculations, on the other hand, are not constrained, and thus become much more computationally demanding.
The SCAILD method requires the calculation of forces on atoms that are thermally displaced from the equilibrium positions defined by the crystal structure that must be represented by a supercell. Graphite consists of graphene sheets that are monolayers of carbon atoms arranged in a hexagonal fashion, and these are stacked in a layered graphite structure with four atoms per unit cell and a relatively large c/a axial ratio of 2.725.
For SCAILD, we applied a (5 × 5 × 2) supercell of the 4-atom hexagonal structure for a total of 200 atoms in the supercell. Comparative studies of a smaller 96-atom supercell showed less than a 5% difference in the lattice free energy at the highest temperatures, suggesting that the computational supercell error is less than 5% in the free energy. Three strategic atomic volumes were chosen for each temperature, T = 500, 1000, 1500, 2000, 2500, 3000, 3500, and 4000 K, and about 150 iterations of the self-consistent cycle of SCAILD converged the free energies and heat capacities. For each of these iterations, a DFT calculation was performed to determine the required atomic forces. Hence, in terms of the lattice dynamics, around 3600 DFT calculations were completed for graphite.
The weak bonding between the graphene sheets, resulting in a large c/a axial ratio in graphite, is not very accurately handled by the conventional DFT approach. The reason for this is that van der Waals interactions are known to play a role in these bonds, and they are not included in DFT. These interactions can, in principle, be calculated at significant computational expense and complexity [5,18]. Fortunately, a reasonable and very efficient workaround solution is available that mitigates the need to include van der Waals contributions to DFT for graphite. The procedure is to apply a local density approximation (LDA) of the DFT electron exchange and correlation function in conjunction with freezing the c/a axial ratio in the graphite structure to its experimental value [7]. This strategy is also consistent with the SCAILD methodology because the SCAILD supercell box is assumed to be rigid with a frozen c/a axial ratio.
For the Gibbs (or Helmholtz at ambient pressure) free energy, F(V,T), we formulated the following equation: F lat represents the contribution from the lattice (phonons) and F el the electronic total energy, excitation energy, and entropy. Analogous to this equation, we express the specific heat at constant pressure, C p , adding these contributions: The first (lattice) term is calculated within the SCAILD approach, while the other term is directly evaluated from the DFT electronic structure of the four-atom graphite unit cell. For more complete and explicit descriptions of these terms, we refer to our recent publication on thermodynamics for actinide systems [16]. We utilize two DFT electronic-structure methods to calculate the contributions needed for Equations (1) and (2). The lattice pieces (F lat and C p lat ) are determined from the lattice dynamics, and the necessary atomic force calculations are performed by employing the pseudopotential plane-wave Vienna ab initio simulation package (VASP) and the projectoraugmented-wave (PAW) method [19][20][21]. The computational constraints for VASP are defined by an energy cut-off of 400 eV and a (2 × 2 × 2) grid for a total of 12 k points for the 200-atom supercell. The remaining terms in Equations (1) and (2) are calculated with an all-electron approach, namely, the full-potential linear muffin-tin orbitals (FPLMTO) method [22]. In this case, all six carbon electrons are included in the band structure, which is obtained from a set up with a double basis set, i.e., with two energy parameters for each orbital. For the 4-atom graphite cell, a (12 × 12 × 4) grid and a total of 69 k points are included. All electronic-structure calculations (VASP and FPLMTO) assume the LDA for electron exchange and correlation, a fixed c/a axial ratio (2.725) for the hexagonal graphite structure, and no spin-orbit coupling. The only reason we divided our DFT calculations between the two methods was because forces on a 200-atom cell can be far more effectively calculated within VASP. The FPLMTO, however, is more accurate and for the small four-atom unit-cell employed for FPLMTO, computational effort is not an issue.
Previous DFT investigations, for example [7,10,23,24], have shown that this theory is able to relatively accurately reproduce room-temperature phonon dispersions for graphite. Figure 2 shows our calculated phonon dispersions for graphite. For the phonons, we applied the small displacement method [25] (which is also the first step of a self-consistent SCAILD calculation) with four configurations of atoms, displacing a magnitude of about 3% of the lattice constant. The atomic volume is determined by the minimum total energy (8.70 Å 3 /atom), subject to the constraint on the c/a axial ratio. The resulting theoretical volume is relatively close to the known room-temperature atomic volume (8.746 Å 3 /atom) of graphite.

issue.
Previous DFT investigations, for example [7,10,23,24], have shown that this theory is able to relatively accurately reproduce room-temperature phonon dispersions for graphite. Figure 2 shows our calculated phonon dispersions for graphite. For the phonons, we applied the small displacement method [25] (which is also the first step of a self-consistent SCAILD calculation) with four configurations of atoms, displacing a magnitude of about 3% of the lattice constant. The atomic volume is determined by the minimum total energy (8.70 Å 3 /atom), subject to the constraint on the c/a axial ratio. The resulting theoretical volume is relatively close to the known room-temperature atomic volume (8.746 Å 3 /atom) of graphite. The rather good agreement with experiments for phonons at low temperature, displayed in Figure 2, suggests that our computational approach to lattice dynamics is appropriate for investigating the thermodynamics of graphite.
In addition to our DFT model, we compared our results with those from the CAL-PHAD approach [28,29]. CALPHAD is a phase-based computational thermodynamics method that incorporates a phase diagram and thermodynamic data from experiments and first-principles calculations when available [30,31], to model the evolution of the Gibbs free energy of a phase as function of composition, temperature, and pressure. In the present study, we used the Thermo-Calc software with a thermodynamic database [32] to calculate the Gibbs energy and heat capacity of graphite.

Results
By combining our calculated electronic and lattice dynamics contributions, we were able to construct the Gibbs free energy (or Helmoltz free energy at zero pressure), see Equation (1). This free energy is a fundamental quantity in thermodynamics, but it is typically not directly measured. Thus, rather than validating our ab initio free energy explicitly with experiments, we compared our results with other theories and modeling that have been developed for graphite. Figure 3 shows our first-principles free energy together with quasi-harmonic (QH) theory, CALPHAD [32], and an equation of state (EOS) model [33]. The rather good agreement with experiments for phonons at low temperature, displayed in Figure 2, suggests that our computational approach to lattice dynamics is appropriate for investigating the thermodynamics of graphite.
In addition to our DFT model, we compared our results with those from the CALPHAD approach [28,29]. CALPHAD is a phase-based computational thermodynamics method that incorporates a phase diagram and thermodynamic data from experiments and firstprinciples calculations when available [30,31], to model the evolution of the Gibbs free energy of a phase as function of composition, temperature, and pressure. In the present study, we used the Thermo-Calc software with a thermodynamic database [32] to calculate the Gibbs energy and heat capacity of graphite.

Results
By combining our calculated electronic and lattice dynamics contributions, we were able to construct the Gibbs free energy (or Helmoltz free energy at zero pressure), see Equation (1). This free energy is a fundamental quantity in thermodynamics, but it is typically not directly measured. Thus, rather than validating our ab initio free energy explicitly with experiments, we compared our results with other theories and modeling that have been developed for graphite. Figure 3 shows our first-principles free energy together with quasi-harmonic (QH) theory, CALPHAD [32], and an equation of state (EOS) model [33].
We calculated the quasi-harmonic free energy in Figure 3 from the GIBBS2 package [34][35][36][37] using the energies determined by FPLMTO calculations of graphite. For all FPLMTO calculations and QH modeling, we made the same approximation to keep the axial c/a ratio constant at an experimental value. Although QH modeling for graphite has been researched, we were not able to find any free-energy data in the literature. Notably, in Figure 3, the CALPHAD and EOS models are in excellent agreement, which was expected considering that both are thermodynamically consistent representations of the experimental data. Furthermore, the figure shows a good agreement between our first-principles anharmonic theory and the two empirical models. Therefore, we conclude that the anharmonic approximation is adequate even at very high temperatures of graphite. The QH treatment, on the other hand, appears reasonable up to approximately 1000 K, but after that, it becomes increasingly inaccurate. This is certainly not unexpected; we found similar behaviors in our thermodynamic study of actinide monocarbides and mononitrides [16]. Appl. Sci. 2022, 12, x FOR PEER REVIEW 5 of 9 We calculated the quasi-harmonic free energy in Figure 3 from the GIBBS2 package [34][35][36][37] using the energies determined by FPLMTO calculations of graphite. For all FPLMTO calculations and QH modeling, we made the same approximation to keep the axial c/a ratio constant at an experimental value. Although QH modeling for graphite has been researched, we were not able to find any free-energy data in the literature. Notably, in Figure 3, the CALPHAD and EOS models are in excellent agreement, which was expected considering that both are thermodynamically consistent representations of the experimental data. Furthermore, the figure shows a good agreement between our first-principles anharmonic theory and the two empirical models. Therefore, we conclude that the anharmonic approximation is adequate even at very high temperatures of graphite. The QH treatment, on the other hand, appears reasonable up to approximately 1000 K, but after that, it becomes increasingly inaccurate. This is certainly not unexpected; we found similar behaviors in our thermodynamic study of actinide monocarbides and mononitrides [16].
Next, we focused on the specific heat at constant pressure. Analogous to the free energy, we obtained Cp as the sum of an electronic and a lattice contribution, see Equation (2). Contrary to the actinide systems, where the density of the actinide 5f states is quite large at excitation energies [16], the electronic contribution for graphite is actually very small because there are very few energy states available for electronic excitations; see the discussion of the graphite electronic density of states in Section 4. Although small, Cp el is still included in the total specific heat shown in Figure 4, but on the scale of the figure, the contribution is not clearly visible.
In addition to our anharmonic theory, Figure 4 shows experimental data [8,38,39], EOS models [33,40], results from quasi-harmonic theory [7] and CALPHAD [32]. The QH model [7] assumes the same fixed c/a value as our own modeling. The EOS models attempt to reproduce the most reliable experimental data for carbon. Kerley and Chhabildas [33] represented the zero-temperature EOS with density and temperature as the independent variables and a Birch-Murnaghan equation, while the lattice vibration terms were formulated using conventional Debye and Einstein models. Contrarily, Fried and Howard [40] proposed an EOS that combines a representation where pressure and temperature are the independent variables within an accurate model modified by Murnaghan, which is able to match all known experimental data for carbon. They included thermal effects through the temperature dependence of the thermal expansion coefficients that were compared with the measured data.
With these EOS model approaches [33,40], one can calculate the Gibbs free energy, from which all thermodynamic properties can be derived. Next, we focused on the specific heat at constant pressure. Analogous to the free energy, we obtained C p as the sum of an electronic and a lattice contribution, see Equation (2). Contrary to the actinide systems, where the density of the actinide 5f states is quite large at excitation energies [16], the electronic contribution for graphite is actually very small because there are very few energy states available for electronic excitations; see the discussion of the graphite electronic density of states in Section 4. Although small, C p el is still included in the total specific heat shown in Figure 4, but on the scale of the figure, the contribution is not clearly visible.

Discussion
When comparing our theory with other experiments, such as low-temperature phonon dispersions, Gibbs free energies (CALPHAD or EOS models representing experimental data), or specific heats, as shown in Figures 2-4, we found a rather satisfactory agreement in all cases. Nevertheless, van der Waals interactions, not included in the present density functional theory, play important roles in graphite. As we mentioned, freezing the axial c/a ratio in a hexagonal graphite structure mitigates some of the deficiencies of this theory and abates the need to include van der Waals contributions, which are computationally expensive and complex non-local electron exchange and correlations. Our theory relies on a far simpler LDA assumption for electron interactions. For the smaller 96-atom supercell, however, we studied the effect of the van der Waals interaction by including it in our VASP calculations within the Grimme DFT-D2 method [18]. As it turns out, it had almost no effect on the SCAILD phonon density of states, presumably because constraining both the volume and the c/a axial ratio suppressed the effect of the van der  [40], and CALPHAD (black dashed line) [32] specific heats.
In addition to our anharmonic theory, Figure 4 shows experimental data [8,38,39], EOS models [33,40], results from quasi-harmonic theory [7] and CALPHAD [32]. The QH model [7] assumes the same fixed c/a value as our own modeling. The EOS models attempt to reproduce the most reliable experimental data for carbon. Kerley and Chhabildas [33] represented the zero-temperature EOS with density and temperature as the independent variables and a Birch-Murnaghan equation, while the lattice vibration terms were formulated using conventional Debye and Einstein models. Contrarily, Fried and Howard [40] proposed an EOS that combines a representation where pressure and temperature are the independent variables within an accurate model modified by Murnaghan, which is able to match all known experimental data for carbon. They included thermal effects through the Appl. Sci. 2022, 12, 7556 6 of 8 temperature dependence of the thermal expansion coefficients that were compared with the measured data.
With these EOS model approaches [33,40], one can calculate the Gibbs free energy, from which all thermodynamic properties can be derived.

Discussion
When comparing our theory with other experiments, such as low-temperature phonon dispersions, Gibbs free energies (CALPHAD or EOS models representing experimental data), or specific heats, as shown in Figures 2-4, we found a rather satisfactory agreement in all cases. Nevertheless, van der Waals interactions, not included in the present density functional theory, play important roles in graphite. As we mentioned, freezing the axial c/a ratio in a hexagonal graphite structure mitigates some of the deficiencies of this theory and abates the need to include van der Waals contributions, which are computationally expensive and complex non-local electron exchange and correlations. Our theory relies on a far simpler LDA assumption for electron interactions. For the smaller 96-atom supercell, however, we studied the effect of the van der Waals interaction by including it in our VASP calculations within the Grimme DFT-D2 method [18]. As it turns out, it had almost no effect on the SCAILD phonon density of states, presumably because constraining both the volume and the c/a axial ratio suppressed the effect of the van der Waals interactions. The advantage of the LDA over the DFT-D2 method is computational expediency. Additionally, for this particular study, DFT-D2 has not yet been implemented in the FPLMTO all-electron code that we employ to calculate the electronic contributions.
To better justify our simplification, we analyzed the calculated electronic structure and compared this with X-ray photoemission spectra.
As shown in Figure 5, we contrast our calculated electronic density of states with X-ray photoemission experiments with 122 eV photons [41]. We notice that the density of states at the highest occupied energy level, the Fermi level, is very small. Consequently, the electronic contribution to the specific heat is also small, as mentioned in the previous section.  Gaussian, together with X-ray photoemission spectra with 122 eV photons (red line) [41]. The Fermi level is at zero energy, and it is indicated with a vertical dashed line.
Although not in agreement on a very detailed level with photoemission, the calculated density of states is consistent with the measurements. As discussed in [42], the appearance of the spectra between −20 and −13 eV, i.e., a stepwise increase from zero to a nearly constant value, indicates two-dimensional behavior. An analogous feature is observed also in the calculated density of states. Both measured and calculated spectra show a similarly significant peak located between −5 and −10 eV. The reasonable comparison between measured and modeled electronic structures, together with the acceptable accuracy of the calculated thermodynamic properties discussed above, give us confidence in the quality of our graphite model.
Lastly, we emphasize that constraining the behavior of the lattice vibrations to a quasi-harmonic model above about 1000 K leads to unacceptable inaccuracies in the thermodynamics of graphite. Gaussian, together with X-ray photoemission spectra with 122 eV photons (red line) [41]. The Fermi level is at zero energy, and it is indicated with a vertical dashed line.
Although not in agreement on a very detailed level with photoemission, the calculated density of states is consistent with the measurements. As discussed in [42], the appearance of the spectra between −20 and −13 eV, i.e., a stepwise increase from zero to a nearly constant value, indicates two-dimensional behavior. An analogous feature is observed also in the calculated density of states. Both measured and calculated spectra show a similarly significant peak located between −5 and −10 eV. The reasonable comparison between measured and modeled electronic structures, together with the acceptable accuracy of the