Assessment of the Heat Capacity by Thermodynamic Approach Based on Density Functional Theory Calculations

Abstract

The theoretical aspects of the thermodynamic calculation of the Gibbs energy and heat capacity of a crystalline system within the frame of the Density Functional Theory (DFT) are introduced in the present chapter. Various approximations of phonon motion (harmonic, quasiharmonic, and anharmonic) and their effects on the thermodynamic properties are discussed. The theoretical basis of the thermodynamic approach of the heat capacity of crystals for given thermodynamic conditions is presented, having as example six polymorphs of the magnesium hydrides.

1. Introduction

Calorimetry is the experimental technique that allows the determination of the heat transferred between two systems with different temperatures. Basically, it can be applied to obtain the specific heats of the substances, as well as the heats of phase transitions and formation/decomposition processes. The accuracy of the assay of the abovementioned parameters depends on the accuracy of the mass and temperature measurements and the purity of the investigated samples. Sometimes it is necessary to estimate the specific heats for hypothetic materials or materials that are expensive, difficult to synthesize or dangerous for humans and environment. For such cases the heat capacity can be estimated by thermodynamic calculations using the enthalpies and entropies from the thermodynamic databases for pure materials, as in the CALPHAD method. An alternative way, especially when the databases do not contain usable information, is to use the computing methods that can predict the total energy and the vibrational frequencies of a given system. The density functional theory (DFT) is an electronic structure method that considers the electron correlation at a low computational cost and provides accurate results. The high-quality calculations of the vibration frequencies for periodic and nonperiodic systems in various approximations (harmonic, quasiharmonic, and anharmonic) support the calculation of the thermodynamic properties (enthalpy, entropy, and Gibbs free energy) for different composition, pressure P, and temperature T conditions, which can be used to calculate the heat capacity.

In this chapter, we present the theoretical basis of the thermodynamics approach of the heat capacity for crystals for given thermodynamic conditions P and T based on the DFT calculations, having as example six polymorphs of the magnesium hydrides, which are reported in Ref. [1] where the assessments of the thermodynamic properties and pressure–temperature phase diagram of the magnesium hydride polymorphs are done.

Several polymorphs of magnesium hydride were experimentally identified and validated by several electronic structure calculations [1, 2, 3]. The knowledge of the various MgH₂ phase stability has led to an increase in research regarding the pressure–temperature phase diagram for the magnesium hydride. Under ambient conditions, the magnesium hydride crystallizes as α-MgH₂ phase, with a rutile-type structure (space group P42/mnm) [4]. Bastide et al. found that under high pressure and temperature conditions, the α-MgH2 structure is transformed into β-MgH₂ (space group Pa-3) and γ-MgH₂ (space group Pbcn); by decreasing the pressure, the β-MgH₂ is transformed into γ-MgH₂ [5].

The magnesium hydride P–T phase diagram based on the thermodynamic calculations [1] shows that in the pressure range 0–1.5 GPa, the α-MgH₂ phase is the most stable and that the γ-MgH₂ is stable below 6.2 GPa. Above this value, the ε-MgH₂ becomes the most stable up to 10 GPa, the maximum pressure considered in the study. However, in the region of 5.5–7.5 GPa, with exception of the cubic phase c-MgH₂, all the investigated magnesium hydrides might coexist, since their enthalpies have similar values. The hypothetic c-MgH₂ polymorph is the less stable phase, excepting the interval of 0.0–0.4 GPa, where it has enthalpy values comparable with the ε-MgH₂ phase. The cubic polymorph was identified in an experiment as metastable nanocrystals, which transform to γ- and α-MgH₂ [6]. In a subsequent study [7], we predicted that the formation/decomposition curve over all polymorphs is starting from 591.1 K for the low-pressure P = 0.03 GPa, growing with the applied pressure.

The theoretical framework for the thermodynamic calculations presented in the present chapter of the book can be extended to nonperiodic systems (defects, surfaces, interfaces, alloys, amorphous, fluids, and isolated molecules) still using the periodic formalism, but modeling the system in the supercell method, or finite systems (molecules and macromolecules, clusters), where the molecular orbital formalism (specific to the quantum chemistry, which is proper to the isolated systems) is used instead of the crystal orbital formalism (specific to the quantum solid state). For the finite systems, the contributions of the rotation and translation freedom degrees have to be added to the partition function and the free energy of the system.

2. Total energy computation methods

The computer-assisted simulations of the particle systems require: (i) a model for the system, which specifies the chemical species, the position in space (given in Cartesian, internal, or redundant coordinates) and, in case of dynamic treatments, the velocity of each particle; (ii) the method that describes the interactions between the particles and different parameters regarding the calculation method; (iii) different parameters that describe the simulation method (threshold parameters, calculation schemes, and eventually the parallelization technique); and (iv) the parameters and properties that have to be reported by the simulation software.

An ideal crystal can be represented as an indefinitely extended lattice that can be obtained by the translation of a repeated parallelepipedic box, called unit cell. The unit cell is populated with a set of atoms (called atomic basis) that may be arranged in some special points characterized by a set of symmetry operations that is specific for each polymorph of a crystalline substance. The unit cell is characterized by the lengths and mutual orientations of three lattice vectors that delimit the unit cell shape. The solid-state physics is the science that characterizes and classifies the crystals and tries to establish a relation between the nature of the atoms, the structures formed by them, and the crystal properties. The infinite crystal is reduced to the study of the properties of the unit cell, the smallest piece of crystal that preserves the properties of the entire system. The crystals have some local or extended defects, but without losing their global ordering. The crystalline materials can be built based on the experimental structural data that are collected in several online databases or trying to predict it by molecular mechanics and ab-initio calculations. The symmetry of the periodic systems (1D for polymers, 2D for surfaces or films, and 3D for solids) can be used to reduce the computation effort [8].

The equations of state of a given system are obtain by successive approximations having as starting point the time-dependent or the—independent Schrödinger equation associated to electrons and nuclei that form the system. The relativistic contributions to total energy are important for heavy chemical elements and must be considered at least as a perturbation. Due to the relative light mass of the electron compared with the mass of the nuclei, movement of the nuclei and electrons is separated in the frame of the Born-Oppenheimer approximation. The approximation is suitable when the electrons wave function gradient depending of nuclei positions has very small values. The Born-Oppenheimer approximation is not valid when energy values of different electronic states are very close in energy at some nuclear configurations.

For the finite systems (atoms, molecules, clusters) the mono-electronic wave function ${\phi }_i\left(\overrightarrow{r}\right)={\sum }_{\mu =1}^{\varpi }{c}_{\mathit{\mu i}}{\chi }_{\mu }\left(\overrightarrow{r}\right)$ is developed in analytic or numeric basis sets ${\chi }_{\mu }\left(\overrightarrow{r}\right)$, and thus the complex equations are transformed in some treatable ones. In the Linear Combination of Atomic Orbitals (LCAO) approach [9], the functions ${\chi }_{\mu }\left(\overrightarrow{r}\right)$ are atomic orbitals, and the amplitude of the coefficients $c_{\mathit{\mu i}}$ can be used to interpret the interaction in the system. The radial part of the AO might have different mathematical representations (Slater-type, Gaussian-type, or numeric atom-centered orbitals). In the LCAO approach, the parameters can be decomposed into atomic orbital contributions that can be used to interpret the interaction in the system. The wave functions are called Molecular Orbitals and the corresponding theory, Quantum Chemistry.

