# High-Temperature Thermodynamics Modeling of Graphite
## 2 Computational Approaches  
Because we attempt to model graphite at high temperatures, we must consider the potential importance of anharmonic atomic vibrations. Although very convenient, quasi-harmonic (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\times 5\times 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  
Figure 1: The phase diagram of carbon [4]. 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.