Open access peer-reviewed chapter

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

Written By

Viorel Chihaia, Valentin Alexiev and Hasan S. AlMatrouk

Submitted: 22 December 2021 Reviewed: 28 February 2022 Published: 18 April 2022

DOI: 10.5772/intechopen.104083

From the Edited Volume

Applications of Calorimetry

Edited by José Luis Rivera-Armenta and Cynthia Graciela Flores-Hernández

Chapter metrics overview

209 Chapter Downloads

View Full Metrics


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.


  • density functional theory
  • phonons
  • thermodynamic calculations
  • Gibbs free energy
  • heat capacity
  • 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 MgH2 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 α-MgH2 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 β-MgH2 (space group Pa-3) and γ-MgH2 (space group Pbcn); by decreasing the pressure, the β-MgH2 is transformed into γ-MgH2 [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 α-MgH2 phase is the most stable and that the γ-MgH2 is stable below 6.2 GPa. Above this value, the ε-MgH2 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-MgH2, all the investigated magnesium hydrides might coexist, since their enthalpies have similar values. The hypothetic c-MgH2 polymorph is the less stable phase, excepting the interval of 0.0–0.4 GPa, where it has enthalpy values comparable with the ε-MgH2 phase. The cubic polymorph was identified in an experiment as metastable nanocrystals, which transform to γ- and α-MgH2 [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 φir=μ=1ϖcμiχμr is developed in analytic or numeric basis sets χμr, and thus the complex equations are transformed in some treatable ones. In the Linear Combination of Atomic Orbitals (LCAO) approach [9], the functions χμr are atomic orbitals, and the amplitude of the coefficients cμ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 φirφnkr=eikrunkr, where 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 k, a set of energies nk 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 nk evolve smoothly with the changes in k, forming a continuum band of states. The electronic wave functions for closed 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 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 denserk-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 unkr=G=1cGinkeiGr 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 φikr=eikrunkr=G=1cGinkeik+Gr. 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 Å−1, as a good balance between the calculation accuracy and the computational effort. The total energy of the polymorphs of the MgH2 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 V0 and the corresponding energy U0 can be determined by a full optimization of the lattice parameters and atom fractional coordinates. Furthermore, the volume dependency of the unit cell UV 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 (V0,U0). The curve can be fitted by the Murnaghan Equation of States (EOS) [46], which relates by a simple polynomial equation the volume Vand its equilibrium value V0, the bulk modulus B=VPVT and its derivative B'=KPT by


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 U0, B0 and B' by fitting the EOS given by Eq. (1) are presented in Table 1. The minimum of the curves gives the equilibrium volume V0 and the equilibrium energy U0, 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 α-MgH2.

Polymorph of MgH2Symmetry groupMgH2 formula unit (f.u.) per unit cellTotal energy U0 (Ry/f.u.)Electronic smearing term (Ry/f.u.)Equili-brium volume V03)Bulk modulus B (GPa)Debye tempe-rature (K)
α (isotropic)P42/mnm2−143.59610.001462.5051.13.59707.4
α (anisotropic)P42/mnm2−143.59610.001462.4858.34.65691.2
cFm-3 m1−143.5734−0.000427.5563.71.00511.2

Table 1.

The properties of the investigated polymorphs of the MgH2 obtained from the Murnaghan equation of state. For the polymorph α-MgH2, 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


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 α-MgH2 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 α-MgH2, the isotropic and anisotropic U(V) essentially have the same behavior around the equilibrium volume V0. The same trend also was observed for the other non-cubic polymorphs of the MgH2.


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α,mI from their equilibrium positions of the atoms I located in the unit cell m, along the Cartesian direction α=x/y/z, as



U0- is the static energy of the system in the equilibrium geometry,

U1=mIαΦmIαuα,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


where ΦmI,lJαβ and ΦmI,lJαβγ 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 HamiltonianHHA=12mIαmIu̇α,mI2+U2, which includes the kinetic energy of the atom I with the mass mI, is the harmonic Hamiltonian of the system. The problem of Na 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


of the dynamical matrix


The solutions ωq,j and eα,Iqj are the frequency and polarization vector that correspond to the phonon normal mode of band index j and wave vector q in the Brillouin zone. The matrix DIJαβq 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α,Iqj 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 nq,j=expωq,j/kBT11. 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 3Na, where Na is the total number of atoms in the unit cell of the crystal.

The relation between ωq,j and q for each mode j, namely ωj=ωjq, is called phonon dispersion. The number of the phonon modes with the frequency between ω and ω+Δω gives the phonon density of states gω=1/Nq,jδωωq,j, where N is the number of unit cells in the crystal. The normalization factor 1/3Na is introduced in order to reduce to 1 the integral of gω over frequency gω=1.

The DFPT method [51], implemented in the module thermo_pw code [44], was used to calculate the phonon spectrum for the MgH2 polymorphs for each of the volume values on a 6 × 6 × 6 grid of q-points and Fourier interpolated in the Brilloun Zone. The phonon density of states (DOSs) are computed by integrating the phonon dispersion in the q-space. We present in Figure 2 the phonon band structure and DOS for the polymorph α-MgH2. 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-MgH2 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 α-MgH2.


4. The canonical partition function and the Helmholtz free energy

The canonical partition function of the system is the defined as ZNVT=ieEiNVkBT, where the summation is done over all the states that are characterized by the energies Ei of the system formed by N particles that occupy the volume V and is kept under the temperature T. kB is the Boltzmann’s constant.

The thermodynamic parameters of the system can be defined as

U=kBT2lnZTN,Vaverage potential energyE8
F=UTSHelmholtz free energyE11

In the relation above kB 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+PV, 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 GPT=F+pV=HTS.

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

  • the bulk modulus BT=V2FV2T, which describes the resistance of a material to the compression,

  • the isochoric heat capacity CV=T2FT2V, 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 CP=HTP=CVTVTP2PVT, and the volumetric thermal expansion as αV=1VVTP. 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 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:


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 is described by the Fermi-Dirac distribution fT=expEFkBT+11, where the EF is the Fermi level, the last occupied electronic level, with the value given by the normalization condition fTd=1. The energy due to the electronic excitation is UelVT=nVfdnVd, where nV is the electronic density of states computed for the volume V. The corresponding electronic entropy is SelV=gkBfVTlnf(VT)+1fVTln1fVT.

where g is equal to 1 for collinear spin polarized or 2 for non-spin-polarized systems. The electronic Helmholtz free energy is FelVT=UelVTTSelVT. The electronic excitation term is calculated based on the Mermin theorem [54]. Felis 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-MgH2, 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 Zvib=q,jn=0en+12ωqjkBT and the vibrational Helmholtz energy in the harmonic approximation is depending on temperature only by the phonon contribution calculated for the equilibrium volume Veq


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


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 ωq,jV and the Helmholtz free energy becomes


or by the grouping of the terms [61].


where UcoldV=U0V+12q,jωq,jV is the cold potential energy (T = 0 K), and FthV=kBTq,jlog1expωq,jVkBT 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 GUcold+PV. The product 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 q can be evaluated from MD simulations by computing the Fourier transforms of velocity–velocity correlation functions gqωT=eiωtI=1NaeiqRIvItvI0TvI0vI0Tdt, where vIt is velocity of each atom I at time t, and RI 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


where P0V0K is the static pressure computed by Eq. (2), PphonVT is phonon contribution to the pressure, and PaeVT is the anharmonic and electronic thermal pressure. The thermal pressure is calculated as the derivative of the thermal free energy


where γq,j=Vωq,jlnωq,jV=V2ωq,j2eq,jDqVeq,j is the so-called mode Grüneisen parameter, which characterizes the volume dependences of the frequency of the mode qj. Similarly to Eq. (2), the pressure PV=dFdVT 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=U0V+PV.

6.2 Debye method

The Debye model [68] considers a simple form for the DOS of the vibration modes gω=Cω2ΘωωD, where Θ is the Heaviside step function, and C=9Na/ωD3 is a constant determined from the condition Cgω=3Na. Thus, the phonon modes are populated just below ωD, which is called Debye frequency. The phonon Helmholtz free energy and the isochoric heat capacity are


where the function DTDT=3TD/T30TD/Tx3ex1dx is the Debye integral and the parameter TD=kB6π2V12Na13fσBsM is the Debye Temperature, which can be interpreted as the temperature at which each mode below the highest-frequency mode ω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


is the total Grüneisen parameter, where


is the contribution of each vibration mode qj to the isochoric heat capacity


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


where BTVT=1V2FVTV2T 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−1 K−1, respectively. The isochoric and isobaric heat capacities have the same values bellow 400 K. The predicted heat capacities for the polymorph α-MgH2 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 α-MgH2 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 Cv versus T obtained by the anisotropic procedure is presented for α-MgH2 (red open circle). The experimental data available for α-MgH2 are indicated by full symbols.

The heat capacity in the Debye model is CVDebye=3NakB4DTDT3TD/TeTD/T1. 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 3NakB 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 CV/P=E2E2kBT2, 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 (Na, V, T constant) or the isothermal-isobaric (Na, P, T constant), the isochoric Cv, and isobaric Cp 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 α-MgH2 is in good agreement with the experimental data. We assess the heat capacity for several other polymorphs of MgH2, for which there is no available experimental or theoretical information.



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.


  1. 1. AlMatrouk HS, Chihaia V, Alexiev V. Density functional study of the thermodynamic properties and phase diagram of the magnesium hydride. Calphad. 2018;60:7-15. DOI: 10.1016/j.calphad.2017.11.001
  2. 2. Vajeeston P, Ravindran P, Hauback BC, Fjellvåg H, Kjekshus A, Furuseth S, et al. Structural stability and pressure-induced phase transitions in MgH2. Physical Review B. 2006;73:224102-224108. DOI: 10.1103/PhysRevB.73.224102
  3. 3. Yartys VA et al. Magnesium based materials for hydrogen based energy storage: Past, present and future. International Journal of Hydrogen Energy. 2019;44:7809-7859. DOI: 10.1016/j.ijhydene.2018.12.212
  4. 4. Hakim K, Rivoldini A, Van Hoolst T, Cottenier S, Jaeken J, Chust T, et al. A new ab initio equation of state of hcp-Fe and its implication on the interior structure and mass-radius relations of rocky super-earths. Icarus. 2018;313:61-78. DOI: 10.1016/j.icarus.2018.05.005
  5. 5. Bastide JP, Bonnetot B, Letoffe JM, Claudy P. Polymorphisme de l'hydrure de magnesium sous haute pression. Materials Research Bulletin. 1980;15:1215-1224. DOI: 10.1016/0025-5408(80)90024-0 ibid. 1980;15:1779-1787. DOI: 10.1016/0025-5408(80)90197-X
  6. 6. El-Eskandarany MS, Banyan M, Al-Ajmia F. Discovering a new MgH2 metastable phase. RSC Advances. 2018;8:32003-32008. DOI: 10.1039/C8RA07068G
  7. 7. AlMatrouk HS, Al-Ajmi F, Do NS, Chihaia V, Alexiev V. The pressure-temperature phase diagram assessment for magnesium hydride formation/decomposition based on DFT and CALPHAD calculations. Modern Approaches on Material Science (MAMS). 2021;4:467-478. DOI: 10.32474/MAMS.2021.04.000180
  8. 8. Orlando R, De La Pierre M, Zicovich-Wilson CM, Erba A, Dovesi R. On the full exploitation of symmetry in periodic (as well as molecular) self-consistent-field ab initio calculations. The Journal of Chemical Physics. 2014;141:104108-104109. DOI: 10.1063/1.4895113
  9. 9. Astier M, Pottier N, Bourgoin JC. Linear-combination-of-atomic-orbitals, self-consistent-field method for the determination of the electronic structure of deep levels in semiconductors. Physical Review B. 1979;19:5265-5276. DOI: 10.1103/PhysRevB.19.5265
  10. 10. Folland NO. Finite-sum approximations to cubic Brillouin-zone integrals. Physical Review B. 1980;22:3669-3677. DOI: 10.1103/PhysRevB.22.3669
  11. 11. Monkhorst HJ, Pack JD. Special points for Brillouin-zone integrations. Physical Review B. 1976;13:5188-5192. DOI: 10.1103/PhysRevB.13.5188
  12. 12. Methfessel M, Paxton AT. High-precision sampling for Brillouin-zone integration in metals. Physical Review B. 1989;40:3616-3621. DOI: 10.1103/PhysRevB.40.3616
  13. 13. Francis GP, Payne MC. Finite basis set corrections to total energy pseudopotential calculations. Journal of Physics: Condensed Matter. 1990;2:4395-4404. DOI: 10.1088/0953-8984/2/19/007
  14. 14. Schwerdtfeger P. The pseudopotential approximation in electronic structure theory. ChemPhysChem. 2011;12:3143-3155. DOI: 10.1002/cphc.201100387
  15. 15. Kresse G, Hafner J. Norm-conserving and ultrasoft pseudopotentials for first-row and transition elements. Journal of Physics: Condensed Matter. 1994;6:8245-8257. DOI: 10.1088/0953-8984/6/40/015
  16. 16. Bachelet GB, Hamann DR, Schlüter M. Pseudopotentials that work: From H to Pu. Physical Review B. 1982;26:4199-4228. DOI: 10.1103/PhysRevB.26.4199. Erratum. ibidem, 1984;29:2309-2309. DOI: 10.1103/PhysRevB.29.2309
  17. 17. Vanderbilt D. Soft self-consistent pseudopotentials in a generalized eigenvalue formalism. Physical Review B, 1990;41:7892-7895(R). DOI: 10.1103/PhysRevB.41.7892
  18. 18. Elsasser C, Takeuchi N, Ho KM, Chan CT, Braun P, Fahnle M. Relativistic effects on ground state properties of 4d and 5d transition metals. Journal of Physics: Condensed Matter. 1990;2:4371-4394. DOI: 10.1088/0953-8984/2/19/006
  19. 19. Pisani C, Dovesi R, Roetti C. In: Berthier G et al., editors. Hartree-Fock Ab Initio Treatment of Crystalline Systems. Part of Lecture Notes in Chemistry. Vol. 48. Berlin Heidelberg: Springer-Verlag; 1988. DOI: 10.1007/978-3-642-93385-1 ISBN-13: 978-3-540-19317-3
  20. 20. Kantorovich L. Quantum Theory of the Solid State: An Introduction. Fundamental Theories of Physics vol. 136. A Van Der Merwe. 2004 Springer Science+Business Media Dordrecht. DOI 10.1007/978-1-4020-2154-1. ISBN 978-1-4020-2153-4
  21. 21. Townsenda J, Kirklanda JK, Vogiatzis KD. Post-Hartree-Fock methods: configuration interaction, many-body perturbation theory, coupled-cluster theory. In: Blinder SM, House JE, editors. Mathematical Physics in Theoretical Chemistry. Amsterdam, Netherlands: Elsevier; 2019. pp. 63-117. ISBN: 978-0-12-813651-5
  22. 22. Csonka GI, Perdew JP, Ruzsinszky A, Philipsen PHT, Lebègue S, Paier J, et al. Assessing the performance of recent density functionals for bulk solids. Physical Review B. 2009;79:155107. DOI: 10.1103/PhysRevB.79.155107
  23. 23. Steinmann SN, Corminboeuf C. Comprehensive benchmarking of a density-dependent dispersion correction. Journal of Chemical Theory and Computation. 2011;7:3567-3577. DOI: 10.1021/ct200602x
  24. 24. Parr RG. A Method for estimating electronic repulsion integrals over LCAO MOs in complex unsaturated molecules. The Journal of Chemical Physics. 1952;20:1499-1499. DOI: 10.1063/1.1700802
  25. 25. Gale JD, Wright K. Lattice dynamics from force-fields as a technique for mineral physics. Reviews in Mineralogy and Geochemistry. 2010;71:391-411. DOI: 10.2138/rmg.2010.71.18
  26. 26. Rohskopf A, Seyf HR, Gordiz K, Tadano T, Henry A. Empirical interatomic potentials optimized for phonon properties. npj Computational Materials. 2017;3:27-27. DOI: 10.1038/s41524-017-0026-y
  27. 27. Marian CM, Heil A, Kleinschmidt M. The DFT/MRCI method. WIREs Computational Molecular Science. 2019;9(e1394):1-31. DOI: 10.1002/wcms.1394
  28. 28. Ziegler T, Krykunov M, Autschbach J. Derivation of the RPA (random phase approximation) equation of ATDDFT (adiabatic time dependent density functional ground state response theory) from an excited state variational approach based on the ground state functional. Journal of Chemical Theory and Computation. 2014;10:3980-3986. DOI: 10.1021/ct500385a
  29. 29. Westermayr J, Marquetand P. Machine learning for electronically excited states of molecules. Chemical Reviews. 2021;121:9873-9926. DOI: 10.1021/acs.chemrev.0c00749
  30. 30. Runge E, Gross EKU. Density-functional theory for time-dependent systems. Physical Review Letters. 1984;52:997-1000. DOI: 10.1103/PhysRevLett.52.997
  31. 31. Helbig N, Fuks JI, Casula M, Verstraete MJ, Marques MAL, Tokatly IV, et al. Density functional theory beyond the linear regime: Validating an adiabatic local density approximation. Physical Review A. 2011;83:032503-032505. DOI: 10.1103/PhysRevA.83.032503
  32. 32. Hedin L. New method for calculating the one-particle Green’s function with application to the electron-gas problem. Physics Review. 1965;139:A796-A823. DOI: 10.1103/PhysRev.139.A796
  33. 33. Leng X, Jin F, Wei M, Ma Y. GW method and Bethe–Salpeter equation for calculating electronic excitations. WIREs Computational Molecular Science. 2016;6:532-550. DOI: 10.1002/wcms.1265
  34. 34. Reining L. The GW approximation: content, successes and limitations. WIREs Computational Molecular Science. 2018;8:e1344-e1326. DOI: 10.1002/wcms.1344
  35. 35. Christensen AS, Kubař T, Cui Q, Elstner M. Semiempirical quantum mechanical methods for noncovalent interactions for chemical and biochemical applications. Chemical Reviews. 2016;116:5301-5337. DOI: 10.1021/acs.chemrev.5b00584
  36. 36. Dral PO, Wu X, Spörkel L, Koslowski A, Thiel W. Semiempirical Quantum-chemical orthogonalization-corrected methods: Benchmarks for ground-state properties. Journal of Chemical Theory and Computation. 2016;12:1097-1120. DOI: 10.1021/acs.jctc.5b01047
  37. 37. Tuna D, Lu Y, Koslowski A, Thiel W. Semiempirical Quantum-chemical orthogonalization-corrected methods: Benchmarks of electronically excited states. Journal of Chemical Theory and Computation. 2016;12:4400-4422. DOI: 10.1021/acs.jctc.6b00403
  38. 38. Sami S, Menger MFSJ, Faraji S, Broer R, Havenith RWA. Q-Force: Quantum Mechanically Augmented Molecular Force Fields. Journal of Chemical Theory and Computation. 2021;17:4946-4960. DOI: 10.1021/acs.jctc.1c00195
  39. 39. Dove MT. Introduction to Lattice Dynamics. Cambridge: Cambridge University Press; 1993. pp. 179-194. ISBN: 9780521392938
  40. 40. Wentzcovitch RM, Yu YG, Wu Z. Thermodynamic properties and phase relations in mantle minerals investigated by first principles quasiharmonic theory. In: Theoretical and Computational Methods in Mineral Physics: Geophysical Applications, Mineral. Chantilly, VA: Mineralogical Society of America; 2010. pp. 59-98. DOI: 10.1515/9781501508448 ISBN: 9780939950850
  41. 41. Errandonea D, Ferrer-Roca C, Martinez-García S, Segura A, Gomis A, Muňoz A, et al. High-pressure x-ray diffraction and ab initio study of Ni2Mo3N, Pd2Mo3N, Pt2Mo3N, Co3Mo3N, and Fe3Mo3N: Two families of ultra-incompressible bimetallic interstitial nitrides. Physical Review B. 2010;82:174105. DOI: 10.1103/PhysRevB.82.174105
  42. 42. Porezag D, Pederson MR, Liu AY. Importance of nonlinear core corrections for density-functional based pseudopotential calculations. Physical Review B. 1999;60:14132-14139. DOI: 10.1103/PhysRevB.60.14132
  43. 43. Giannozzi P et al. QUANTUM ESPRESSO: A modular and open-source software project for quantum simulations of materials. Journal of Physics: Condensed Matter. 2009;21:395502. DOI: 10.1088/0953-8984/21/39/395502. [Accessed on 21 November, 2021]
  44. 44. Available from: [Accessed on 21 November, 2021]
  45. 45. Perdew JP, Burke K, Ernzerhof M. Generalized gradient approximation made simple. Physical Review Letters. 1996;77:3865-3868. DOI: 10.1103/PhysRevLett.77.3865 Erratum: ibidem 1997;78:1396. DOI: 10.1103/PhysRevLett.78.1396
  46. 46. Murnaghan FD. The compressibility of media under extreme pressures. Proceedings of the National Academy of Sciences of the United States of America. 1944;30:244-247. DOI: 10.1073/pnas.30.9.244
  47. 47. Cohen RE, Gülseren O, Hemley RJ. Accuracy of equation-of-state formulations. American Mineralogist. 2000;85:338-344. DOI: 10.2138/am-2000-2-312
  48. 48. Otero-de-la-Roza A, Luaña V. Gibbs2: A new version of the quasi-harmonic model code. I. Robust treatment of the static data. Computer Physics Communications. 2011;182:1708-1720. DOI: 10.1016/j.cpc.2011.04.016
  49. 49. Togo A, Tanaka I. First principles phonon calculations in materials science. Scripta Materialia. 2015;108:1-5. DOI: 10.1016/j.scriptamat.2015.07.021
  50. 50. Kresse G, Furthmuller J, Hafner J. Ab initio force constant approach to phonon dispersion relations of diamond and graphite. Europhysics Letters. 1995;32:729-734. DOI: 10.1209/0295-5075/32/9/005
  51. 51. Baroni S, de Gironcoli S, Corso AD. Phonons and related crystal properties from density-functional perturbation theory. Reviews of Modern Physics. 2001;73:515-562. DOI: 10.1103/RevModPhys.73.515
  52. 52. Kolesnikov AI, Antonov VE, Efimchenko VS, Granrotha G, Klyamkin SN, Levchenko AV, et al. Neutron spectroscopy of magnesium dihydride. The Journal of Alloys and Compounds. 2011;509:S599-S603. DOI: 10.1016/j.jallcom.2010.10.156
  53. 53. Lasave J, Dominguez F, Koval S, Stachiotti M, Migoni RL. Shell-model description of lattice dynamical properties of MgH2. Journal of Physics: Condensed Matter. 2005;17:7133-7141. DOI: 10.1088/0953-8984/17/44/006
  54. 54. Mermin ND. Thermal properties of the inhomogeneous electron gas. Physics Review. 1965;137:A1441-A1443. DOI: 10.1103/PhysRev.137.A1441
  55. 55. Hofmann OT, Zojer E, Hörmann L, Jeindl A, Maurer RJ. First-principles calculations of hybrid inorganic–organic interfaces: From state-of-the-art to best practice. Physical Chemistry Chemical Physics. 2021;23:8132-8180. DOI: 10.1039/d0cp06605b
  56. 56. Ashcroft NW, Mermin ND. Solid State Physics. Philadelphia: Saunders College Publishing; 1976. ISBN: 0-03-083993-9, 9780030839931
  57. 57. Marx D, Hutter J. Ab Initio Molecular Dynamics: Basic Theory and Advanced Methods. Cambridge: Cambridge University Press; 2009. ISBN: 978-1-107-66353-4
  58. 58. Shchur LN. Journal of Physics: Conference Series. Vol. 1252. Bristol, United Kingdom: IOP Publishing Ltd; 2019. pp. 012010-012017. DOI: 10.1088/1742-6596/1252/1/012010
  59. 59. Born M, Huang K. Dynamical Theory of Crystal Lattices. New York: Oxford University Press; 1988. ISBN: 978-0198503699
  60. 60. Chen XJ, Zhang C, Meng Y, Zhang RQ, Lin HQ, Struzhkin VV. Mao HK, β-tin→Imma →sh phase transitions of germanium. Physical Review Letters. 2011;106:135502-135504. DOI: 10.1103/PhysRevLett.106.135502
  61. 61. Zhang H, Shang SL, Wang Y, Saengdeejing A, Chen LQ, Liu ZK. First-principles calculations of the elastic, phonon and thermodynamic properties of Al12Mg17. Acta Materialia. 2010;58:4012-4018. DOI: 10.1016/j.actamat.2010.03.020
  62. 62. Bian Q, Bose SK, Shukla RC. Vibrational and thermodynamic properties of metals from a model embedded-atom potential. Journal of Physics and Chemistry of Solids. 2008;69:168-181. DOI: 10.1016/j.jpcs.2007.08.046
  63. 63. Moser D, Baldissin G, Bull DJ, Riley DJ, Morrison I, Ross DK, et al. The pressure–temperature phase diagram of MgH2 and isotopic substitution. Journal of Physics: Condensed Matter. 2011;23:305403-305800. DOI: 10.1088/0953-8984/23/30/305403
  64. 64. Turney JE, Landry ES, McGaughey AJH, Amon CH. Predicting phonon properties and thermal conductivity from anharmonic lattice dynamics calculations and molecular dynamics simulations. Physical Review B. 2009;79:064301. DOI: 10.1103/PhysRevB.79.064301
  65. 65. Kohanoff J. Phonon spectra from short non-thermally equilibrated molecular dynamics simulations. Computational Materials Science. 1994;2:221-232. DOI: 10.1016/0927-0256(94)90103-1
  66. 66. Cao L, Stoltz G, Lelièvre T, Marinica MC, Athènes M. Free energy calculations from adaptive molecular dynamics simulations with adiabatic reweighting. The Journal of Chemical Physics. 2014;140:104108. DOI: 10.1063/1.4866811
  67. 67. Wang CZ, Chan CT, Ho KM. Tight-binding molecular-dynamics study of phonon anharmonic effects in silicon and diamond. Physical Review B. 1990;42:11276-11283. DOI: 10.1103/PhysRevB.42.11276
  68. 68. Debye P. Zur Theorie der spezifischen Wärmen. Annals of Physics. 1912;39:789-839. DOI: 10.1002/andp.19123441404
  69. 69. Anderson OL. A simplified method for calculating the Debye temperature from elastic constants. Journal of Physics and Chemistry of Solids. 1963;24:909-917. DOI: 10.1016/0022-3697(63)90067-2
  70. 70. Deus P, Schneider HA, Voland U. Estimation of the Debye temperature of diamond-like semiconducting compounds by means of the Lindemann rule. The Journal Crystal Research and Technology. 1981;16:951-948. DOI: 10.1002/crat.19810160814
  71. 71. Wolf U, Bohmhammel K, Wolf G. Supports open access. Thermochim Acta. 1998;310:37-42. DOI: 10.1016/S0040-6031(97)00382-1
  72. 72. NIST-JANAF. Thermochemical Tables. Last Update to Data Content. College Park, Maryland, United States Available from: American Institute of Physics; 1998 [Accessed on 21 November, 2021]
  73. 73. Barril X, Orozco M, Luque FJ. Predicting relative binding free energies of tacrine-Huperzine a hybrids as inhibitors of acetylcholinesterase. Journal of Medicinal Chemistry. 1999;42:5110-5119. DOI: 10.1021/jm990371u
  74. 74. Frenkel D, Smit B. Understanding Molecular Simulation: From Algorithms to Applications. Cambridge, Massachusetts, United States: Academic Press; 2002. DOI: 10.1016/B978-0-12-267351-1.X5000-7 ISBN: 978-0-12-267351-1
  75. 75. Gowdini E, Ahmad AA, Mabudi A, Hadipour NL, Kharazian B. A molecular dynamics study on the thermal properties of carbon-based gold nanoparticles. The Journal of Molecular Modeling. 2020;26:307-309. DOI: 10.1007/s00894-020-04559-2

Written By

Viorel Chihaia, Valentin Alexiev and Hasan S. AlMatrouk

Submitted: 22 December 2021 Reviewed: 28 February 2022 Published: 18 April 2022