For a periodic system (crystal, slabs, surfaces, wires, and tubes), the Born-von Karman boundary condition introduces the expansion of physical quantities to Fourier series. To simplify the solving of the mono-electron Schrödinger equation, the Bloch theorem exploits the translation symmetry of the system and implicitly of the potential, by factorizing the mono-electron wave function as ${\phi }_i\left(\overrightarrow{r}\right)\to {\phi }_{n\overrightarrow{k}}\left(\overrightarrow{r}\right)=e^{i\overrightarrow{k}\overrightarrow{r}}{u}_{n\overrightarrow{k}}\left(\overrightarrow{r}\right)$, where $\overrightarrow{k}$ is a vector defined in the reciprocal space, and n is the band index. Thus, the Bloch theorem indicates the way to reduce the computing of an infinite number of electronic wave functions to a finite number of electron wave functions, as well as the indexation of the electron wave functions by the band index n. The wave functions are called Crystal Orbitals (CO) and the corresponding theory, Quantum Solid-State.

By solving the reduced mono-electronic equation for each $\overrightarrow{k}$, a set of energies ${\in }_{n\overrightarrow{k}}$ is obtained. By using the periodicity of the reciprocal space, a polyhedron called Brillouin zone (BZ) may be defined in the reciprocal space. The wave vectors outside the Brillouin zone simply correspond to states that are physically identical to those states within the BZ. For each band n, the energy levels ${\in }_{n\overrightarrow{k}}$ evolve smoothly with the changes in $\overrightarrow{k}$, forming a continuum band of states. The electronic wave functions for closed $\overrightarrow{k}$-points are very similar, and the integration in the reciprocal space can be reduced to a summation over a grid of k-points [10]. The integrals over k-space converge exponentially with the number of sampling k-points and several recipes are available to compute the sets of spatial $\overrightarrow{k}$-points for different symmetries in order to accelerate the convergence of BZ integrations were developed [11]. Due to the partial filling of the energy bands in the case of the metals, the BZ is discontinuous. In this case, the calculations of the integrals over BZ with a denser$\overrightarrow{k}$-grid and the broadening of the electronic levels may reduce the magnitude of the errors [12].

The development of the electronic functions using a plane wave (PW) basis set $u_{n\overrightarrow{k}}\left(\overrightarrow{r}\right)={\sum }_{\overrightarrow{G}=1}^{\infty }{c}_{\overrightarrow{G}i}\left(n,\overrightarrow{k}\right)e^{i\overrightarrow{G}\overrightarrow{r}}$ is a natural choice for the crystals, as the equations obtained are very similar with those of Nearly Free Electron model [13]. Therefore, the mono-electronic wave function may be written as ${\phi }_{i\overrightarrow{k}}\left(\overrightarrow{r}\right)=e^{i\overrightarrow{k}\overrightarrow{r}}{u}_{n\overrightarrow{k}}\left(\overrightarrow{r}\right)={\sum }_{\overrightarrow{G}=1}^{\infty }{c}_{\overrightarrow{G}i}\left(n,\overrightarrow{k}\right)e^{i\left(\overrightarrow{k}+\overrightarrow{G}\right)\overrightarrow{r}}$. The different terms of the total energy are written as Fourier transforms, thus simplifying the numerical treatments. For the periodic systems, the use of a plane wave basis set in the description of the COs offers a number of advantages, including the simplicity of the basis functions, the absence of basis set superposition error, and the ability to calculate efficiently the forces acting on the atoms. In case of the ionic crystals, the PW method is not efficient because of the high number of plane waves that are required for an accurate description of the wave functions near the ionic core.

Only the electrons that occupy the high energy levels (called valence electrons) are responsible for formation and breaking of the chemical bonds and for the interaction with the low-energy radiation. The rest of the electrons (called core electrons) generally are not affected by the chemical environment and are not as significant as the valence electrons. Thus, treating explicitly only the valence electrons, a further decrease of the computational effort can be achived. The core electrons are emulated by effective core potentials (ECPs) in quantum chemistry LCAO methods or by pseudopotentials (PPs) in the solid-state PW methods [14]. Thus, only the valence electrons are considered in the electronic equations [15]. The core potentials are developed so as to reproduce the energies, as well as the wave function amplitude (outside of a given cutoff radius) of the atomic core wave functions. The PPs that satisfy the condition of the normalization outside of the cutoff radius are called norm-conserving pseudopotentials (NCPP) [16], and those for which this condition is relaxed are called ultrasoft pseudopotentials (USPPs) [17]. The PPs allow one to perform the calculations at a lower energy cutoff. The USPPs are less computationally expensive in comparison with NCPP. For the heavy elements, the relativistic effects can be incorporated in the ECP/PP [18].

2.1 The simulation methods for the electronic ground states

From a broad palette of electronic structure methods, the quantum methods based on the Hartree-Fock Theory (HFT) [19, 20] and the Density Functional Theory (DFT) [20] are the most popular and mostly used. The neglect of the electron correlation in HFT affects the quality of the results, especially the energy-derived properties. The influence of the electron correlation on the band structure and on the cohesion energy in the oxide materials has been the object of a huge number of articles. The electron correlation effects are important for the ionic crystals, as the electron density is localized in a reduced domain around the ionic core. The existing post-Hartee-Fock correction schemes based on the Configuration Interaction, Coupled Cluster Theory, or Møller-Plesset Perturbation Theory are very accurate [21], but too expensive from a computational point of view to be applied to large systems. Some technical difficulties make them very rare in the solid-state software. Fortunately, the methods based on the DFT are good alternatives. In DFT, the total energy is expressed in terms of total electron density, rather than the many-electron wave function specific to HFT. The choice of the exchange-correlation potentials is a matter of trial and error in the DFT, as the method itself does not provide an explicit dependency of the exchange-correlation potentials on the electron density. Developed initially based on the model of the uniform electron distribution in the Local Density Approximation (LDA), the accuracy of the DFT was increased after including additional corrections that consider the variation of the density in the Generalized Gradient Approximation (GGA). Despite the numerous improvements of the DFT methods on the top of GGA (Self-Interaction and Hubbard U corrections, meta-GGA, and hybrid functionals) [22] regarding the overestimation of the band gap for the semiconductors, there are still efforts to find transferable correlation-exchange potentials. The various proposed exchange-correlation potentials do not describe well the weak van der Waals interactions between the two chemical systems. The simplest solution is to treat the dispersion interaction adding an analytic empirical term of London type to the total energy [23].

