Specific heat at constant pressure from first principles: contributions from fully anharmonic vibrations

Specific heat at constant pressure (Cp) for diamond Si and wurtzite GaN is calculated using a novel, first principles method based on density functional theory. This method, termed the Beyond Quasi-Harmonic method, completely takes into account all anharmonic vibrations–both from changes in volume and phonon interactions. Our calculated values for Si show excellent agreement with the generally accepted experimental values. For GaN, however, there is a large discrepancy between two sets of experiments. From an analysis of the available literature, and our own data we develop a set of recommended values.


Introduction
Specific heat at constant pressure (C p ) is a fundamental thermodynamic quantity with applications to protein folding [1,2], geology and earth science [3,4], semiconductors [5], high-temperature ceramics [6], thermal conductivity [7], and nanotechnology [8]. In addition, it is of paramount importance in the calculation of other fundamental thermodynamic quantities including enthalpy and Gibbs free energy.
Currently, however, there is a gap in the theory between the quasi-harmonic approximation, and the UP-TILD [9] 1 and TU-TILD [6] 2 approaches. In the quasi-harmonic approximation, a new set of phonon states is calculated for a range of volumes. The data are fitted with a model assumption in order to obtain the temperature-dependence of the Bulk Modulus, free energy, and coefficient of thermal expansion [10]. The use of empirically-derived equations is also possible [11]. The primary limitations are that the quasi-harmonic approximation only partially takes into account anharmonic effects 3 and semi-empirical methods necessitate performing an experiment. The other main option is the relatively new UP-TILD and the faster TU-TILD method. While quite robust and rich in information these methods suffer from being fairly complicated precisely because they make so few approximations and naturally are somewhat time-consuming to set up and execute.
Thus, there is a gap in the current theory. One can use the quasi-harmonic method, together with model assumptions and/or empirical parameters and limited anharmonic contributions. Or, one can use the TILD methods with virtually no approximations and there's nothing in between. Our method fills this gap. It completely takes into account all anharmonic phonon contributions, does not rely on any experimental results or model assumptions, is no more complicated than the quasi-harmonic approximation and is almost entirely computationally parallel.
The method described in this paper is based on the previously-developed supercell preparation technique [12] and density functional theory (DFT). Taken together, supercell preparation and DFT can be applied to extract the anharmonic vibrational energy for a given temperature. This process can then be used to extract C p (T).
We note here that it is only the phonon contribution that we calculate, and that the electron contribution is beyond the scope of this work. Although, techniques to include the electron contributions are well-established [13,14] and could be used in conjunction with the method presented here.
We calculate C p for two important materials-diamond Si and wurtzite GaN. The experimental values for Si are well-established (see for example figure 5 in [15]), and when we compare our results to the accepted values for Si [5,16] we find excellent agreement.
However, despite considerable experimental work [11,[17][18][19][20][21][22][23][24] there is no clear consensus for what the correct values are for GaN. Using our new method we calculate C p (T) for GaN up to 1500 K and compare our results with the aforementioned experiments. We find that our results show excellent agreement with one set of experiments, while differing from others. With a critical analysis of the literature and our own calculation we determine which GaN values ought to be preferred.
We organize our paper in the following way: section 2 describes our theoretical approach, how we calculate C v and the supercell preparation method. Section 3 describes the anharmonic energy calculations, and the final two sections present and discuss our results.
2. Computational and theoretical approach 2.1. Density functional theory The theoretical results contained here were obtained from self-consistent, first-principles density functional theory [25] with the SIESTA code [26,27]. The exchange-correlation potential is that of Ceperley-Alder [28] as parametrized by Perdew and Zunger [29]. Norm-conserving pseudopotentials in the Kleinman-Bylander [30] form are used to remove the core regions from the calculations. The basis sets for the valence states are linear combinations of numerical atomic orbitals of the Sankey type [31][32][33], generalized to be arbitrarily complete with the inclusion of multiple zeta orbitals and polarization states. In the present calculations, double zeta basis sets are used for the Si and N atoms and double zeta polarized (d-orbitals) for the Ga atoms. These settings have been shown to produce good results for the present materials for a large number of thermodynamic properties [5,26,27]. Both Si and GaN are represented by a 216 atom cell, with a 1×1×1 Monkhorst-Pack [34] k-point sampling, and repeated boundary conditions. Si is in the diamond crystal structure and GaN is in the wurtzite crystal structure. The lattice constants that we found for Si and GaN are 5.446 Åfor Si and a=3.215 Åwith c/a=1.617 for GaN. These are close to the experimental values of 5.431 Å [35] for Si and a=3.189Åand c/a=1.62 [36] for GaN.

Specific heat at constant volume
The matrix elements of the harmonic dynamical matrix are calculated at 0 K using the frozen phonon method, where the matrix is calculated from the small displacements and energy changes of the atoms. The eigenvalues of this matrix give the phonon density of states at the Γ point. In order to get a more complete picture, we sampled an additional 30 k-points in the Brillouin Zone. The density of states (g(ω)) was normalized so that where N is the number of atoms. g(ω) is shown in figure 1. We calculated the specific heat at constant volume using the Helmholtz free energy calculated using the phonon density of states [37]: The results of this calculation are shown in figure 2. Graphs for the low-T values for C v showing the ∼T 3 dependence are shown in figure 3.

Supercell preparation
The supercell preparation method was first developed to calculate the temperature dependence of vibrational mode lifetimes [38], and was later expanded to study heat flow in molecular dynamics simulations [12,39]. Example calculations from previous work include Kapitza resistance [40], isotope effects [41], and effects from other defects including those on the surface [42]. Ergo, it has been previously described in detail and so will only be summarized here-Although a review and a complete description of the method can be found in [12]. Its purpose is to prepare a supercell in a linear combination of normal modes with random mode phases and random mode energies in order to produce a specific microstate which corresponds to a desired initial T distribution. The energies of each mode are determined via an inverse-transform method with a Boltzmann distribution (equation (3)), and the phases are calculated from a simple random number generation, from 0 to 2π: In the above equation γ s is a random number between 0 and 1, E s is the energy of mode s, and b = , where T eq is the desired equilibrium temperature. The normal vibrational modes are calculated using the harmonic approximation. The relative position of each atom for a given microstate can be calculated using the following equations:  Where q s is generalized position coordinate of normal mode s, A s and ω s are the amplitude and frequency of mode s, respectively. r α and a r eq are the final position and equilibrium position of atom α, respectively. The velocity of each atom is obtained by taking the time derivative of r α . Thus, one microstate is completely determined once each mode s is assigned one γ s , and one f s . The harmonic vibrational energy of the resulting microstate (E harm ) is given by Hookʼs law: (3) for E s and plugging that into equation (6) above, we can solve for the amplitude in terms of the random parameter γ s which gives:

Solving equation
The temperature of the material is then calculated from the kinetic energy and confirmed to be what it should be There is one slight difference between the current method and past preparation methods. In the past, the temperature of a distribution (macrostate) would be averaged over a sufficient number of microstates in order to produce the correct result. But here the amplitude of each mode is scaled so that the temperature of each microstate is the same as the mean. This makes the calculation of C p a little bit cleaner. To accomplish this, we use the same method used previously to create a hot slice [12]. Later, we will see that what we create is a hotblock. The amplitude of each mode is multiplied by a scaling factor, so that the microstate has exactly the correct mean temperature. This scaling factor (S) is given by: Where T eq is the desired equilibrium temperature and T 0 is the mode temperature which resulted from equations (3) and (9). Put simply, we create a microcanonical ensemble from a canonical one. It is well-known that the same results will be obtained in either case [44]. From this point forward, we refrain from using the term 'microstate' and instead use 'configuration' to refer to a system in which each atom has been given a specific position and velocity. Of course this scaling is not necessary, but it makes the calculation of the anharmonic energy easier and reduces noise. This is due to the fact that the kinetic energy is no longer calculated from an ensemble average, but rather an exact value (see Results and Discussion section for further details, in particular equations (13) and (15)).

Anharmonic energy
3.1. Negative control Within the supercell preparation method, normal modes are the eigenvectors of the dynamical matrix in the harmonic approximation mentioned previously (see section 2.2). As such, the system is given an initial energy and temperature while assuming only the harmonic approximation, which at low temperatures is very often quite a good approximation (see figure 2 of [12]). For higher temperatures however it is not, and since DFT inherently takes all anharmonic terms into account, a configuration with energy in excess of what would be expected from the harmonic energy equation (equation (6)) will frequently be constructed. In the past these effects were seen as a source of error [12], but here they are seen as a useful phenomenon from which we can extract information.
We have observed a significant amount of data in a variety of systems and materials in which a high temperature will often produce a system with excess energy [39,40], but also observed that with low temperatures this phenomenon becomes much more rare (see again figure 2 [12]). Because these errors only occurred at high temperatures we strongly suspected that they were caused by anharmonic effects. To confirm this was actually true and not simply an error in our calculations that just so happened to occur only at high temperatures, we designed a negative control in which the temperature would be quite high, but the anharmonic effects would be minimal. It is well-known that changes in volume are the cause of a majority of anharmonic effects in many materials [45]. So this control came in the form of a system in which the volume was kept constant, but the temperature was near-melting. We set up 50 configurations of a 216-atom Si supercell with a temperature of 1450 K using the supercell preparation method. 1450 K is close to the melting point of Si at 1685 K [35]. At temperatures comparable to this (and even much lower), there were difficulties with constructed configurations which were simply too hot [39,40]. The volume was held strictly constant via repeated boundary conditions. If the effects that we saw previously were indeed caused by anharmonic effects, then those excess energy 'errors' should disappear. The result was a system with an average temperature of 1443 +/− 15 K (data not shown). This is in good agreement (within 0.5%) with the 1450 K that would be expected from a ideal, perfectly harmonic system.
Thus, the supercell preparation technique works perfectly well in any situation in which anharmonic effects are low-constant volume, or low temperature [12]-only producing excess energy in situations where anharmonic effects are expected to be significant-high temperature. And since the supercell preparation method assumes all anharmonic contributions to be zero, but is run in a DFT-environment in which all anharmonic terms are taken into account; we can be fairly certain that the excess energy seen in our DFT calculations is not an error but must be the anharmonic vibrational energy.

Calculating anharmonic energy
The anharmonic energy is calculated for each material at 50 K to 100 K temperature intervals. In order to allow for volume expansion while also using repeated boundary conditions we create a hot block inside our supercell using the supercell preparation method, which is free to expand into its surroundings. For Si this took the form of a near-perfectly cubic block consisting of 28 atoms. For GaN, this was a hexagonal tube consisting of 24 atoms, 12 Ga and 12 N. A cubic block was used with our cubic crystal system, and a hexagonal block with the hexagonal system, in order to ensure that each of the normal vibrational modes would be proportionally represented within the hot block. Each T distribution-approximated to be a set of 100 configurations-is constructed in such a way so as to ensure that the temperature the system would have in equilibrium is small so that we can be sure that volume changes of the entire supercell will be small. We chose 200 K, but we note here that this is the temperature calculated from the harmonic approximation so the actual final temperature is usually higher, but does not go above 315 K for any of the systems discussed here. From previous C v calculations and C p experiments, this seemed to be a reasonable compromise between minimizing the effects volume changes of the supercell while remaining close to standard conditions [5]. This means that the difference between each T distribution is not the final temperature or energy of the supercell, but rather the initial T of the hot block. So, with each successive increase of 50 K (figure 4), the hot block increases 50 K and the remaining part of the supercell has its temperature reduced appropriately so as to make the final temperature of the entire supercell 200 K. The anharmonic energy is given by whatever is left over after we subtract off the harmonic energy from the total energy, which is given by Where E total refers to the total vibrational energy and the potential energy is taken to be equal to be the Kohn-Sham energy with the energy of the relaxed coordinates subtracted off, represented as E Kohn-Sham 4 [46][47][48]. Equivalently, the total energy can be written in terms of the sum of the harmonic (E harm ) and anharmonic (E anharm ) energies: Setting the right-hand-side of equation (11) equal to the right-hand-side of equation (10)

( )
Where on the right hand side of equation (13) E harm has been replaced by equation (6). Here we see why it was convenient to scale the amplitudes. Now the anharmonic energy is given by the portion of the Kohn-Sham energy which is in excess of what would be expected from the harmonic approximation, since the kinetic energy was scaled to be what one would expect from the harmonic approximation.

Results
Once the anharmonic energy is obtained, it is plotted against the temperature (figure 4) and then a polynomial fit and its temperature derivative are calculated and added as a premium to C v (figure 5).
C p (T) is given by the following equation: where C v (T) is calculated using the methods described in section 2. Combining the two above equations we find ( ) It should be noted that in equation (15) that we equate the derivatives of E total and H (enthalpy) but we do not equate the energies themselves. One primary reason for this is that DFT does not properly take into account the zero-point energy and so these two quantities differ by a constant. The anharmonic energy (equation (18)) was not added to C v for any temperatures below 100 K. For theses temperatures it was assumed that C v ≈C p is a good approximation. Furthermore, equation (18) does not include any points below 100 K, and so would not be valid for these temperatures.

Discussion
We limit our discussion to temperatures above 300 K, since the values below this temperature are well established for many materials, including the GaN [19][20][21] and Si [5,16] used here. Furthermore, the theoretical calculations are not difficult because » C C v p below room temperature. The results of our calculations are shown in figure 5. The individual points are experiment with the solid and dotted lines being our calculations. Our values for Si average within 1% of the generally accepted values given in [16]. Of the values available in the literature for GaN there seems to be two main groups Leitner [23], Itagaki [24] and Jacob [17] on one side with Sanati [11] and Zieborak [18] on the other. In [17], Jacob developed a set of recommended values while analyzing the data sets from the first group and these are the data plotted in figure 5 for GaN. Therefore, for simplicity, the first group's data can be thought of as represented-at least qualitativelyby the hollow black squares in figure 5. The second groups' experiments are plotted individually. The data Figure 5. The specific heat C p as calculated in this work is given by the solid red lines and C v the dotted blue line in both the top and bottom graphs. Top: 216 atom diamond Si supercell. The hollow black squares are experimental data taken from [16]. Bottom: 216 atom wurtzite GaN supercell. The hollow black squares are taken from [17], the hollow black circles are taken from [18], the black dashed line is the quasi-harmonic calculations taken from [11] and the hollow green triangles are the experimental values, also taken from [11].
obtained by [22] and [19] differ so greatly from all other experiments that they won't be considered going forward.
Our values for GaN disagree with the first group [17,23,24], but agree quite well with the second [11,18], averaging within 1%. We feel that the second group is likely the more accurate for a number of reasons. First, Sanati and Zieborak are the only two experiments where the measurements are performed directly on solid crystalline GaN which is what our calculations were based on. The second group of experiments used polycrystalline and/or powdered samples [23,24] or indirect calculations [17]. Second, we would note that the data in [11] agrees very well with the quasi harmonic calculation from the same publication. This calculation was semi-empirical, drawing from data obtained from other experiments. The quasi-harmonic curve that they obtain is slightly below the data, as is expected from a quasi-harmonic calculation [49]. Lastly, there is some evidence that using polycrystalline GaN samples will cause sublimation at high temperatures and artificially increase the specific heat [19]. This could explain the effective split that we see in the data. But we don't want to over-emphasize this because the results of that experiment are not being considered here, as they are considerably different from all other experiments. Nonetheless, we feel it is worth pointing out that in this regard Lee's data agrees qualitatively with our own and the rest of the literature. For these reasons, we believe that our results along with those of Sanati and Zieborak give the most reliable values.

Limitations
The final order of the polynomial fit from figure 4 was determined when successive terms failed to change the values of C p appreciably. For both Si and GaN this was the 3rd-order term, since adding the 4th-order term changes very little (see figure 6).
The greatest amount of error for our calculations will always occur at the highest-T endpoints-1500 K for both materials. These points are technically undefined since they do not have an adjacent point from which a derivative can be calculated, and so they are formally an extrapolation.
The anharmonic energy, as calculated, is technically a lower bound. Pinning down exactly which atoms can be rightly said to have anharmonic energy is extremely difficult. This is partly because the nature of DFT itself-it has only one Slater determinate, and the Kohn-Sham energy is a function of the electron density only. But it is also a result of very common approximations within DFT, like the Local Density Approximation (used here) which assumes that the exchange-correlation potential is the same as that of a homogeneous electron gas. In addition the volume expansion that the hot block experiences presses on the rest of the supercell so some of the energy is necessarily transferred.
This means that the anharmonic energy ( figure 4) is assumed to be evenly distributed over the entire supercell, and not just in the hot block. This means that the anharmonic energy is technically a lower-bound, but as long as a sufficiently large portion of the supercell is chosen the values are very close to correct. We use roughly 1/8 of the whole system as the hot block-28 and 24 atoms for Si and GaN, respectively.
If a block which is too large is chosen, then a phenomenon similar to the negative control that was mentioned previously in section 3.1 will be obtained with the results for C p tending to an artificially small value. This scenario is depicted in figure 7. We used a 30-atom hot block, with a constant equilibrium temperature of Figure 6. The results for C p are shown when a fit of the anharmonic energy is taken to the 4th-order (green dashed line) are compared to the 3rd-order fit (solid, red line), and to the 2nd-order fit (dash-dotted black line). The 3rd-order fit was used in the paper. The data are for the Si system. 300 K. The anharmonic energy largely stopped increasing once the hot block reached 600 K. Interestingly, using just the three points that we were able to obtain (300 K, 400 K, and 500 K) and extrapolating out to 1500 K and back to 200 K yields practically the same values for C p as the full calculations with a hot block of 24 atoms. The data for this 30-atom hot block is shown in figure 7. The data for the 30 atom block only has 3 points and the y-intercept of the fit was set to zero, so the 24 atom curve shown in figure 5 should be considered the more reliable of the two. But we note here that the results are nearly identical for the temperature ranges where both should be considered reliable (through 600 K), and are easily within 1% for higher temperatures.

Conclusions
A new, completely first-principles method for calculating C p has been developed and shown to produce reliable and computationally efficient values. The presented method is computationally efficient with the only nonparallel element being the single self-consistent DFT step needed to calculate the Kohn-Sham energy.
This method can be used on its own as it was here, to study materials with moderate to wide band gaps, or it can be used to tandem with more complete methods [14] if electron contributions are needed.
In addition, we find that the heat capacity values given in [11,18] for GaN ought to be preferred over other experiments. However, we do not claim to be the final word on this, and recognize that other interpretations of the data are conceivable. More experiments need to be performed on solid wurtzite GaN samples in order for our preference to become certainty, but until then the balance of evidence is in favor of the data given in [11,18].

ORCID iDs
Christopher M Stanley https:/ /orcid.org/0000-0002-0595-3877 Figure 7. Shown above are the data for GaN when a 30-atom block was used with a 300 K equilibrium temperature. Top: the hollow red circles are the data where anharmonic energy is allowed to grow uninhibited. The hollow blue square represents the point where the anharmonic energy was inhibited, and would have given an artificially low value for C p . The dotted black line is a fit to the data, which is used in the calculation shown in the bottom graph. Bottom: the black dashed line is C p calculated using the fit in the top graph, which used a 30 atom block, T eq =300 K. The solid red line is the same line used in figure 5, a 24 atom block, T eq =200 K. C v is included for reference as the dotted blue line.