# High-Temperature Thermodynamics Modeling of Graphite
## 2 Computational Approaches
For the Gibbs (or Helmholtz at ambient pressure) free energy, \(F(V,T)\), we formulated the following equation:  
\[F(V,T)=F_{lat}(V,T)+F_{el}(V,T). \tag{1}\]  
\(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:  
\[C_{p}(T)=C_{p}^{lat}(T)+C_{p}^{el}(T). \tag{2}\]  
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 projector-augmented-wave (PAW) method [19; 20; 21]. The computational constraints for VASP are defined by an energy cut-off of 400 eV and a (2 \(\times\) 2 \(\times\) 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 \(\times\) 12 \(\times\) 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 A\({}^{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 A\({}^{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 CALPHAD 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.