Due to the similarity of the Hartree-Fock and Kohn-Sham equations, several electronic structure codes work within the framework of both HFT and DFT and can treat the exchange interaction in a hybrid scheme, incorporating the full or partial Fock exchange in the DFT calculations. Both methods can separately treat the two spin orientations up and down of the electrons, in so-called Unrestricted HF in HFT and spin-polarized or spin-density calculations in DFT. This double framework approach allows the study of the electron correlation effects in 0d-3d periodical systems in a unified way.

The several approaches can be applied in order to reduce the calculation effort by approximation of the integrals by simpler formulas or just parameterization as in the Extended Hückel and various Zero Differential Overlap methods [24] in the frame of the HFT or Density Functional Tight Binding (DFTB) in DFT. These methods are called semiempirical methods, as they contain some empirical determined parameters. Generally, the semiempirical methods consider only the valence electrons, the effects of the core electrons being included in the parameterization of the method. The computational effort required by the semiempirical methods is significantly reduced comparing with the ab-initio methods due to the drastically simplification of numeric calculations and of the great reduction of the considered number of electrons. The prices that must be paid are the reduced accuracy and the reduced transferability of the parameters to other chemical systems that those used for parameterization. However, the semiempirical methods can be used as tool for the pre-selection of the materials that can be investigated in the high-throughput computational screening techniques.

In order to further reduce the calculation efforts, the interaction between atoms can be represented by analytical formula developing so-called the empirical force fields (EFFs), where usually the electrons are not explicitly considered. Thus, the electronic freedom degrees are eliminated, and the computing effort is drastically reduced to square of the total number of particles. Further reduction of the number of the freedom degrees can be done by partly or total freezing of the internal geometry of molecules. The parameters of EFF can be obtained by fitting the total energies and forces calculated by ab-initio electronic structure methods [25]. Also, the phonon properties have to be included into the parameterization procedure as reference data for accurate calculations of the thermodynamic properties [26].

2.2 The simulation methods for the electronic excited states

Neither HFT nor DFT is able to treat the electronic excitation and to characterize the electrons in the excited states. The HFT applied to the excited states is equivalent to the electron configuration methods for the ground state. The DFT is constructed based on the ground electronic density and cannot be directly extended to the excited electronic states. The use of multi-reference and the perturbative correlation interaction methods over the Hartree-Fock wave functions is able to characterize the excited states, but they cannot be applied to large systems because of the high computational efforts. Moreover, such methods require very large basis sets, which drastically increase the size of the CPU memory that is necessary to store very large matrices. There are trials to correct for the excited states, the DFT-determined states with the Configuration Interaction [27], Random-Phase Approximation [28], or Machine Learning [29], but such methods are not implemented in the available quantum chemistry or solid-state software.

A more natural approach is to start from the excitation process by the electromagnetic field as an external field. Thus, the HFT and DFT must be reconsidered in the frame of the time-dependent Schrödinger equation. In the case of DFT, an analogous equation to the static Kohn-Sham theorem that states that any expectation value is a functional of the density, and the initial state is established. This formalism is called Time-Dependent Density Functional Theory (TDDFT) [30]. TDDFT has the same problem like the static DFT as the exchange-correlation potential is not defined. The adiabatic approximation treats the exchange-correlation kernel as static, which permits its evaluation from the derivative of the ground state exchange-correlation potential with respect to the density. The simplest choice is the Adiabatic Local Density Approximation, in which the exchange-correlation kernel is calculated from the ground-state LDA functional [31].

The excited electron and the local environment of the electron (hole) behave like a collective excitation called pseudoparticle, which could be treated in the Many-Body Perturbation Theory (MBPT) and the Green’s function formalism. The GW formalism uses an similar equation to Kohn-Sham equation that governs the energy and the wave functions of the quasiparticle, only that the exchange-correlation potential is replaced by an integral over the self-energy operator that incorporates all the electron–electron interactions [32]. The mean field is determined from DFT calculations and is used to calculate the GW interaction terms. When in GW only the energy of the quasiparticle is modified, but the wave function is kept unchanged, then the method is simplified (so-called G0W0), and the energy of the quasiparticle is obtained as a first-order HF or KS energy. The corrected energies in the G0W0 method give accurate band gaps for semiconductors or insulators. The introduction of the self-consistency within GW gives better band gaps than G0W0, but with a much higher computational effort. The use of the hybrid functionals for calculation of the initial wave functions is desired, but the GW calculations are very expensive, and thus, the application is limited to the small systems. There is some improvement of the algorithms that allow application of GW method to large systems. The GW approaches can be applied to the neutral excitations, but are not able to consider the charged excitations. The Bethe-Salpeter equation (BSE), which is based on the two-particle Green’s functions and the effective two-particle interaction kernel, solves the shortcoming [33]. The kernel can be expressed as sum of the derivatives of the Hartree potential function on the self-energy calculation in GW method. The GW and BSE methods predict more accurately the excitation energies and absorption spectra, compared with DFT methods [34]. Unfortunately, the computational effort increases drastically in the order DFT < GW < BSE. For large systems (more than few hundreds of atoms), even the DFT calculations are prohibited. The semiempirical methods with an accurate parameterization [35] can be applied at reduced computational efforts for large systems in the ground [36] or excited states [37].

2.3 Static and dynamic properties

Total energy is dependent on the relative arrangement of the atoms. Changing continuously the position of the atoms, the energy is continuously modified, and this function energy – coordinate is called potential energy surface (PES). The configuration of atoms that are characterized by minimum or maximum values of the system energy corresponds to an equilibrium or transition states, respectively. Such atomic configurations can be determined by using some mathematical methods that modify the position of the atoms in order to minimize or maximize the total energy, preferentially using the first and second derivatives of the total energy. The derivatives are obtained analytically or numerically. The energy and its derivatives calculated by the electronic structure methods can be used as reference data in the empirical force fields parameterization, in order to reproduce the static and dynamics properties of the atomic systems [38]. The static calculations are very useful to characterize the stability, elastic and electronic properties, the vibration spectra and the thermodynamic properties, and the way of transition from a structure to another.

Ehrenfest theorem establishes the theoretical basis for time evolution of a quantum system. Due to nuclei large masses, the Ehrenfest theorem can be reduce to Newton equation, which rules the atoms movement on the PES. The theory that describes the system time evolution is called Molecular Dynamics (MD). The dynamic properties of the investigated system can be accessed through the MD simulations, which consist of the integration of the Newton equation of each atom of the system. The forces that act on the atoms can be evaluated from electronic structure calculation or can be evaluated by much cheaper empirical force fields. Besides total energy minimization, MD simulations allows the simulation: (i) of the behavior of atoms that form the system (ii) of any kind of chemical species using empirical and quantum force field (iii) of the different statistical ensembles. After saving the trajectories in files and data process, we can: (i) visualize the dynamic of the system; (ii) plot and analyze temperature, pressure, distances, volume, tensile forces, and cell parameters; (iii) calculate velocity correlation factor and vibration intensities; (iv) calculate radial distribution functions and structure factors; (v) calculate various thermodynamics parameters in different thermodynamic ensembles [39].

PES structure and system thermodynamics can be characterized by the heuristic methods and data averaging, which describe the system by very large number of attempts. The most common method is Monte Carlo (MC) based on random-walk movement on PES. The sequences of the events are not related with the time evolution of the system as in case of the MD; they are rather dependent on the chosen algorithm. Because the time evolution of the system is not considered, the MC simulation cannot describe the nonequilibrium processes. Time independence can be an advantage for processes that are taking place on a large timescale or in case of PES with a high “roughness.”

2.4 Equation of states and pressure

The calculation of accurate equations of state (EOS), thermodynamic properties, and phase diagrams is a very useful tool in a number of fields, including geophysics [40] and materials research [41]. One of the main advantages of the calculation of pressure and temperature-dependent crystal properties is the ability to investigate the extreme thermodynamic conditions, unattainable by experimental means. Indeed, provided that there are no technical difficulties (e.g., pseudopotential transferability issues [42]), pressure effects can be accounted for simply by compression/expansion of the calculated crystal geometry to smaller/larger volumes compared with the equilibrium volume.

In the present study, the DFT calculations are based on plane waves basis sets techniques combined with ultrasoft pseudopotentials and were carried out with the Quantum Espresso (QE) package [43]. The thermo_pw package [44] was used to obtain the thermodynamic properties within harmonic and quasi-harmonic approximations. The PBE-GGA [45] exchange-correlation potentials were used in the calculations. We have chosen computational settings to ensure that all investigated properties are well converged. The kinetic energy cutoff was set to 55–60 Ry while the charge density cutoff was set to 300 Ry. The integration over the Brillouin zone (BZ) was performed employing a Monkhorst–Pack mesh with a spacing of 0.04 Å⁻¹, as a good balance between the calculation accuracy and the computational effort. The total energy of the polymorphs of the MgH₂ has been minimized function with respect to the unit cell parameters and fractional coordinates of the atoms, preserving the symmetry of the crystal.

The equilibrium volume $V_0$ and the corresponding energy $U_0$ can be determined by a full optimization of the lattice parameters and atom fractional coordinates. Furthermore, the volume dependency of the unit cell $U\left(V\right)$ can be obtained by scaling the lattice constants with the same factor (isotropic procedure) and calculating the corresponding energy U(V). The minimum of the curve corresponds to the equilibrium point ($V_0,U_0$). The curve can be fitted by the Murnaghan Equation of States (EOS) [46], which relates by a simple polynomial equation the volume $V$and its equilibrium value $V_{\mathit{0}}$, the bulk modulus $B=-V{\left(\frac{\partial P}{\partial V}\right)}_T$ and its derivative $B\text{'}={\left(\frac{\partial K}{\partial P}\right)}_T$ by

$$U\left(V\right)=U_{\mathit{0}}+\frac{B_{0}V}{B\text{'}}\left[\frac{{\left(V_0/V\right)}^{B\text{'}}}{B\text{'}-1}+1\right]-\frac{V_0{B}_0}{B\text{'}-1}$$E1

More complex EOSs are available in the literature (see the citations in Refs. [47, 48]), but their qualities are similar. The dependence of the energy on the unit cell volume obtained for the different magnesium hydride structures is presented in Figure 1. The obtained values of the parameters $U_{\mathit{0}}$, $B_0$ and $B\text{'}$ by fitting the EOS given by Eq. (1) are presented in Table 1. The minimum of the curves gives the equilibrium volume V₀ and the equilibrium energy U₀, which are identical with those calculated by a full optimization (lattice parameters and fraction coordinates of the atoms) and are in very good agreement with the experimental values (see Table 1 in Ref. [1]).

Figure 1.

The energy vs. volume for the polymorphs α, β, γ, δ, ε and cubic of the magnesium hydride obtained by an isotropic procedure. The dependency energy-volume determined by an anisotropic procedure is given for the polymorph α-MgH₂.

Polymorph of MgH2	Symmetry group	MgH2 formula unit (f.u.) per unit cell	Total energy U0 (Ry/f.u.)	Electronic smearing term (Ry/f.u.)	Equili-brium volume V0 (Å^3)	Bulk modulus B (GPa)	B´	Debye tempe-rature (K)
α (isotropic)	P42/mnm	2	−143.5961	0.0014	62.50	51.1	3.59	707.4
α (anisotropic)	P42/mnm	2	−143.5961	0.0014	62.48	58.3	4.65	691.2
β	Pa-3	4	−143.5902	0.0002	113.67	54.6	2.48	699.9
c	Fm-3 m	1	−143.5734	−0.0004	27.55	63.7	1.00	511.2
δ	Pbca	8	−143.5880	0.0000	222.23	57.8	1.60	739.9
ε	Pnma	4	−143.5801	0.0000	101.48	61.12	4.10	627.5
γ	Pbcn	4	−143.5966	0.0001	123.15	47.5	3.48	717.7

Table 1.

The properties of the investigated polymorphs of the MgH₂ obtained from the Murnaghan equation of state. For the polymorph α-MgH₂, the values are presented both for the isotropic and anisotropic volume adjustments. The total energy was corrected by cold smearing procedure, included in Quantum espresso code.

At low temperatures, the system energy is only dependent on the system volume, and the pressure can be estimated by

$$P\left(V\right)=-\frac{\mathit{dE}}{\mathit{dV}}=\frac{B_0}{B\text{'}}\left[{\left(\frac{V_0}{V}\right)}^{B\text{'}}-1\right]$$E2

A special care has to be paid when the dependency of energy versus volume is determined in the case of the non-cubic systems, where the lattice parameters and fractional coordinates of the atoms must be optimized for each value of the volume, in order to assure that the energy has the minimum value. Such calculations are done on a grid of lattice parameter under the constraint of a given value of the volume, and those lattice parameters that assure minimum energies are identified. In order to check the effects of the lattice relaxation, we determined the U(V) dependency for the polymorph α-MgH₂ by scaling the lattice constants with the same factor (isotropic procedure) and by a lattice-grid calculation (anisotropic procedure). In the isotropic procedure, nine equidistant scaling factors of the volume between 0.80 and 1.20 were considered. In the anisotropic case, a grid 5 × 5 of a and c/a, centered at the equilibrium values, with steps of 0.05 Å and 0.02, respectively, was considered. The energy was fitted with quartic polynomials as a function of a and c/a and the pairs (a,c/a) for which the energies have a minimum were identified. The results of the Murnaghan fit are given in Table 1 and Figure 1. It can be seen that in case of α-MgH₂, the isotropic and anisotropic U(V) essentially have the same behavior around the equilibrium volume V₀. The same trend also was observed for the other non-cubic polymorphs of the MgH₂.

3. Phonon calculations in lattice dynamics method: Harmonic approximation

The crystal potential energy can be expanded as a Taylor series [49] by small displacements $u_{\alpha ,\mathit{mI}}$ from their equilibrium positions of the atoms I located in the unit cell m, along the Cartesian direction $\alpha =x/y/z$, as

$$U=U_o+U_1+U_2+U_3+\dots$$E3

where:

$U_0$- is the static energy of the system in the equilibrium geometry,

$U_1={\sum }_{\mathit{mI\alpha }}{\mathrm{\Phi }}_{\mathit{mI}}^{\alpha }{u}_{\alpha ,\mathit{mI}}=0$- is zero as the system is in the equilibrium geometry.

The other terms are the n-body crystal potentials (n = 2, 3, …). The second- and the third-order potentials are

$$U_2=\frac{1}{2}\sum _{\mathit{mI\alpha }}\sum _{\mathit{lJ\beta }}{\mathrm{\Phi }}_{\mathit{mI},\mathit{lJ}}^{\mathit{\alpha \beta }}{u}_{\alpha ,\mathit{mI}}{u}_{\beta ,\mathit{lJ}}$$E4

$$U_3=\frac{1}{6}\sum _{\mathit{mI\alpha }}\sum _{\mathit{lJ\beta }}\sum _{\mathit{nK\gamma }}{\mathrm{\Phi }}_{\mathit{mI},\mathit{lJ}}^{\mathit{\alpha \beta \gamma }}{u}_{\alpha ,\mathit{mI}}{u}_{\beta ,\mathit{lJ}}{u}_{\gamma ,\mathit{nK}}$$E5

where ${\mathrm{\Phi }}_{\mathit{mI},\mathit{lJ}}^{\mathit{\alpha \beta }}$ and ${\mathrm{\Phi }}_{\mathit{mI},\mathit{lJ}}^{\mathit{\alpha \beta \gamma }}$ are the harmonic and cubic anharmonic force constants, respectively.

Limiting the expansion to second term, we have the harmonic approximation, which is the fundament of the vibrational frequencies calculations. The reduced Hamiltonian${\mathit{H}}_{\mathrm{HA}}=\frac{1}{2}{\sum }_{\mathit{mI\alpha }}{m}_I{\dot{u}}_{\alpha ,\mathit{mI}}^2+U_2$, which includes the kinetic energy of the atom I with the mass $m_I$, is the harmonic Hamiltonian of the system. The problem of N_a atoms per unit cell that moves in a periodic potential is separable and can be solved exactly by diagonalization of the equation of the eigenvalues and eigenvectors

$$\sum _{\beta ,J}{D}_{\mathit{IJ}}^{\mathit{\alpha \beta }}\left(\overrightarrow{q}\right)e_{\beta ,J}\left(\overrightarrow{q},j\right)={\omega }_{\overrightarrow{q},j}^2{e}_{\alpha ,I}\left(\overrightarrow{q},j\right)$$E6

of the dynamical matrix

$$D_{\mathit{IJ}}^{\mathit{\alpha \beta }}\left(\overrightarrow{q}\right)=\frac{1}{{\left(m_I{m}_J\right)}^{\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}}\sum _l{\mathrm{\Phi }}_{0I,\mathit{lJ}}^{\mathit{\alpha \beta }}{e}^{i\overrightarrow{q}\left({\overrightarrow{r}}_{\mathit{lJ}}-{\overrightarrow{r}}_{0I}\right)}$$E7

The solutions ${\omega }_{\overrightarrow{q},j}$ and $e_{\alpha ,I}\left(\overrightarrow{q},j\right)$ are the frequency and polarization vector that correspond to the phonon normal mode of band index j and wave vector $\overrightarrow{q}$ in the Brillouin zone. The matrix $D_{\mathit{IJ}}^{\mathit{\alpha \beta }}\left(\overrightarrow{q}\right)$ is a hermitic one and its eigenvalues have to be positive. However, sometimes the normal modes might have negative eigenvalues, which correspond to imaginary frequencies. In such a case, the energy decreases along the eigenvector $e_{\alpha ,I}\left(\overrightarrow{q},j\right)$ and the system becomes unstable.

The phonon frequencies for a unit cell in equilibrium (i.e., the energy is minimized and there are no forces on atoms and no stress in the unit cell) can be calculated from the derivative of the energy (phonon freeze method) [50] or considering the forces of each atom in the frame of the linear response method (Density Functional Perturbation Theory—DFPT) [51]. The phonon occupation number at the equilibrium is given by the Bose-Einstein distribution $n_{\overrightarrow{q},j}={\left(exp\left(\mathrm{\hslash }{\omega }_{\overrightarrow{q},j}/k_{B}T\right)-1\right)}^{-1}$. At absolute zero temperature, the phonon population is zero, at low temperatures, there is a small probability for the phonons to exist, and at high temperatures, the number of phonons increases with temperature. The maximum number of phonon modes is given by the maximum number of freedom degrees $3\cdot {\mathrm{N}}_a$, where N_a is the total number of atoms in the unit cell of the crystal.

The relation between ${\omega }_{\overrightarrow{q},j}$ and $\overrightarrow{q}$ for each mode j, namely ${\omega }_j={\omega }_j\left(\overrightarrow{q}\right)$, is called phonon dispersion. The number of the phonon modes with the frequency between $\omega$ and $\omega +\mathit{\Delta \omega }$ gives the phonon density of states $g\left(\omega \right)=1/N{\sum }_{\overrightarrow{q},j}\delta \left(\omega -{\omega }_{\overrightarrow{q},j}\right)$, where N is the number of unit cells in the crystal. The normalization factor $1/\left(3\cdot {\mathrm{N}}_a\right)$ is introduced in order to reduce to 1 the integral of $g\left(\omega \right)$ over frequency $\int g\left(\omega \right)\mathit{d\omega }=1$.

The DFPT method [51], implemented in the module thermo_pw code [44], was used to calculate the phonon spectrum for the MgH₂ polymorphs for each of the volume values on a 6 × 6 × 6 grid of $\overrightarrow{q}$-points and Fourier interpolated in the Brilloun Zone. The phonon density of states (DOSs) are computed by integrating the phonon dispersion in the $\overrightarrow{q}$-space. We present in Figure 2 the phonon band structure and DOS for the polymorph α-MgH₂. The phonon spectra are similar with the one obtained by inelastic neutron scattering [52] and by shell-model EFF lattice dynamics calculations [53]. The phonon spectra for the cubic polymorph c-MgH₂ show a large number of imaginary frequencies for all nine different volumes. Therefore, this polymorph is not considered for the heat capacity calculations.

Figure 2.

The phonon dispersion (left) and density of states (right) computed for the polymorph α-MgH₂.

4. The canonical partition function and the Helmholtz free energy

The canonical partition function of the system is the defined as $Z\left(N,V,T\right)={\sum }_i{\mathrm{e}}^{-\frac{E_i\left(N,V\right)}{k_{B}T}}$, where the summation is done over all the states that are characterized by the energies E_i of the system formed by N particles that occupy the volume V and is kept under the temperature T. k_B is the Boltzmann’s constant.

The thermodynamic parameters of the system can be defined as

$$U=k_B{T}^2{\left(\frac{\partial lnZ}{\partial T}\right)}_{N,V}-\text{average potential energy}$$E8

$$S=k_{B}lnZ+U/T-\text{entropy}$$E9

$$H=U+P\cdot V-\text{enthalpy}$$E10

$$F=U-T\cdot S-\text{Helmholtz free energy}$$E11

In the relation above $k_B$ is the Boltzmann constant, T is the temperature, P is the pressure, and V is the volume of the unit cell. Enthalpy, defined as $H=U+P\cdot V$, is a measure of a system capacity to release heat as a nonmechanical work. The most fundamental thermodynamic parameter is the Gibbs free energy, which is defined as $G\left(P,T\right)=F+p\cdot V=H-T\cdot S$.

Based on the Helmholtz free energy, other important thermodynamic parameters can be defined:

• the bulk modulus $B_T=V{\left(\frac{{\partial }^{2}F}{\partial V^2}\right)}_T$, which describes the resistance of a material to the compression,

• the isochoric heat capacity $C_V=-T{\left(\frac{{\partial }^{2}F}{\partial T^2}\right)}_V$, which is the amount of heat per unit mass of material that is required to raise the temperature by one unit under the condition of the constant volume,

• the isobaric heat capacity $C_P={\left(\frac{\partial H}{\partial T}\right)}_P=C_V-T{\left(\frac{\partial V}{\partial T}\right)}_P^2{\left(\frac{\partial P}{\partial V}\right)}_T$, and the volumetric thermal expansion as ${\alpha }_V=\frac{1}{V}{\left(\frac{\partial V}{\partial T}\right)}_P$. Generally, the energy of the system can be decomposed based on the different degrees of freedom as the electronic and nuclear (translation, rotation, and vibration) motions, domains contributions (bulk, surface, interface, and defects), or phenomena (magnetism, irradiation). Some hybrid contributions have to be added in case of not completely separation of the different freedom degrees (as example, the electron–phonon coupling). In the case of the crystalline systems, the translation and rotation motions of the atomic basis are not present, and their contributions are not included in the structure of the Helmholtz free energy. Furthermore, the product $\mathit{PV}$ is very small in comparison with the other terms, and usually it is neglected. Therefore, usually the stability of the crystals is characterized by the Helmholtz free energy, rather than the Gibbs energy.

For a nonmagnetic ideal crystal, the Helmholtz free energy can be written as:

$$F\left(V,T\right)=U_0\left(V\right)+F_{\mathit{el}}\left(V,T\right)+F_{\mathit{\text{phon}}}\left(V,T\right)$$E12

where the first term in Eq. (12) corresponds to the total energy, the second term is the contribution of the electronic excitation, and the last term is given by the nuclear vibrational motion. The first and third terms can be calculated both by the electronic structure or empirical force fields methods and the second one just by electronic structure methods.

The electronic structure methods typically do not take explicitly into account the temperature effects on the atomic ground state and the computed properties correspond to T = 0 K. However, for the incomplete occupied bands, like in case of the metals, the occupation probability of the electronic levels $\in$ is described by the Fermi-Dirac distribution $f\left(\in ,T\right)={\left(exp\left(\frac{\in -E_F}{k_{B}T}\right)+1\right)}^{-1}$, where the $E_F$ is the Fermi level, the last occupied electronic level, with the value given by the normalization condition $\int f\left(\in ,T\right)d\in =1$. The energy due to the electronic excitation is $U_{\mathit{el}}\left(V,T\right)=\int n\left(V,\in \right)f\in d\in -\int n\left(V,\in \right)\in d\in$, where $n\left(V,\in \right)$ is the electronic density of states computed for the volume V. The corresponding electronic entropy is $S_{\mathit{el}}\left(V\right)=-{\mathit{gk}}_B\int \left[f\left(V,T\right)lnf(V,T)+\left(1-f\left(V,T\right)\right)ln\left(1-f\left(V,T\right)\right)\right]$.

where g is equal to 1 for collinear spin polarized or 2 for non-spin-polarized systems. The electronic Helmholtz free energy is $F_{\mathit{el}}\left(V,T\right)=U_{\mathit{el}}\left(V,T\right)-{\mathit{TS}}_{\mathit{el}}\left(V,T\right)$. The electronic excitation term is calculated based on the Mermin theorem [54]. $F_{\mathit{el}}$is important for the metallic systems, especially at high temperatures. Usually it is neglected in the other cases. However, the partial occupation of the electronic level can be used also for the nonmetallic systems in the smearing scheme [55] that is useful to accelerate the self-consistent field convergence, mixing the occupied and unoccupied states that are clearly separated at T = 0 K. The contribution of the smearing to the electronic energy is refereed as Mermin free energy. The smearing scheme was applied to perform the calculations for the investigated polymorphs of the magnesium hydride. It is efficient in the accelerating the SCF convergence, inducing a very low correction to the total energies (see Table 1). The electronic contribution to the electronic Helmholtz energy is negligible except for the cubic polymorph c-MgH₂, which has metallic properties.

The third term in Eq. (12) is the contribution of the phonon vibration and is dependent on temperature. The partition function corresponding to the phonon is $Z_{\mathit{vib}}=\prod _{\overrightarrow{q},j}\left({\sum }_{n=0}^{\infty }{e}^{-\left(n+\frac{1}{2}\right)\frac{\mathrm{\hslash }{\omega }_{\overrightarrow{q}j}}{k_{B}T}}\right)$ and the vibrational Helmholtz energy in the harmonic approximation is depending on temperature only by the phonon contribution calculated for the equilibrium volume $V_{\mathit{eq}}$

$$F_{\mathit{HA}}\left(V_{\mathit{eq}},T\right)=U_0\left(V\right)+\frac{1}{2}\sum _{\overrightarrow{q},j}\mathrm{\hslash }{\omega }_{\overrightarrow{q},j}+k_{B}T\sum _{\overrightarrow{q},j}log\left(1-e^{-\frac{\mathrm{\hslash }{\omega }_{\overrightarrow{q},j}}{k_{B}T}}\right)$$E13

For low temperatures, the amplitude of the vibrations is reduced and the harmonic approximation is valid for almost all the cases. For $N_a$ atoms per unit cell, the isochoric heat capacity is

$$C_V=-T{\left(\frac{{\partial }^{2}F}{\partial T^2}\right)}_V=k_B{N}_a\sum g\left(\omega \right){\left(\frac{\mathrm{\hslash }\omega }{k_{B}T}\right)}^2\frac{exp\left(\mathrm{\hslash }{\omega }_{\overrightarrow{q},j}/k_{B}T\right)}{exp\left(\mathrm{\hslash }{\omega }_{\overrightarrow{q},j}/k_{B}T\right)-1}$$E14

5. Computational thermodynamics in quasi-harmonic approximation

The neglect of anharmonicity by the truncation of the third term in the development of the total energy leads to well-known unphysical behavior [56]: the zero thermal expansion, infinite thermal conductivity, and phonon lifetime. The inclusion of temperature effects, primarily related to the vibrational degrees of freedom inside the crystal, is more delicate. There are essentially two mainstream ways of incorporating temperature in a theoretical calculation: Molecular Dynamics [57] and Monte Carlo [58] simulations, and the quasiharmonic approximation (QHA) [59]. The former techniques are ideally suited for situations close to the classical limit, at temperatures close to or including the melting temperature. The latter is assuming the harmonic approximation at any given crystal geometry, even if it does not correspond to the equilibrium structure. Plenty of examples of the success of QHA in the prediction of thermodynamic properties and phase stability of solids can be found in the literature [60].

In QHA, the phonon calculations are done for several volumes ${\omega }_{\overrightarrow{q},j}\left(V\right)$ and the Helmholtz free energy becomes

$$F_{\mathit{QHA}}\left(V,T\right)=U_0\left(V\right)+\frac{1}{2}\sum _{\overrightarrow{q},j}\mathrm{\hslash }{\omega }_{\overrightarrow{q},j}\left(V\right)+k_{B}T\sum _{\overrightarrow{q},j}log\left[1-exp\left(-\frac{\mathrm{\hslash }{\omega }_{\overrightarrow{q},j}\left(V\right)}{k_{B}T}\right)\right]$$E15

or by the grouping of the terms [61].

$$F_{\mathit{QHA}}\left(V,T\right)=U_{\mathit{\text{cold}}}\left(V\right)+F_{\mathit{th}}\left(V\right)$$E16

where $U_{\mathit{\text{cold}}}\left(V\right)=U_0\left(V\right)+\frac{1}{2}{\sum }_{\overrightarrow{q},j}\mathrm{\hslash }{\omega }_{\overrightarrow{q},j}\left(V\right)$ is the cold potential energy (T = 0 K), and $F_{\mathit{th}}\left(V\right)=k_{B}T{\sum }_{\overrightarrow{q},j}log\left(1-exp\left(-\frac{\mathrm{\hslash }{\omega }_{\overrightarrow{q},j}\left(V\right)}{k_{B}T}\right)\right)$ is the thermal component of the Helmholtz free energy given by the phonons.

The accurate calculation of the phonon frequency requires the energy evaluation by electronic structure methods, but also the usage of some empirical force fields might give good results [62]. At very low temperature, the thermal part of the phonon contribution is negligible and the entropy does not contribute to the Gibbs energy as the product TS is also negligible. In these conditions, the Gibbs energy of the crystal becomes $G\cong U_{\mathit{\text{cold}}}+\mathit{PV}$. The product $\mathit{PV}$ is very small comparing with the total energy and the Gibbs energy is approximated with the cold energy. The anharmonic effects become significant for the magnesium hydride just above 800 K [63]. Therefore, we do not consider here the anharmonic effects on the thermodynamic properties.

Alternatively to the static calculations of the phonons in the harmonic approximation, the anharmonic phonon spectra can be calculated by Molecular Dynamics simulations [64, 65]. The phonon dispersion of a given phonon wave vector $\overrightarrow{q}$ can be evaluated from MD simulations by computing the Fourier transforms of velocity–velocity correlation functions $g_{\overrightarrow{q}}\left(\omega ,T\right)=\int e^{\mathit{i\omega t}}{\sum }_{I=1}^{N_a}{e}^{i\overrightarrow{q}{\overrightarrow{R}}_I}\frac{{⟨{\overrightarrow{v}}_I\left(t\right){\overrightarrow{v}}_I\left(0\right)⟩}_T}{{⟨{\overrightarrow{v}}_I\left(0\right){\overrightarrow{v}}_I\left(0\right)⟩}_T}\mathit{dt}$, where ${\overrightarrow{v}}_I\left(t\right)$ is velocity of each atom I at time t, and ${\overrightarrow{R}}_I$ is the lattice position of the same atom. The angular brackets represent the time average for a given temperature T considered during the MD simulations. The phonon dispersions are explicitly temperature-dependent and include also the anharmonic contributions. The MD calculations fully consider the anharmonic effects, but do not consider the ZPE effects, as they follow the classical statistical mechanics [66]. Therefore, a quantum correction has to be considered, especially for low temperatures [67].

6. Computational thermodynamics

6.1 Thermodynamic modeling

Considering the temperature effects on the volume the pressure can be formulated as

$$P\left(V,T\right)=-{\left(\frac{\partial F}{\partial V}\right)}_T=P_{\mathit{\text{stat}}}\left(V,0\mathrm{K}\right)+P_{\mathit{\text{phon}}}\left(V,T\right)+P_{\mathit{ae}}\left(V,T\right)$$E17

where $P_0\left(V,0\mathrm{K}\right)$ is the static pressure computed by Eq. (2), $P_{\mathit{\text{phon}}}\left(V,T\right)$ is phonon contribution to the pressure, and $P_{\mathit{ae}}\left(V,T\right)$ is the anharmonic and electronic thermal pressure. The thermal pressure is calculated as the derivative of the thermal free energy

$$P_{\mathit{\text{phon}}}=-{\left(\frac{\partial F_{\mathit{\text{phon}}}}{\partial V}\right)}_T=-\sum _{\overrightarrow{q}j}\left(\frac{1}{2}+\frac{\mathrm{\hslash }{\omega }_{\overrightarrow{q}j}}{e^{\frac{\mathrm{\hslash }{\omega }_{\overrightarrow{q}j}}{k_{B}T}}-1}\right)\mathrm{\hslash }\frac{d{\omega }_{\overrightarrow{q}j}}{\mathit{dV}}=-\sum _{\overrightarrow{q}j}{\gamma }_{\overrightarrow{q}j}\left(\frac{1}{2}+\frac{\mathrm{\hslash }{\omega }_{\overrightarrow{q}j}}{e^{\frac{\mathrm{\hslash }{\omega }_{\overrightarrow{q}j}}{k_{B}T}}-1}\right)$$E18

where ${\gamma }_{\overrightarrow{q},j}=-\frac{V}{{\omega }_{\overrightarrow{q},j}}\frac{\partial ln{\omega }_{\overrightarrow{q},j}}{\partial V}=-\frac{V}{2{\omega }_{\overrightarrow{q},j}^2}⟨e_{\overrightarrow{q},j}\left|\frac{\partial D\left(\overrightarrow{q}\right)}{\partial V}\right|e_{\overrightarrow{q},j}⟩$ is the so-called mode Grüneisen parameter, which characterizes the volume dependences of the frequency of the mode $\left(\overrightarrow{q},j\right)$. Similarly to Eq. (2), the pressure $P\left(V\right)=-{\left(\frac{\mathit{dF}}{\mathit{dV}}\right)}_T$ can be calculated for each considered volume of the unit cell based on the equations of state of Murnaghan type [46]. The enthalpy can be calculated as $H=U_0\left(V\right)+\mathit{PV}$.

6.2 Debye method

The Debye model [68] considers a simple form for the DOS of the vibration modes $g\left(\omega \right)=C\cdot {\omega }^2\cdot \mathrm{\Theta }\left(\omega -{\omega }_D\right)$, where $\mathrm{\Theta }$ is the Heaviside step function, and $C=9N_a/{\omega }_D^3$ is a constant determined from the condition $C\int g\left(\omega \right)\mathit{d\omega }=3N_a$. Thus, the phonon modes are populated just below ${\omega }_D$, which is called Debye frequency. The phonon Helmholtz free energy and the isochoric heat capacity are

$$F_{\mathit{\text{phon}}}^{\mathit{\text{Debye}}}\left(T\right)=N_a{k}_{B}T\left(\frac{9T_D}{8T}+3ln\left(1-e^{-\frac{T_D}{T}}\right)-D\left(\frac{T_D}{T}\right)\right)$$E19

where the function $D\left(\frac{T_D}{T}\right)=\frac{3}{{\left(T_D/T\right)}^3}\underset{0}{\overset{T_D/T}{\int }}\frac{x^3}{e^x-1}\mathit{dx}$ is the Debye integral and the parameter $T_D=\frac{\mathrm{\hslash }}{k_B}{\left[6{\pi }^2{V}^{\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}{N}_a\right]}^{\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$3$}\right.}f\left(\sigma \right)\sqrt{\frac{B_s}{M}}$ is the Debye Temperature, which can be interpreted as the temperature at which each mode below the highest-frequency mode ${\omega }_D$ is excited. The Debye model allows the calculation of the thermodynamic parameters avoiding the phonon calculations if the Debye temperature can be determined from experimental or theoretical elastic constants [69] or from the value of the melting temperature [70]. The estimated Debye temperatures for T = 300 K are given in the Table 1, for all the polymorphs.

6.3 Heat capacity: Grüneisen parameter

The weighted average heat capacity of the individual phonon modes

$$\gamma =\frac{{\sum }_{\overrightarrow{q},j}{\gamma }_{\overrightarrow{q},j}{C}_V\left(\overrightarrow{q},j\right)}{{\sum }_{\overrightarrow{q},j}{C}_V\left(\overrightarrow{q},j\right)}$$E20

is the total Grüneisen parameter, where

$$C_V\left(\overrightarrow{q},j\right)=k_B{\left(\frac{\mathrm{\hslash }{\omega }_{\overrightarrow{q},j}}{k_{B}T}\right)}^2{\left[exp\left(\frac{\mathrm{\hslash }{\omega }_{\overrightarrow{q},j}}{k_{B}T}\right)-1\right]}^{-2}exp\left(\frac{\mathrm{\hslash }{\omega }_{\overrightarrow{q},j}}{k_{B}T}\right)$$E21

is the contribution of each vibration mode $\left(\overrightarrow{q},j\right)$ to the isochoric heat capacity

$$C_V=\sum _{\overrightarrow{q},j}{C}_V\left(\overrightarrow{q},j\right)$$E22

The constant volume (isochoric) heat capacity was obtained from the quasi-harmonic phonon frequencies calculated at each fixed volume. The isochoric heat capacity at each temperature is then obtained interpolating at the temperature-dependent volume values obtained at each temperature from the minimization of the free energy.

The constant pressure (isobaric) heat capacity is obtained as

$$C_P=C_V+\mathit{TV}{\beta }^2{B}_T$$E23

where $B_T\left(V,T\right)=\frac{1}{V}{\left(\frac{{\partial }^{2}F\left(V,T\right)}{\partial V^2}\right)}_T$ is the isothermal bulk modulus, which can be calculated as a function of temperature [61].

In Figure 3, the predicted isochoric and isobaric heat capacities of the magnesium hydride polymorphs are shown to have similar behavior, with deviations at 800 K below 3.81 and 3.44 J mol⁻¹ K⁻¹, respectively. The isochoric and isobaric heat capacities have the same values bellow 400 K. The predicted heat capacities for the polymorph α-MgH₂ are in excellent agreement with the experimental data [52, 71, 72] despite that the calculations are done for ideal crystalline systems, real sample are generally polycrystalline and contain a lot of defects depending on the preparation procedure. The iso- and anisotropic calculated isochoric heat capacities for α-MgH₂ are identical while those for the isobaric heat capacity are almost identical up to 550 K, which is the decomposition temperature of the crystal; even at the higher temperatures, the iso- and aniso- differences are marginal. Similar results regarding the isotropic heat capacities are obtained for the other non-cubic polymorphs. Therefore, we conclude that the isotropic treatment by uniform contraction/expansion of the unit cell of the considered magnesium hydride polymorphs introduces negligible effects on the calculated isochoric heat capacities, below their decomposition temperatures. The electronic contribution to the heat capacity was found to be negligible.

Figure 3.

The isochoric (left) and the isobaric (right) heat capacities versus temperature, computed for the five polymorphs of the magnesium hydride α, β, δ, ε, and γ by the isotropic procedure. In addition, the C_v versus T obtained by the anisotropic procedure is presented for α-MgH₂ (red open circle). The experimental data available for α-MgH₂ are indicated by full symbols.

The heat capacity in the Debye model is $C_V^{\mathit{\text{Debye}}}=3N_a{k}_B\left(4D\left(\frac{T_D}{T}\right)-\frac{3T_D/T}{e^{T_D/T}-1}\right)$. The formula is valid at low temperatures, predicting correctly the temperature dependence of the heat capacity as being proportional to the temperature at power 3, and recovers at high temperatures the Dulong–Petit law, which states that the heat capacity is proportional with the number of atoms per unit cell $3N_a{k}_B$ at high temperatures.

The Molecular Dynamics and Monte Carlo simulations in tandem with the thermodynamic integration [73] and free energy perturbation [74] methods can be used to calculate the Helmholtz free energy and other thermodynamic properties. The electronic structure methods or well-parameterized EFF can be used in MD simulations to the vibration DOS and to calculate the thermodynamic properties such as heat capacity [75]. The heat capacity can be directly evaluated in MD or MC simulations as $C_{V/P}=\frac{⟨E^2⟩-{⟨E⟩}^2}{k_B{T}^2}$, where E is total energy calculated at each simulation step, and the brackets indicate the average over the sequence of configurations produced during the MD or MC steps. Depending on the used thermodynamic ensemble, canonic (N_a, V, T constant) or the isothermal-isobaric (N_a, P, T constant), the isochoric C_v, and isobaric C_p heat capacity, respectively, are calculated. The MD and MC simulations involve the evaluation of the energy and interatomic forces for a very large number of atomic configurations and are not practical to involve electronic structure methods, especially in the case of large systems.

7. Conclusions

The electronic structure and empirical force fields methods are discussed shortly, emphasizing the calculation of the total energy and its derivatives, as well as the phonon spectra. The methods for the calculation of the Helmholtz, Enthalpy, and Gibbs energy are described. The heat capacities of several polymorphs of the magnesium hydride are investigated by density functional theory within the frame of the quasi-harmonic approximation. The predicted temperature variation of the heat capacity for the polymorph of the magnesium hydride α-MgH₂ is in good agreement with the experimental data. We assess the heat capacity for several other polymorphs of MgH₂, for which there is no available experimental or theoretical information.

Acknowledgments

The work was supported by a grant of the Romanian Ministry of Education and Research, CCCDI – UEFISCDI, project number PN-III-P2-2.1-PED-2019-4816, within PNCDI III. CV and AV acknowledge the financial support for mutual visits provided on the basis of the inter-academic exchange agreement between the Bulgarian Academy of Science and the Romanian Academy.

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.