Applications in General

#### **Chapter 1**

## Isothermal Calorimetry: Molecular Interactions between Small Molecules in Organic Solvents

*Raquel Gutiérrez-Climente, Elise Prost, Aude Cordin, Carlos Chesta and Luminita Duma*

#### **Abstract**

Isothermal titration calorimetry (ITC) is widely used to study protein-ligand, DNA-drug and/or protein-protein interactions but its application for small molecule complexation remains limited namely when the titration is performed in organic solvents. Compared to other dedicated spectroscopic techniques like nuclear magnetic resonance, infrared spectrometry or fluorimetry, which require a series of experiments to extract site-specific stoichiometry and affinity information, ITC provides in a single experiment a complete thermodynamic picture of the overall interaction mechanism. This chapter presents examples that support the high potential of ITC to probe interactions between small molecules in methanol, acetonitrile and methanol/ water mixture on a Nano ITC Low Volume device (TA Instruments), with an emphasis on both simple (1:1) and more complex (1:1 and 1:2) interaction mechanisms.

**Keywords:** isothermal calorimetry, thermodynamics, association constant, stoichiometry, small molecules, non-aqueous solvents

#### **1. Introduction**

Molecular recognition processes, omnipresent in nature, are of crucial importance in all living species. The magnitude of any molecular interaction, which can be translated in terms of heat released or adsorbed, can vary depending on the chemical nature of the interacting partners, on their concentration and the solvation environment in which the process takes place. Developed originally in the middle of the 1960s [1] for studying chemical reactions [2], ITC can measure with high precision and accuracy the heat energy associated with intermolecular reactions but also solvation and dilution experiments. Over the years, ITC became gradually widespread and popular to characterize the thermodynamics signature of molecular interactions in drug design [3, 4], in which the knowledge of the thermodynamic parameters in combination with the structural and kinetic information is decisive during the hit-tolead optimization, an early step in the drug design process. At this optimization stage, hundreds of compounds with promising affinity against the protein target are screened to identify the best one or two candidate molecules, usually from different

chemical series [5]. Compared to other methods employed in the field of drug discovery, ITC offers two clear advantages that facilitates the selection of the candidate molecules: On one side, ITC is the only method that directly measures the reaction enthalpy change [6], which can be considered as an interaction descriptor [7] and therefore an extremely useful parameter for structure thermodynamics correlations. It offers a valuable information notably for the compounds having similar binding affinities but different thermodynamic parameters [7, 8]. Traditionally, drugs characterized by an enthalpy-driven binding were preferred [7, 9]. Nevertheless, the nature of the binding site deserves to be considered because, when an apolar part of the drug interacts with an apolar region of the protein target, the entropic contribution will be favored and the reaction becomes entropy-driven [10]. On the other side, thanks to the detection of the heat change in water molecules and the transfer of protons of the drug molecule during the solvation and dilution experiments, ITC illustrates the differences between polar and apolar interactions that are invisible using techniques such as X-ray crystallography and surface plasmon resonance (SPR) [10]. All these advantages are also highly valuable for other ITC applications such as enzymecatalyzed reactions [11] and host-guest supramolecular complexation [12]. In supramolecular chemistry for example, ITC combined with supramolecular structure information can provide deeper information on the energies associated with noncovalent interactions (hydrogen bonding, electrostatic, π-π stacking, cation-π and anion-π interactions) and the hydrophobic effect induced by the displacement of water molecules [13].

Recently, the technique appeared also particularly useful and versatile in kinetic assays [14] where the direct measurement of a catalytic reaction [15, 16] was possible. Other applications include the monitoring of microbial activity and dynamics [17, 18], the stability assessment of (bio)pharmaceuticals [19], etc. The development of Low Volume (LV) Nano ITC calorimeters should extend the application fields of this technique not only to biomolecules available in small amounts but also to the study of complexation reactions in organic solvents thanks for example to the availability, on the same calorimeter, of a standard buret handle for aqueous solutions and organic solvents compatible buret handle.

The main advantage of ITC, when compared to other interaction-study dedicated spectroscopies, like nuclear magnetic resonance (NMR), Fourier transform infrared (FTIR) spectrometry or fluorimetry, is the ability to provide in a single experiment the entire thermodynamic profile of the investigated interaction process. In practice, an ITC experiment measures accurately the heat released or adsorbed when a molecule solution is titrated into another in a given aqueous or non-aqueous solution. The large majority of ITC measurements is conducted in aqueous solutions to study protein-ligand, DNA-drug and/or protein-protein interactions [10] whereas most of the ITC investigations in organic solvents or non-aqueous/aqueous mixtures focus on the solvation or dissolution thermodynamic study of various small molecules [20–23], drugs [24], single amino acids [25–27], small peptides [28], metal ions [29], etc. It is worth noting also some complexation studies for copper ions with β-alanine in ethanol [30], 15-crown-5 ether with Na<sup>+</sup> in water-ethanol [31], 18-crown-6 with triglycine in water-acetone and water-dimethyl sulfoxide [32], β-cyclodextrin with benzoic acid in water-ethanol [32], fluorescein isothiocyanate with polymers in water-methanol [33], and the association analysis of urea-based supramolecular polymers in different solvents [34].

The scarce resources in past calorimetry or contemporary ITC literature about intermolecular interactions in organic solvents prompted us to write this chapter.

*Isothermal Calorimetry: Molecular Interactions between Small Molecules in Organic Solvents DOI: http://dx.doi.org/10.5772/intechopen.104756*

After a brief introduction of the main basic concepts, we will show the ITC characterization of several 1:1 and 1:2 complexes in organic solvents together with the analysis of the data and the results which provide the thermodynamic profile, stoichiometry and association constant of the reaction. The correlation of ITC results with association or structural information from other spectroscopies may help understanding the formation of even more complex binding events.

#### **2. Isothermal calorimetry: general principles and experimental details**

Calorimetry measures the changes in heat released or absorbed during a chemical reaction or a physical process. Heats are measured with instruments called calorimeters classified as adiabatic or isothermal [35]. Heat measured on adiabatic devices show a permanent increase or decrease in temperature. The temperature changes thus obtained are directly related with the heat capacity of the instrument which therefore requires the acquisition of distinct experiments for calibration. In isothermal calorimeters the heat is allowed to flow between the reaction cell and a heat sink surrounding the two cells (see **Figure 1a**) and is actively regulated to maintain a constant level by power compensation. The development of micro- and nanocalorimetry rendered possible the detection of very small heat changes in small volumes of samples. Nano calorimeters have a detection limit in the nanowatt range [35]. Historically, Nano calorimeters have been applied to the study of the reactants available in reduced amounts like biomolecules.

The design of the NanoITC from TA Instruments relies on the differential power compensation technology (see **Figure 1a**) which allows to optimize sensitivity and responsiveness. Nanowatt sensitivity is achieved thanks to an internal reference. Two cylinder-shaped identical chambers (also called sample and reference cells) of 170 μL are located in a compartment which works as a thermal barrier. As described in the user manual [36], semiconducting thermoelectric devices (or TED) control and detect temperature differences between the sample and the reference chambers. In a titration experiment, both cells are entirely filled: the reference cell with the pure solvent, the sample cell with one reactant (or titrand) and the syringe with the other reactant (or titrant) of the reaction under study, both in the same solvent. The titrant is usually prepared at 10-fold higher concentration for a 1:1 binding model [37]. ITC titration experiments are implemented by incremental injection of a precise volume of titrant into the solution of titrand at discrete time intervals (see **Figure 1c**). Typically, about 25 injections per experiment are performed using a motor-driven syringe capable to deliver defined volume within 1–10 μL per injection. The syringe is coaxially introduced in the sample cell through a long access tube. The stepping motor precisely controls not only the injection volume but also the stirring speed of the reactants in the sample cell. Therefore, for each injection, the interaction between the two reactants releases (or adsorbs) heat that increases (or decreases) the sample cell temperature. This temperature change will activate the feedback heat controller (power compensation) on the sample cell such as to maintain a zero temperature difference between the two cells. For each heat variation during a stepwise titration, the feedback regulator will compensate this difference by decreasing (or increasing) the heat of the sample cell by the amount of heat supplied by the reaction and, at each variation associated to each injection step, will lead to a peak in the thermogram. **Figure 1c** illustrates the construction of a representative thermogram as the stepwise titration proceeds. During the titration experiment, the reactant in the sample cell is gradually

**Figure 1.**

*a) Drawing of the Nano ITC measuring unit (TA instruments) and its basic elements. b) Typical stepwise ITC raw thermogram and the corresponding integrated data. The parameters obtained after adjusting the data with the one-site binding model are also highlighted in green. c) Chemical reaction and illustration of the evolution of the stepwise raw thermograph as the titrant is injected in the sample cell.*

transformed into the molecular complex. The injection of the titrant is conducted until the titrand in the sample cell is fully saturated and the heat signal becomes equal to the background heat generally equal to the dilution heat of the titrant [37].

The normalized integrated area of each peak is next approached with the appropriate model to estimate the affinity, enthalpy and stoichiometry of the interaction. The first recorded experiment can be further analyzed with the Experiment Design tool of NanoAnalyze to find out the optimal concentration conditions leading to a "S-shape"

#### *Isothermal Calorimetry: Molecular Interactions between Small Molecules in Organic Solvents DOI: http://dx.doi.org/10.5772/intechopen.104756*

thermogram (see **Figure 1b**). The titration experiment is then repeated with these optimized concentrations and the data analyzed. The integration of each peak area in the thermogram (**Figure 1b** top) gives the amount of heat exchanged while injecting a known amount of the reaction partner into the sample cell. When the concentrations in the syringe and the sample cell were appropriately chosen, a display of the integrated heat signal as a function of mole ratio (of the injected compound over the one contained in the sample cell) reveals a typical sigmoidal shape (**Figure 1b** bottom). An inspection of the integrated titration curve gives information on the molar enthalpy (ΔH) from the height of the curve, the number of binding sites (or stoichiometry) of the reaction from the position of the inflection point on the mole ratio axis, and the association constant (Ka) from the slope of the curve. Nowadays, the companies provide an analysis software that incorporates the specifics of the instrument. It is therefore relatively easy to make a first estimation of the ΔH, n and Ka quantities thanks to various models that approach the experimental data using a nonlinear routine to find the most probable thermodynamic parameters describing the interaction process.

In the absence of any solubility or aggregation issues, direct and reverse titration [38, 39] can be considered when studying the interaction of small molecules in a given solvent. For two A and B interacting molecules (**Figure 1c**), the injection of aliquots of A into B is called direct titration whereas the injection of aliquots of B into A is termed reverse titration. As previously mentioned, the titration experiment is well performed when the concentrations are chosen such as the molecule in the sample cell is fully saturated. Assuming that the molecular mass of reactants does not affect the equilibrium binding equations (microscopy reversibility), direct and reverse titrations should be equivalent. They can be modeled by the same set of equations and should lead to the same thermodynamic parameters for a reversible reaction.

The reaction in **Figure 1c** corresponds to a complex where a single molecule A interacts with a single molecule B (i.e., the one to one model which corresponds to a stoichiometry of 1). Since the reaction enthalpy change (ΔH) and association constant (Ka) are directly measured from an ITC titration experiment, under constant temperature and pressure, the Gibbs free energy (ΔG) and the entropy (ΔS) can obtained using the following relationships:

$$\begin{cases} \triangle \mathbf{G} = \triangle \mathbf{H} - \mathbf{T} \triangle \mathbf{S} \\ \triangle \mathbf{G} = -\mathbf{R} \mathbf{T} \ln \mathbf{K\_a} \end{cases} \tag{1}$$

with R the ideal gas constant (R = 8.314 J/mol�K) and T the temperature in K. Equilibrium association constants can be obtained with acceptable statistical precision if the following condition is fulfilled:

$$1 < c < 1000 \tag{2}$$

where *c* (also known as *c*-value) is given by the relationship *c* ¼ nKa½ � B , [B]: concentration of B molecule [40]. The condition in Eq. (2) is directly related to the shape of the binding curve (also termed thermogram). For protein-ligand complexes generally characterized by an affinity in the μM range (i.e. 106 M�<sup>1</sup> Ka), a *c*-value within 20 and 200 is recommended in order to minimize Ka errors and thermodynamic parameters as shown by comparing ITC reports from different laboratories on the same benchmark protein-ligand complex in buffer [41–43].

By performing the ITC titration experiment at different temperatures, the change in heat capacity, ΔCp, for the A�B complex formation can be estimated according to the following equation:

$$
\Delta \mathbf{C\_p} = \frac{\mathbf{d} \Delta \mathbf{H}}{\mathbf{d} \mathbf{T}} \tag{3}
$$

with ΔCp in J/mol�K.

The enthalpy represents the energy change of the system when the molecule A interacts with B in a given solvent. Different type of noncovalent interactions (hydrogen bonds, ions pairs, van der Waals forces, etc.) can take place at the binding interface and therefore affect the enthalpy change. For example, the formation of noncovalent interactions between atoms is an exothermic process characterized by a negative enthalpy change whereas their breaking is an endothermic one with a positive enthalpy change. The heat released during a binding process describes the entire system under study with individual contributions from the interacting partners but also the solvent. In reality, the measured enthalpy change upon binding is the sum of many positive and negative ΔH contributions resulting from the simultaneous formation and disruption of noncovalent interactions [44]. Like in protein-ligand interactions in an aqueous medium, the observed ΔH of binding in an organic solvent is a global property which reflects the partial loss of solvent contacts of the interacting partners, the formation of complex noncovalent interactions and the solvent rearrangement near the complex surface.

As main direct experimental observable in calorimetry, the measured heat is correlated with the reaction taking place at the molecular level and the aim of the calorimetry is to provide reliable heat data capable to characterize molecular interactions. In the literature, the enthalpies directly obtained by calorimetry have been correlated with the values of binding enthalpy derived from the van't Hoff relationships. Experimental and simulation studies have shown that statistically relevant discrepancies are notably found when the experimental setup or data analysis are not correctly performed [45, 46].

Incremental titration described previously is the most common titration method used. Continuous titrations, which consists of constantly injecting the titrant into the calorimeter vessel while monitoring the thermal power, are shorter than the incremental ones and therefore, they can be of great interest for unstable samples. The development of the continuous ITC (cITC) method for micro- and nano-calorimeters [47] rendered the technique even faster and more versatile for the study of thermodynamic processes in a complex interaction. Interestingly, it can also represent a quick alterative to find out the concentration conditions leading to an exploitable thermogram. The screening of the optimum conditions (Eq. 2) may be speeded up even more if the cITC thermograms can be exploited by the Experiment Design tool in NanoAnalyze software, as mentioned for the classical ITC data. Unfortunately, this is not yet the case with the currently available NanoAnalyze software (3.12.0).

Another advantage of the cITC, is the potential expansion of the equilibrium constants accessible by ITC. Indeed, Markova and Hallén [47] have shown by computer simulations that cITC expands by 3 orders of magnitude the range of Ka achievable by ITC. Therefore, for cITC the Eq. (2) becomes:

$$1 < \mathbf{c} < \mathbf{3x10}^{\delta} \tag{4}$$

*Isothermal Calorimetry: Molecular Interactions between Small Molecules in Organic Solvents DOI: http://dx.doi.org/10.5772/intechopen.104756*

with 10<sup>12</sup> M<sup>1</sup> being the highest equilibrium constant reachable by cITC. This is rendered possible by the increased data points density in a continuous titration experiment which better defines the slope region in a 1:1 binding curve.

#### **3. Experimental details**

The reagents used in this work are summarized in **Table 1**. All measurements were performed on a differential power compensation Nano ITC Low Volume (Nano ITC LV) calorimeter from TA Instruments (Waters, France), using a 50 μL injection syringe while stirring at 400 rpm (the maximum stirring value). The Nano ITC LV has two gold 170 μL reaction vessels and a buret assembly holding a stainless-needle syringe with a twisted paddle at the tip and titrant exit at the bottom. Therefore, the syringe serves not only for the delivery of the titrant but also for stirring. To avoid the formation of bubbles in the cells and syringe, the samples were degassed in a vacuum degassing station [48] for 15 minutes immediately before use. Injections were started after achievement of baseline stability (using the automatic equilibration mode for "small Heats" of ITCrun program controlling the calorimeter with the following criteria: absolute acceptable slope: ΔH 0.1 μW/h; acceptable absolute standard deviation: 0.01 μW). In addition, an equilibration time of 300 s has been considered before the first and after the last injection to assess the quality of the baseline. The experimental parameters were: 350 μL in the sample cell, 50 μL in the syringe, 350 μL solvent in the reference cell (changed weekly), 25 injections of 2.02 μL except for the first injection which was of 0.48 μL. The integrated heat effects of each injection were corrected by subtraction of the corresponding integrated heat effects associated with titrant dilution into the solvent. The experimental data obtained from the corrected calorimetric titration were analyzed on the basis of different interaction models with the NanoAnalyze software. The first injection was not taken into consideration for data analysis.

The electrical calibration of the calorimeter was performed according to the manufacturer's instructions. Water in water dilution experiments are regularly performed to check and validate the initial criteria of the manufacturer. Cells and syringe correct cleaning is essential to avoid artifacts and produce good quality data. Additionally, it is crucial to keep the needle of the syringe perfectly straight. Cleaning of the sample cell can be performed automatically using the vacuum of the degassing device and is generally done with 1 L of 2.5% DECON followed by 1 L of MilliQ water. An ITC


N*-[3-(dimethylamino)propyl]methacrylamide.*

#### **Table 1.** *Reagents used in this work.*

experiment takes between 1 and 2 hours depending on the injection delay and the equilibration duration necessary to fulfill the heat stability statistical criteria. ITCRun, the software which controls the Nano ITC calorimeter, does not record the evolution of the heat values during the equilibration delay but only its duration.

Despite the relative simplicity of ITC experiments, the selection of the right binding model for the fitting of the experimental data in order to estimate the thermodynamic parameters can be challenging, especially for complexes where one of the reactants present multiple sites and there is not previous information about the stoichiometry or binding mechanism [49]. Among the various softwares currently available for the titration data analysis, the softwares provided by the ITC manufacturers, i.e. Origin from MicroCal/Malvern and NanoAnalyze in the case of TA instruments, the interpretation of the binding isotherms can be done for each data set individually using either classical models such as one-site independent or two-sites sequential models, multiple sites, dimer dissociation, cooperative and competitive replacement models. The possibility of dilution subtraction or the selection of a control model together with the binding model remains at the user's choice.

Over the last decade, alternative softwares have been developed for more complex processes, e.g. AFFINImeter [50], pytc [51], HypΔH [52], CHASM [53] or SEDPHAT [54]. The free platform SEDPHAT gives the possibility to combine several experimental data (calorimetry, spectrophotometry, sedimentation and surface binding assays) in a single global analysis. The main purpose of the platform development was to reduce the discrepancies between the thermodynamic parameters obtained using different devices and setups (e.g. the use of different sample volume, concentrations, immobilization of the template, etc...). The SEDPHAT results presented in the present chapter concern only calorimetry data in order to perform a global analysis of several measurements, including repetitions of the performed direct and reverse titrations, increasing therefore the confidence in the binding parameters obtained individually with NanoAnalyze software.

#### **4. ITC data and their evaluation**

The following subsections illustrate the extraction of thermodynamic parameters and association constant from the experimental raw data for one-site and two-sites binding complexes in different non-aqueous solvents.

#### **4.1 1:1 Complexes in different solvents and calorimetric heat capacity**

Molecular recognition processes are the archetypal reactions in various domains going from life science to technology and the design of a suitable target which binds the other partner with specificity remains challenging nowadays. The study of the binding complexes presented herein was originally motivated by the need to optimize the design of molecularly imprinted polymers [55] by characterizing the affinity and the complete thermodynamic profile of the monomer-target interaction in the solvent used for the synthesis of the final polymer. A better comprehension of the interactions and energies involved in the preorganization of the monomers around the target could improve their synthesis protocols by selecting monomers with high affinity for the target.

The feasibility of elucidating the thermodynamic parameters for a 1:1 binding model using stepwise titration experiments is first demonstrated with the interaction *Isothermal Calorimetry: Molecular Interactions between Small Molecules in Organic Solvents DOI: http://dx.doi.org/10.5772/intechopen.104756*

#### **Figure 2.**

*ITC titration of 2 into 1 in MeOH at 15 (a) and 35°C (b) as filled black circles. The dilution of 2 into MeOH is also shown as gray filled circles. The middle thermogram shows the integrated heat data as filled black circles and their nonlinear fitting using the 1:1 binding model as red continuous line. c) Integrated corrected data at 10 and 35°C. d) Linear fit of the calorimetric ΔH values versus temperature which gives access to the change in heat capacity, ΔCp 0.5 kJ/molK.*

between benzoic acid **(1)**, a food antimicrobial, and *N*-[3-(dimethylamino)propyl] methacrylamide (**2**), a monomer, in methanol (**Figure 2**). In water at 25°C, **1** shows a pKa of 4.2 [49].

whereas **2** is a typical amine base characterized by a pKa of 9.2. In methanol and other organic solvents, the pKa cannot be defined but it is expected that **1** and **2** behave as a strong Lewis acid and base, respectively. Qualitative evidence of the interaction between benzoic acid and 2-dimethylaminoethyl methacrylate-based copolymers has been provided by infrared (IR) and <sup>1</sup> H NMR spectroscopies [50].

1 H NMR spectra recorded for different molar compositions of benzoic acid:copolymer containing amine groups in acetone showed that carboxyl-amine interaction causes a shift of the aromatic proton resonances to high field (small ppm values) whereas protons near the amine group shift to low field (big ppm values). These chemical shifts suggest that the carboxyl-amine complexation is associated with an increase in electron density for the carboxyl and a decrease of the electron density for the amine group, results consistent with the formation of a contact ion pair stabilized by attractive electrostatic forces and H-bonding. The same conclusions are reached by analyzing the IR spectra of the mixtures in acetone. For these complexes in chloroform and acetone solvents, it was stated that the excess of benzoic acid tends to self-association producing therefore benzoic acid dimers in addition to the benzoic acid-amine complex.

**Figures 2** and **3** display the calorimetry data for the interaction between benzoic acid (**1**) and *N*-[3-(dimethylamino)propyl]methacrylamide (**2**) in methanol and acetonitrile, respectively. In order to estimate the heat capacity changes, the heat reaction response has been measured in methanol over a temperature range between 10 and 35°C. The top panels in **Figure 2a** and **b** indicate that the complexation is exothermic over the full temperature range spanned. The middle and the bottom panels show the integrated experimental heat data fitted with the 1:1 independent binding model and the residuals which describe the differences between the interaction model and the measured data. In all experiments, 60 mM **2** were titrated into 8 mM **1** except

#### **Figure 3.**

*ITC titration of 2 into 1 in acetonitrile at 25°C as filled black circles. The dilution of 2 into acetonitrile is overlaid as gray filled circles. The middle thermogram shows the integrated heat data as filled black circles and the nonlinear fitting of the data using the 1:1 binding model as red continuous line.*

*Isothermal Calorimetry: Molecular Interactions between Small Molecules in Organic Solvents DOI: http://dx.doi.org/10.5772/intechopen.104756*


*a kJ/mol, <sup>Δ</sup>S is given in J/mol K. <sup>b</sup>*

*calculated* <sup>c</sup>*-value. <sup>c</sup>*

*100 mM <sup>2</sup> into 10 mM <sup>1</sup>. <sup>d</sup>*

*60 mM <sup>1</sup> into 8 mM <sup>2</sup> (reverse titration). <sup>e</sup>*

*global analysis done on direct and reverse datasets and performed with SEDPHAT platform [54], n value is calculated based on the values obtained in the individual analysis with NITPIC [56] as SEDPHAT employs a different parameter in its analysis. f*

*in acetonitrile.*

#### **Table 2.**

*Best-fit thermodynamic parameters from ITC measurements for the 1–2 complex in methanol.*

otherwise indicated. **Table 2** summarizes the ensemble of thermodynamic parameters together with the stoichiometry and the association constant derived using a 1:1 independent binding model at each temperature. The calculated *c*-value is also shown and falls in the range 2–20. The concentration conditions have been optimized using the tool "ITC Experiment Design" from NanoAnalyze at 25°C for a *c*-value 10. A greater *c*-value would have required higher concentrations for benzoic acid and therefore would have prevented us from carrying out the reverse titration owing to benzoic acid reduced solubility in MeOH above 60 mM. In parallel, to check the accuracy of the **1**–**2** affinity and thermodynamics parameters obtained by ITC, a global fit analysis for both direct (60 mM **2** into 8 mM **1**) and reverse (60 mM **1** into 8 mM **2**) titrations in MeOH at 25°C has also been performed with the open source SEDPHAT platform. First, the ITC data have been integrated using the NITPIC software which also allows to subtract the corresponding dilution for each direct and reverse titration dataset and gives access to the stoichiometry parameter. Second, the direct and indirect NITPIC datasets were integrated and saved in a SEDPHAT configuration file. The data were then fitted with the one-site binding model and the statistics of the thermodynamic parameters calculated for a confidence level at 95%. The results obtained with SEDPHAT (**Table 2**) are almost identical to the ones obtained for the direct titration with NanoAnalyze. SEDPHAT contains explicit factors which can account for errors in the active concentrations. Their inspection for the individual analysis of both and direct titration data suggests an error in concentrations notably for the indirect titration. This might explain the bigger estimations of the stoichiometry and enthalpy change in the case of reverse titration.

The results summarized in **Table 2** show that the **1**–**2** complexation is an exothermic (ΔH < 0) and enthalpy-driven (|ΔH| > |TΔS|) process at all temperatures and experimental conditions studied. Ka decreases with increasing temperature while ΔH becomes more exothermic as the temperature rises. The |ΔS| show also a tendency to increase with temperature. Within the experimental uncertainties, the calculated |ΔS| values, which range within 20 and 50 kJ/molK, are in good agreement with the expected entropy change characterizing the ion pair formation from neutral reactants in polar solvents [57, 58]. The decrease of the association constant by about one order of magnitude between 10 and 35°C is particularly relevant for the design of molecularly imprinted polymers because their synthesis is usually performed at temperatures higher than the ones used here (50 to 70°C). It therefore suggests the importance of probing the interaction between the target and the functional monomer at the synthesis temperature.

In order to estimate the heat capacity changes, the heat reaction response has been measured over a temperature range between 10 and 35°C. The analysis of the variation of ΔH as function of T (**Figure 2d**) allows to estimate a value for ΔCp 0.5 kJ/ molK. Formally, for a simple equilibrium studied over a relatively narrow temperature range, ΔCp should be zero. Thus, these results suggest that the complexation process may be more complicated than initially assumed. Taken together, these results suggest that the **1**–**2** complex formation mechanism probably involves two or more coupled equilibria whose relative importance depends on temperature and reactant concentrations. For example, if at the concentrations used, a small fraction of the benzoic acid is in the dimeric form, the **1**–**2** acid/base complex formation requires the dissociation of the dimer and therefore a coupling of the two equilibria (complexation and dimerization). This can explain the discrepancies observed in n (which is clearly different of 1), ΔH and ΔS calculated at different temperatures (**Figure 2a** and **b**). Coupled equilibria may also explain the differences between "direct" and "reverse" titrations.

**Figure 3** shows ITC titration of 60 mM of **2** in acetonitrile in 8 mM of **1** in the same solvent at 25°C. The thermodynamic parameters obtained from fitting the experimental data are listed in **Table 2**. This study was performed to investigate the role of H-bond interactions in the thermodynamics of the reaction. Acetonitrile and methanol show similar dielectric constants (36.6 and 33.0, respectively) and thus should exhibit similar abilities to stabilize the ionic product (relative to the reactants). However, both solvents differ markedly in their ability to form H-bonds. Methanol, a protic solvent, should further stabilize reactants and reaction product through the formation of such bonds. However, a comparative analysis of the ΔG values (ΔH and ΔS) obtained in acetonitrile and methanol at 25°C (see **Table 2**) shows that, within experimental uncertainties, they are practically identical. These results suggest that H-bonds do not preferentially stabilize neither the reactants nor the complex, and therefore, they are not determinant in the thermodynamics of the complexation process.

The second complexation process explored by stepwise ITC titration concerns the interaction between the glucuronic acid (**4**), a sugar acid derived from glucose, and (4-acrylamidophenyl)amino methaniminium acetate (**5**), a polymerizable benzamidine salt, which can form with carboxylates stoichiometric non-covalent complexes characterized by affinities higher than 103 M<sup>1</sup> [59]. We previously studied the **4**–**5** complexation by <sup>1</sup> H NMR spectroscopy in DMSO-*d6* [60]. Job's plot [61, 62] and titration experiments demonstrated the formation of a 1:1 complex with an affinity of 7.1 10<sup>3</sup> M<sup>1</sup> . The complex formation in MeOD/D2O (4/1 v/v) gives slightly smaller chemical shift differences and an association constant of 4.4 10<sup>3</sup> M<sup>1</sup> .

**Figure 4** displays the calorimetry data obtained at 25°C in MeOH/H2O (4/1 v/v) for direct (i.e., monomer into the glucuronic acid) and reverse (i.e., glucuronic acid *Isothermal Calorimetry: Molecular Interactions between Small Molecules in Organic Solvents DOI: http://dx.doi.org/10.5772/intechopen.104756*

#### **Figure 4.**

*Direct (a) and reverse (b) titrations of AB<sup>+</sup> Ac (5) into GlucH (4) in MeOH/H2O (4/1 v/v) at 25°C as filled black circles. The dilutions of 4 and 5, respectively, into MeOH/H2O (4/1 v/v) are also shown as gray filled circles. The middle thermogram shows the integrated heat data as filled black circles and the nonlinear fitting of the data using the 1:1 independent binding model as red continuous line.*

into the monomer) titrations. An inspection of the reaction scheme (top of **Figure 4**) suggests that the process can be considered as a displacement reaction (i.e., a reaction where glucuronic acid replaces acetate as a counter ion of the salt). The titration thermograms in top panels of **Figure 4a** and **b** show that the macroscopic experimental heat either changes from endothermic to exothermic (direct titration) or remains


*a kJ/mol, <sup>Δ</sup>S is given in J/mol K. <sup>b</sup>*

*calculated* <sup>c</sup>*-value. <sup>c</sup>*

*21 mM <sup>5</sup> into 3 mM <sup>4</sup> (direct titration). <sup>d</sup>*

*25 mM 4 into 3 mM 5 (reverse titration).*

#### **Table 3.**

*Thermodynamic parameters of the 4–5 complex in MeOH/H2O (4/1 v/v) at 25°C.*

**Figure 5.** *Global displacement reaction described as a 3 steps chemical process.*

endothermic (reverse titration). The dilution of **4** and **5** is exothermic. The nonlinear fitting of the macroscopic heat data after dilution subtraction gives the thermodynamic parameters, the stoichiometry and the association constant reported in **Table 3**. The relatively small errors can be correlated with the calculated *c*-value, which is an order of magnitude higher than for the **1**–**2** complex described previously.

For both (direct and reverse) titrations, the enthalpy change (ΔH) is close to zero whereas the Gibbs free energy change (ΔG) is negative, which indicates a spontaneous process. Since T|ΔS| > > |ΔH|, the reaction is entropy driven. As we show below, the negative driving force of the displacement reaction is mainly due to the observed differences in the association constant (K<sup>0</sup> ) of the two acids. In water, the reported pKa are 4.75 [63] and 2.88 [64] for acetic and glucuronic acids, respectively. This difference in acidity is explained by the presence of one oxygen atom in the carbon α to the glucuronic -COOH group, which significantly weakens the O-H bond. The thermodynamics of the displacement reaction can be analyzed as a reaction occurring in consecutive stages. This is possible because Gibbs free energy is a state function of the system (i.e., it is independent pathway taken but only on the initial and final states) and therefore, the global displacement reaction (top of **Figure 5**) can be written as the sum of 3 hypothetical steps reaction characterized by their Gibbs free energies: ΔG1, ΔG2 and ΔG3, so that ΔG = ΔG1 + ΔG2 + ΔG3.

The first step describes to the dissociation of the ABþAc� into the corresponding solvated free ions characterized by ΔG1 > 0. The second step is an acid-base reaction in which the acetate and the glucuronic acid exchange a proton for which the ΔG can be calculated as explained latter. Finally, the third step corresponds to the association of the glucoronate ion with the monomer counter ion to complete the displacement process and has a ΔG3 < 0. Although ΔG1 and ΔG3 are not identical, as they characterize the dissociation/association of different ion pairs, they will have similar orders of magnitude and we can assume that: ΔG1 + ΔG3 � 0. Based on this hypothesis, the main contribution to the global reaction should come from ΔG2. The equilibrium constant for this process (K2) can been determined from the pKa values to be �80. Using the relationships of Eq. (1), we obtain ΔG2 � � 10 kJ/mol. Although this value is half of that obtained experimentally, it should be remembered that the K2 was roughly calculated from the pKa obtained in water. The decreased dielectric constant in the MeOH/H2O (4/1 v/v) mixture used in this study should enhance the differences in acidity between the two acids, making the process even more spontaneous.

*Isothermal Calorimetry: Molecular Interactions between Small Molecules in Organic Solvents DOI: http://dx.doi.org/10.5772/intechopen.104756*

**Figure 6.**

*ITC titration for 1:1 and 1:2 complexes in methanol at 25°C. a) 2 into 1 and b) 2 into (3) titrations as black filled circles. The dilution of 2 into methanol is shown as grey filled circles on the top graphs as well. The middle thermogram shows the integrated heat data as filled black circles and the nonlinear fitting of the data using the onesite binding model as red continuous line or the sequential two-sites binding model for the interaction with isophthalic acid (right side) as green continuous line.*

#### **4.2 Two-sites complexes**

The involvement of more than one molecular process into a reaction is not always directly visible on the measured heat profile. The casual case where the two interaction mechanisms have an enthalpy that nicely emphasizes each inflection of 1:2 (or 1:3) complexation is relatively rare. **Figure 6** shows the example of a two-sites assembly with the interaction between the dimethylamino-based monomer and isophthalic acid, a target with two carboxylic groups, by putting it in mirror of the complex formed between the same monomer with benzoic acid. Both titrations have been carried out at 25°C in methanol.

The experimental heat response for a 2:1 complexation could be in theory approached with three different binding models (one-site independent, two-sites independent and two-sites sequential). Nevertheless, different arguments may help to restrict the plausible models: close inspection of the data to see if the two individual binding processes are clearly visible on the heat response, choice of the model with the smallest number of adjustable parameters, analysis of the nature of the functional groups of the interacting partners, interaction knowledge from other techniques, etc.

For all the thermograms in **Figure 6**, the binding isotherms were obtained as previously described by integrating the individual heat variation which results from the injection of the titrant. In both cases, initially a nonlinear least-squares regression analysis considering a single binding site was employed for both titrations. As represented in **Figure 6a**, the model fits well the data for the complex with benzoic acid, allowing to estimate the dissociation constant and the binding enthalpy while the stoichiometry value close to 1 confirmed the one-site interaction. In the case of the complex with the isophthalic acid, the one-site binding model is the simplest one to test with only 3 variable parameters and could represent a first acceptable choice if we assume the two carboxyl groups as chemically equivalent. However, this assumption is not strictly true, since it has been reported that the pKa values of the two carboxyl groups of **3** at 25°C in water are: pKa1 = 3.46, pKa2 = 4.46 (ref). Therefore, this model is expected to estimate a stoichiometry of 2 together with a binding enthalpy (ΔH) and an association constant (Ka) as an average of the two complexation events. **Table 4** collects the thermodynamic parameters


*a kJ/mol except for <sup>Δ</sup>S which is given in J/mol K. <sup>b</sup>*

*calculated* <sup>c</sup>*-value. <sup>c</sup>*

*one-site independent binding model.*

*d sequential two-sites binding model.*

#### **Table 4.**

*Thermodynamic parameters obtained by nonlinear fitting of heat data corresponding to the titration of 60 mM 2 into 8 mM 1 (first row) and 150 mM 2 into 10 mM 3 (last two rows) in MeOH at 25°C.*

obtained by nonlinear fitting the data with the one-site binding model. It is important to note that the thermodynamic parameters shown in **Table 4** are per mole of binding sites (n = 2), so the global association constant (which considers full site occupancy) is (Ka) <sup>2</sup> 105 <sup>M</sup><sup>2</sup> and the global values of <sup>Δ</sup>G, <sup>Δ</sup>H and <sup>Δ</sup>S are twice those reported, i.e.: (28 4) and -(34 2) kJ/mol and (20 20) kJ/mol K, respectively.

The sequential two-sites model was also employed to approach the experimental ITC data. The evolution of the residuals and the decrease of the standard deviation suggest that the sequential two-sites binding model better approximates the experimental titration data. This model assumes the existence of two distinct binding sites in **3** with 1:1 stoichiometry each (i.e., the population of type 1 sites and type 2 sites in the sample is identical). The fitting of the experimental data to the model therefore provides 4 thermodynamic parameters, 2 for the initial stage (occupancy of site 1) and 3 for the second stage (occupancy of site 2). Results are summarized in **Table 4**. Again, the global equilibrium constant is the product of the Ka obtained at each step (Ka1xKa2) 105 <sup>M</sup><sup>2</sup> , which is in full agreement with the result obtained using de 1:1 model. Similarly, the overall ΔG, ΔH and ΔS are the sum of the individual contribution to give: (38 6) and -(42 4) kJ/mol and (50 20) kJ/mol K, respectively. Within the experimental uncertainties, these values are similar to the ones obtained using the simplest one-site independent model. When the results determined with the sequential model for the **3**–**2** reaction are compared with those obtained for the **1**–**2** reaction, modeled with the onesite independent model, some interesting aspects emerge. For example, the first reaction step between **2** and **3** is spontaneous and its ΔG 17 kJ/mol is very similar to that obtained for the reaction of **2** with **1**, ΔG 16 kJ/mol. However, such coincidence occurs due to the compensation between ΔH and ΔS values (see **Table 4**), the entropic change for the first step of the reaction being for the reaction **3**–**2** positive. The results also suggest that the second reaction step is slightly less spontaneous, but is accompanied by a very large negative enthalpy and entropy change. Additional experimental data are required to better understand this system.

We should also mention that the two-sites independent model has also been tested but it has been discarded in the end as it contains the highest number of adjustable parameters (6 in total) and was giving unrealistic stoichiometries for the two sites.

#### **5. Conclusion and perspectives**

Isothermal calorimetry is a label free non-destructive technique that became over the years the method of choice for most binding studies in solution. Compared to other

#### *Isothermal Calorimetry: Molecular Interactions between Small Molecules in Organic Solvents DOI: http://dx.doi.org/10.5772/intechopen.104756*

dedicated spectroscopic techniques like nuclear magnetic resonance, infrared spectrometry or fluorimetry, which require a series of experiments to extract site-specific stoichiometry and affinity information, ITC gives access - in a single experiment – to a complete thermodynamic picture of the overall interaction mechanism. The continuous development of the technique led to highly sensitive instruments capable to detect small amounts of samples. For example, isothermal micro- and nano-calorimetry can measure tiny heat changes and became an important method for the thermodynamics profile analysis of chemical and biochemical reactions, solvation and dissolution processes. The automation of the entire experimental process with the new generation of calorimeters optimizes not only the use of the calorimeter but also increases the number of samples that can be studied. In parallel, considerable efforts have been made to propose various binding models and softwares with powerful routines for the optimization of experimental conditions (notably in terms of concentrations).

Despite the widespread application of ITC to probe interactions between biomolecules, the technique is hardly ever used for the characterization of the complexation of small molecules, especially when the titration is performed in organic solvents. Herein, we presented few examples that support the high potential of ITC for the study of interactions between small molecules in methanol, acetonitrile and methanol/water mixture on a Nano ITC Low Volume device (TA Instruments), with an emphasis on both simple (1:1) and more complex (1:1 and 1:2) interaction mechanisms.

In addition to the binding studies at equilibrium, ITC provided promising results for the investigation of reaction kinetics, irreversible reactions, reactions under pressure, etc. The field of application of isothermal calorimetry is continuously expanding from pharmacology to life science, clinical medicine, environmental science, biotechnology, ecology, etc.

#### **Acknowledgements**

The authors thank Jacques Loubens at TA Instruments, Waters, France for his kind assistance with the ITCRun and data analysis NanoAnalyze softwares. We are grateful to Bernadette Tse Sum Bui for kindly providing the monomer (4-acrylamidophenyl) aminomethaminium acetate. We thank the European Regional Development Fund and the Region of Picardy (CPER 2007-2020). L.D. acknowledges financial support from the European Commission (project NOSY for New Operational Sensing sYstems, Grant agreement ID 653839, H2020-EU.3.7), the Hauts-de-France Region, the European Regional Development Fund (ERDF) 2014/2020 and Idex Sorbonne Université Investissements d'Avenir (2019/2020 EMERGENCE program).

#### **Conflict of interest**

The authors declare no conflict of interest.

*Applications of Calorimetry*

### **Author details**

Raquel Gutiérrez-Climente1†, Elise Prost2†, Aude Cordin<sup>2</sup> , Carlos Chesta<sup>3</sup> and Luminita Duma<sup>4</sup> \*

1 University Montpellier, IBMM, CNRS, ENSCM, Montpellier, France

2 University of Technology of Compiègne, UPJV, UMR CNRS 7025, Enzyme and Cell Engineering, Research Centre Royallieu, Compiègne, France

3 National Universiy of Río Cuarto, IITEMA-CONICET, Río Cuarto, Argentina

4 Champagne-Ardenne University, CNRS, ICMR UMR 7312, Reims, France

\*Address all correspondence to: luminita.duma@univ-reims.fr

† These authors contributed equally.

© 2022 The Author(s). Licensee IntechOpen. This chapter is distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/3.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

*Isothermal Calorimetry: Molecular Interactions between Small Molecules in Organic Solvents DOI: http://dx.doi.org/10.5772/intechopen.104756*

#### **References**

[1] Christensen JJ, Izatt RM, Hansen LD. New precision thermometric titration calorimeter. The Review of Scientific Instruments. 1965;**36**(6):779-783. DOI: 10.1063/1.1719702

[2] Hansen LD, Christensen JJ, Izatt RM. Entropy titration. A calorimetric method for the determination of ΔG, ΔH and ΔS. Chemical Communications. 1965;**70**(3):36-38. DOI: 10.1039/C1965 0000036

[3] Chaires JB. Calorimetry and thermodynamics in drug design. Annual Review of Biophysics. 2008;**37**:135-151. DOI: 10.1146/annurev. biophys.36.040306.132812

[4] Rizzuti B, Lan W, Santofimia-Castaño P, Zhou Z, Velázquez-Campoy A, Abián O, et al. Design of inhibitors of the intrinsically disordered protein nupr1: Balance between drug affinity and target function. Biomolecules. 2021;**11**(10):1-19. DOI: 10.3390/ biom11101453

[5] Hughes JP, Rees S, Kalindjian SB, Philpott KL. Principles of early drug discovery. British Journal of Pharmacology. 2011;**162**(6):1239-1249. DOI: 10.1111/j.1476-5381.2010.01127.x

[6] Linkuvienė V, Krainer G, Chen W-Y, Matulis D. Isothermal titration calorimetry for drug design: Precision of the enthalpy and binding constant measurements and comparison of the instruments. Analytical Biochemistry. 2016;**515**:61-64. DOI: 10.1016/j.ab. 2016.10.005

[7] Ladbury JE, Klebe G, Freire E. Adding calorimetric data to decision making in lead discovery: A hot tip. Nature Reviews. Drug Discovery. 2010;**9**(1): 23-27. DOI: 10.1038/nrd3054

[8] Velazquez-Campoy A, Luque I, Freire E. The use of isothermal titration calorimetry in drug design: Applications to high affinity binding and protonation/ deprotonation coupling. Netsu Sokutei. 2001;**28**(2):68-73. DOI: 10.11311/ jscta1974.28.68

[9] Freire E. Do enthalpy and entropy distinguish first in class from best in class? Drug Discovery Today. 2008; **13**(19):869-874. DOI: 10.1016/ j.drudis.2008.07.005

[10] Falconer RJ, Schuur B, Mittermaier AK. Applications of isothermal titration calorimetry in pure and applied research from 2016 to 2020. Journal of Molecular Recognition. 2021;**34**(10):e2901. DOI: 10.1002/jmr.2901

[11] Vander Meulen KA, Horowitz S, Trievel RC, Butcher SE. Measuring the kinetics of molecular association by isothermal titration calorimetry. Methods in Enzymology. 2016;**567**:181-213. DOI: 10.1016/bs.mie.2015.08.012

[12] Schmidtchen FP. Isothermal titration calorimetry in supramolecular chemistry. In: Anal. Methods Supramol. Chem. New Jersey, USA, Weinheim, Germany: John Wiley & Sons, Ltd, Wiley-VCH Verlag GmbH & Co. KGaA; 2006. pp. 55-78. DOI: 10.1002/ 9783527610273.ch3

[13] Biedermann F, Schneider H, Platz H. Experimental binding energies in supramolecular complexes. Chemical Reviews. 2016;**116**:5216-5300. DOI: 10.1021/acs.chemrev.5b00583

[14] Wang Y, Wang G, Moitessier N, Mittermaier AK. Enzyme kinetics by isothermal titration calorimetry: Allostery, inhibition, and dynamics.

Frontiers in Molecular Biosciences. 2020; **7**:1-19. DOI: 10.3389/fmolb.2020.583826

[15] Di Trani JM, Moitessier N, Mittermaier AK. Measuring rapid timescale reaction kinetics using isothermal titration calorimetry. Analytical Chemistry. 2017;**89**(13):7022-7030

[16] Di Trani JM, De Cesco S, O'Leary R, Plescia J, Do Nascimento CJ, Moitessier N, et al. Rapid measurement of inhibitor binding kinetics by isothermal titration calorimetry. Nature Communications. 2018;**9**(1):893. DOI: 10.1038/s41467-018-03263-3

[17] Braissant O, Wirz D, Göpfert B, Daniels AU. Use of isothermal microcalorimetry to monitor microbial activities. FEMS Microbiology Letters. 2010;**303**(1):1-8. DOI: 10.1111/ j.1574-6968.2009.01819.x

[18] Medina S, Raviv M, Saadi I, Laor Y. Methodological aspects of microcalorimetry used to assess the dynamics of microbial activity during composting. Bioresource Technology. 2009;**100**(20):4814-4820. DOI: 10.1016/ j.biortech.2009.05.015

[19] Gaisford S. Stability assessment of pharmaceuticals and biopharmaceuticals by isothermal calorimetry. Current Pharmaceutical Biotechnology. 2005; **6**(3):181-191. DOI: 10.2174/ 1389201054022913

[20] Kimura T, Matsushita T, Kamiyama T. Enthalpies of solution of aliphatic compounds in dimethyl sulfoxide. Thermochimica Acta. 2004; **416**(1–2):129-134. DOI: 10.1016/j.tca. 2003.02.001

[21] Sugiura T, Ogawa H. Thermodynamic properties of solvation of aromatic compounds in cyclohexane, heptane, benzene, 1,4-dioxane, and

chloroform at 298.15 K. The Journal of Chemical Thermodynamics. 2009; **41**(11):1297-1302. DOI: 10.1016/j.jct. 2009.06.001

[22] Korolev VP, Antonova OA, Smirnova NL, Kustov AV. Thermal properties of tetraalkylammonium bromides in several solvents. Journal of Thermal Analysis and Calorimetry. 2011; **103**(1):401-407. DOI: 10.1007/s10973- 009-0639-6

[23] Riveros DC, Martínez F, Vargas EF. Enthalpies of solution of methylcalix[4] resorcinarene in non-aqueous solvents as a function of concentration and temperature. Thermochimica Acta. 2012; **548**:13-16. DOI: 10.1016/j.tca.2012. 08.022

[24] Alves N, Bai G, Bastos M. Enthalpies of solution of paracetamol and sodium diclofenac in phosphate buffer and in DMSO at 298.15 K. Thermochimica Acta. 2006;**441**(1):16-19. DOI: 10.1016/ j.tca.2005.11.031

[25] Badelin VG, Mezhevoi IN. The thermochemical characteristics of solution of L-cysteine and L-asparagine in aqueous 1,4-dioxane and acetone. Russian Journal of Physical Chemistry A. 2009;**83**(7):1121-1124. DOI: 10.1134/ S0036024409070139

[26] Mezhevoi IN, Badelin VG. Standard enthalpies of dissolution of L-alanine in the water solutions of glycerol, ethylene glycol, and 1,2-propylene glycol at 298.15 K. Russ. The Journal of Physical Chemistry. A. 2010;**84**(4):607-610. DOI: 10.1134/S0036024410040151

[27] Badelin VG, Smirnov VI. The dependence of the enthalpy of solution of L-methionine on the composition of water-alcohol binary solvents at 298.15 K. Russ. The Journal of Physical

*Isothermal Calorimetry: Molecular Interactions between Small Molecules in Organic Solvents DOI: http://dx.doi.org/10.5772/intechopen.104756*

Chemistry. A. 2010;**84**(7):1163-1168. DOI: 10.1134/S0036024410070150

[28] Smirnov VI, Badelin VG. Enthalpies of solution of glycylglycine and diglycylglycine in aqueous alcohols at 298.15 K. Thermochimica Acta. 2008;**471** (1–2):97-99. DOI: 10.1016/j.tca.2008. 01.020

[29] Guseva GB, Antina EV, Berezin MB, V'yugin AI. Enthalpies of dissolution of cobalt(II) and copper(H) acetylacetonates in organic solvents. Russ. The Journal of Physical Chemistry. 2005;**79**(6):908-911

[30] Vandyshev VN, Ledenkov SF. Thermochemical study of complexation and solvation in β-alanine- Cu(NO3) 2-water-ethanol system: 1. Enthalpy characteristics of complexing copper(II) ion with β-alanine. Russian Journal of Inorganic Chemistry. 2011;**56**(3): 479-483. DOI: 10.1134/S0036023 611030259

[31] Jóźwiak M, Madej-Kiełbik L. Effect of temperature on the process of complex formation crown ether 15C5 with Na+ in the (water + ethanol) mixture at temperatures from (293.15 to 308.15) K. Thermochimica Acta. 2014; **580**:13-19. DOI: 10.1016/j. tca.2014.01.019

[32] Usacheva TR, Pham Thi L, Terekhova IV, Kumeev RS, Sharnin VA. Application of isothermal titration calorimetry for evaluation of wateracetone and water-dimethylsulfoxide solvent influence on the molecular complex formation between 18-crown-6 and triglycine at 298.15 K. Journal of Thermal Analysis and Calorimetry. 2015; **121**(3):975-981. DOI: 10.1007/ s10973-015-4630-0

[33] Bernaczek K, Mielańczyk A, Grzywna ZJ, Neugebauer D. Interactions between fluorescein isothiocyanate and star-shaped polymer carriers studied by isothermal titration calorimetry (ITC). Thermochimica Acta. 2016;**641**:8-13. DOI: 10.1016/j.tca.2016.08.007

[34] Arnaud A, Bouteiller L. Isothermal titration calorimetry of supramolecular polymers. Langmuir. 2004;**20**(16): 6858-6863. DOI: 10.1021/la049365d

[35] Wadsö I, Goldberg RN. Standards in isothermal microcalorimetry: (IUPAC technical report). Pure and Applied Chemistry. 2001;**73**(10):1625-1639. DOI: 10.1351/pac200173101625

[36] TA Instruments, Nano ITC. Available from: https://www.tainstrume nts.com/nano-itc/ [Accessed: January 15, 2022]

[37] Velázquez-Campoy A, Ohtaka H, Nezami A, Muzammil S, Freire E. Isothermal titration calorimetry. Current Protocols in Cell Biology. 2004;**23**:1-24. DOI: 10.1002/0471143030.cb1708s23

[38] Quin CF. Advanced ITC techniques: The reverse titration. 2013. http://www. tainstruments.com/pdf/literature/ MCAPN-2013-01\_Advanced\_ITC\_Tech niques-The\_Reverse\_Titration.pdf

[39] Spuches AM, Kruszyna HG, Rich AM, Wilcox DE. Thermodynamics of the as(III)-thiol interaction: Arsenite and monomethylarsenite complexes with glutathione, dihydrolipoic acid, and other thiol ligands. Inorganic Chemistry. 2005;**44**(8):2964-2972. DOI: 10.1021/ ic048694q

[40] Wiseman T, Williston S, Brandts JF, Lin LN. Rapid measurement of binding constants and heats of binding using a new titration calorimeter. Analytical Biochemistry. 1989;**179**(1):131-137. DOI: 10.1016/0003-2697(89)90213-3

[41] Myszka DG, Abdiche YN, Arisaka F, Byron O, Eisenstein E, Hensley P, et al. The ABRF-MIRG'02 study: Assembly state, thermodynamic, and kinetic analysis of an enzyme/inhibitor interaction. Journal of Biomolecular Techniques. 2003;**14**(4):247-269

[42] Bastos M, Velazquez-Campoy A. Isothermal titration calorimetry (ITC): A standard operating procedure (SOP). European Biophysics Journal. 2021;**50** (3–4):363-371. DOI: 10.1007/s00249- 021-01509-5

[43] Campoy AV, Claro B, Abian O, Höring J, Bourlon L. A multi-laboratory benchmark study of isothermal titration calorimetry (ITC) using Ca2+ and Mg2+ binding to EDTA. European Biophysics Journal. 2021;**50**:429-451. DOI: 10.1007/s00249-021-01523-7

[44] Fisher HF, Singh N. Calorimetric determination of binding constants. Energy Biological Macromolecules. 1995; **259**:194-221. DOI: 10.1016/0076-6879 (95)59045-5

[45] Horn JR, Russell D, Lewis EA, Murphy KP. Van't Hoff and calorimetric enthalpies from isothermal titration calorimetry: Are there significant discrepancies? Biochemistry. 2001; **40**(6):1774-1778. DOI: 10.1021/ bi002408e

[46] Mizoue LS, Tellinghuisen J. Calorimetric vs. van't Hoff binding enthalpies from isothermal titration calorimetry: Ba2+-crown ether complexation. Biophysical Chemistry. 2004;**110**(1–2):15-24. DOI: 10.1016/j. bpc.2003.12.011

[47] Markova N, Hallén D. The development of a continuous isothermal titration calorimetric method for equilibrium studies. Analytical

Biochemistry. 2004;**331**(1):77-88. DOI: 10.1016/j.ab.2004.03.022

[48] Degassing station. https://www.ta instruments.com/degassing-station/. [Accessed: February 4, 2022]

[49] Herrera I, Winnik MA. Differential binding models for direct and reverse isothermal titration calorimetry. The Journal of Physical Chemistry. A. 2016; **120**:2077-2086. DOI: 10.1021/acs. jpcb.5b09202

[50] Piñeiro Á, Muñoz E, Sabín J, Costas M, Bastos M, Velázquez-Campoy A, et al. AFFINImeter: A software to analyze molecular recognition processes from experimental data. Analytical Biochemistry. 2019;**577**:117-134. DOI: 10.1016/j.ab.2019.02.031

[51] Duvvuri H, Wheeler LC, Harms MJ. Pytc: Open-source python software for global analyses of isothermal titration calorimetry data. Biochemistry. 2018; **57**(18):2578-2583. DOI: 10.1021/acs. biochem.7b01264

[52] Gans P, Sabatini A, Vacca A. Simultaneous calculation of equilibrium constants and standard formation enthalpies from calorimetric data for systems with multiple equilibria in solution. Journal of Solution Chemistry. 2008;**37**(4):467-476. DOI: 10.1007/ s10953-008-9246-6

[53] Le VH, Buscaglia R, Chaires JB, Lewis EA. Modeling complex equilibria in isothermal titration calorimetry experiments: Thermodynamic parameters estimation for a threebinding-site model. Analytical Biochemistry. 2013;**434**(2):233-241. DOI: 10.1016/j.ab.2012.11.030

[54] Zhao H, Piszczek G, Schuck P. SEDPHAT – A platform for global ITC analysis and global multi-method

*Isothermal Calorimetry: Molecular Interactions between Small Molecules in Organic Solvents DOI: http://dx.doi.org/10.5772/intechopen.104756*

analysis of molecular interactions. Methods. 2015;**76**:137-148. DOI: 10.1016/ j.ymeth.2014.11.012

[55] Nicholls IA. Thermodynamic considerations for the design of and ligand recognition by molecularly imprinted polymers. Chemistry Letters. 1995;**24**(11):1035-1036. DOI: 10.1246/ cl.1995.1035

[56] Scheuermann TH, Brautigam CA. High-precision, automated integration of multiple isothermal titration calorimetric thermograms: New features of NITPIC. Methods. 2015;**76**:87-98. DOI: 10.1016/j. ymeth.2014.11.024

[57] Dogadkin DN, Soboleva IV, Kuz'min MG. Formation enthalpy and entropy of exciplexes with variable extent of charge transfer in solvents of different polarity. High Energy Chemistry. 2001;**35**(4):251-257. DOI: 10.1023/A:1017636612185

[58] Solís C, Grosso V, Faggioli N, Cosa G, Romero M, Previtali C, et al. Estimation of the solvent reorganization energy and the absolute energy of solvation of charge-transfer states from their emission spectra. Photochemical & Photobiological Sciences. 2010;**9**(5):675-686. DOI: 10.1039/b9pp00190e

[59] Nestora S, Merlier F, Beyazit S, Prost E, Duma L, Baril B, et al. Plastic antibodies for cosmetics: Molecularly imprinted polymers scavenge precursors of malodors. Angewandte Chemie, International Edition. 2016;**55**:6252-6256. DOI: 10.1002/anie.201602076

[60] Panagiotopoulou M, Salinas Y, Beyazit S, Kunath S, Duma L, Prost E, et al. Molecularly imprinted polymer coated quantum dots for multiplexed cell targeting and imaging. Angewandte Chemie, International Edition. 2016;

**55**(29):8244-8248. DOI: 10.1002/anie. 201601122

[61] Job P. Formation and stability of inorganic complexes in solution. Annali di Chimica. 1928;**9**:113-203

[62] Renny JS, Tomasevich LL, Tallmadge EH, Collum DB. Method of continuous variations: Applications of job plots to the study of molecular associations in organometallic chemistry. Angewandte Chemie, International Edition. 2013;**52**(46):11998-12013. DOI: 10.1002/anie.201304157

[63] Goldberg RN, Kishore N, Lennen RM. Thermodynamic quantities for the ionization reactions of buffers. Journal of Physical and Chemical Reference Data. 2002;**31**(2):231-370. DOI: 10.1063/1.1416902

[64] Wang HM, Loganathan D, Linhardt RJ. Determination of the pK(a) of glucuronic acid and the carboxy groups of heparin by 13Cnuclear-magnetic-resonance spectroscopy. The Biochemical Journal. 1991;**278**(3):689-695. DOI: 10.1042/ bj2780689

#### **Chapter 2**

## A State of Art Review on Thermodynamics Performance Analysis in Pulse Detonation Combustor

*Pinku Debnath and Krishna Murari Pandey*

#### **Abstract**

Pulse detonation engines (PDEs) are most exciting for future propulsion generation. Detonation combustion in pulse detonation combustor is an energetic combustion process which is differs from other combustion process. The detonation wave propagation in detonation tube is a pulse setting combustion phenomena. Detonation combustion process is thousands times faster than deflagration combustion process. PDE utilizes several pulse of detonation wave to produce propulsive force. The potential applications of PDEs are drastically reduces the cost of orbit transfer vehicle system and flying mode applications. Of course it can be used as ground level applications also. Draw back are DDT in shortest possible time in the combustor. In this regards, worldwide researchers are focusing on scientific and technical issues related to improvement of PDC. The present chapter deals with review study on detonation combustion process, historical overview on chemical kinetics, calorimetric and entropy transport, energy and exergy analysis and factor effecting on deflagration to detonation transition with recommendable future research.

**Keywords:** pulse detonation combustor, CFD, detonation, deflagration, ejector

#### **1. Introduction**

The high speed engine concept was born in the early 1900s, which produced shaft work and designed to drive a variety of vehicles, including ships and locomotives, until further introduction of jet engine on 1930s. The history of pulse detonation engines concept can be traced back to German engineer Hoffmann, H. [1]. In 1941, they tested a prototype engine using acetylene-oxygen and benzene-oxygen mixtures. Earlier in between 1952 and 1956, Nicholls et al. [2] at the University of Michigan have independently come up with the idea of using intermittent detonation for propulsion system and built the first PDE, which utilized detonation of hydrogen-air mixture to produce thrust. When crude oil prices increased significantly in the mid-1980s Eidelman et al. reinitiated research on PDE to overcome these scarcities. Krzycki [3] experimentally studied the propane-air pulse detonation engine which is operating

at 25–50 Hz at naval ordnance test station at China Lake, California on 1962. From this analysis it was observed that minimum thrust is produce at minimum operating frequency. Lynch et al. [4] performed CFD studies on PDREs and air breathing PDEs on 1990. Their analyses forecasted that PDEs would be incorporated with space transport vehicles by the early 2000s. Earlier starting in 1990s, experimental study of single and multi-pulse detonation engine combustor were conducted by Bussing et al. [5] at Adroit systems Inc., a company that was bought up by Pratt and Whitney on 2001. They tested pulse detonation engine using different fuels including hydrogen and ethylene. Most of the PDE research centers are found in Canada, US, Russia, China. There are very less number of PDE research centre in India, so researcher are focusing on pulse detonation research area as it is excited for future propulsion technology [6].

Propulsion applications of detonation can be classified into three categories: standing detonation, pulse detonation and rotating detonation [7–9]. The basic pulse detonation engine has a very simple structure. It consists of a constant area tube. The deflagration to detonation transition is controlled by supplying fuel and oxidizer in detonation tube. The ignition system and nozzle are used for accelerating the flow, which is to be used for propulsion. A practical pulse detonation engine may also have one or more devices to bring about deflagration to detonation transition such device are Shchelkin spiral and blockage [10]. The PDE consists of two or more combustion chamber, which is joined to common plenum chamber. The conditions are applied for accelerating the flow before entering the detonation tube with different nozzle. The can-annular four chambered PDE is can illustrated for propulsion system. In multichambered design, each chamber can be at different stage in the cycle, thus creating a smoother flow through the nozzle [11]. The ejector enhances the deflagration to detonation transition in detonation tube with an array setting in exit section of pulse detonation rocket engine. Another feature of the ejector design is that the detonation waves from the combustors can be used to enhance the propulsive performance, which provides additional thrust enhancement [12].

Detonation is a supersonic mode of combustion process. In combustion process detonation waves are much more energetic process than conventional deflagration combustion and it produces a very strong wave coupled with a chemical reaction zone, propagating at supersonic speed. A detonation wave compresses combustion mixture, increasing the combustion product pressure, density of species mass fraction. It is a subsonic combustion process and fuel air reaction propagates at relatively low speed and reasonably low pressure from a trailing reaction zone. The propagation of deflagration mode of combustion consists of diffusion of unburned gases ahead of flame front and burnt gases behind the combustion flame. Deflagration produces small decrease in pressure and can be modeled as a constant pressure process [13]. One of the primary attributes deflagration flame travels at a speed, which is significantly lower than that the speed of sound (Ma < 1). So it can be identify by subsonic combustion process. Detonation combustion is a constant volume combustion process. The strength of leading shock depends on the detonation wave propagation velocity. A simple planar model for the supersonic detonation shock wave is used for Chapman-Jouguet detonation model analysis [14]. This is a rapid exothermic reaction and instantaneously changes the local pressure and temperature. The ignition of fuelair mixture can produce deflagration flame and later on transition to detonation wave. The different combustor geometry can accelerate the deflagration flame and transition to detonation wave. Several researchers have been studied on PDE with research gape and scope of future research work [15–17]. The applications of RDE chamber

#### *A State of Art Review on Thermodynamics Performance Analysis in Pulse Detonation Combustor DOI: http://dx.doi.org/10.5772/intechopen.103005*

are jet engines, such as turbojet or gas turbine, ramjet or rocket. Continuous detonation wave engine is used for supersonic and hypersonic propulsion applications. The framework of French Research and Development and scientific research also consider these for space applications. However incoming reacting air-mixture is greater than the C-J velocity of fuel-air mixture. Such engine is scramjet engine with an oblique detonation wave at inlet to combustor called the oblique detonation wave engine [18]. Chapman [19] explained on 1889 that the minimum speed of burnt gas is equal to speed of sound in gas mixture. Later on Jouguet on 1905 [20] applied Hugoniot's method to explain the detonation velocity. The explosive mixture can get supported with two modes of combustion. When the flame propagates at slow velocity relative to unburnt gases, it is define as deflagration mode. In detonation mode wave propagates at about 2000 m/s accompanied by an overpressure rise is near about 20 bars [21–23]. They independently developed the basic thermodynamic model behind detonation. Principle operation of standing detonation engine is relatively simple. Fuel is injected into supersonic flow and detonation wave is stabilized inside the engine by wedge or other means and products are expanding inside nozzle. The combustion wave velocity can be propagates at higher the C-J detonation velocity within the Mach number of 5. The principle of rotating detonation engine (RDE) is based on the formation of detonation in a disk type combustion chamber. The shape of combustion chamber is toroidal or ring-like shape [24, 25]. The detonation wave parameters are depending on critical detonation tube diameter and minimum detonation tube diameter. The minimum and critical diameters are important parameters for evaluation of performance of PDE. The detonation will successfully propagate in a tube when the diameter must be larger than λ/3, where λ is the cell size. For square and rectangular ducts, the width and height of the duct must be larger than λ [26]. A review of the gas dynamics and chemistry of real detonation is discovered by Fickett and Davis [27]. They found out initiation of detonation wave, which follows by a series of percentage of fuel-oxidizer mixture in combustion chamber. The detonation wave in a confined tube causes the reaction of fuel-air mixtures, which creates turbulence; as a result "an explosion in an explosion" is takes placed. The two strong shock waves are created in the opposite direction, the forward shock waves are known as retonation. A self-propagating C-J detonation wave is formed at steady state retonation process. The pre-detonation wave velocity is 1000 m/s while the characteristic C-J detonation speed is over 2000 m/s. A large explosion occurs at onset of detonation, resulting in an over-driven detonation wave that decays to the C-J velocity. The wall roughness controls the wave propagation by inducing large-amplitude unsteady and turbulent flow, complex wave interaction processes and high temperature behind shock reflections. These effects represent ways that the flow can generate large-scale turbulence for flame folding and large temperature fluctuations causing detonation initiation [28].

#### **2. Review on thermodynamics cycle analysis**

A pulse detonation engine uses repetitive cycle of detonation waves to combust fuel-oxidizer mixture for producing thrust. PDE operates by propagating detonation wave through a tube filled with a combustible mixture and generates propulsive thrust. This process results are near about constant volume combustion process, which produces high pressures from the leading shock wave. Pulse detonation engine consists of valve less combustor with straight tube, which is closed at one end and open at other end. The pulse detonation engine combustion cycle consist of four basic thermodynamics process. The first process is filling time (*tfill*) of fuel-air mixture for detonation combustion, which is estimated as length of the tube over filling velocity. The second one is detonation combustion. This process takes place within fraction of millisecond. As soon as the detonation wave reaches to the closed end region the pressure and velocity decrease from initial position to exit end region. Fully developed detonation wave travels with the magnitude of Chapman-Jouguet speed. This C-J speed of reacting fuel-air mixture varies between 1400 m/s and 1800 m/s. The detonation time of the wave (*tc*) is therefore similarly estimated by the length of the tube over the C-J wave velocity. The time required for blow down (*tb*) stage can be estimated by the length of the tube to the rarefaction velocity. At last in the purging process the tube is scavenged off hot detonation products by using fresh air. Purging process is necessary to prevent auto ignition of the fresh fuel-oxidizer mixture. The time taken for purging the tube with the fresh air (*tpurge*) is the length of the tube over the purging velocity [29]. So total sum of the time for all the four stages are as follows:

$$\mathbf{T} = \mathbf{t}\_{f\&l} + \mathbf{t}\_c + \mathbf{t}\_b + \mathbf{t}\_{\text{purge}} \tag{1}$$

The PDE can run by any fuel, liquid or gaseous, like natural gas, propane, bio-gas, hydrogen, kerosene, jet fuels and octane etc. From an engineering stand point fuel can be selected for based on heating value, detonability, ignition time, energy release, adiabatic flame temperature and sensivity with air [30]. Povinelli and Yungster [31] studied the thermodynamic cycle of hydrogen-air mixture at static conditions in pulse detonation combustor. The specific thrust, fuel consumptions and impulse of detonation combustion are analyzed by using CFD analysis with finite rate chemistry. Alam et al. [32] studied on Brayton, Humphery and ideal thermodynamics cycle analysis in pulse detonation combustor. They found Humphery cycle efficiency can be increases with higher value of compression ratio. The thermodynamics cycle efficiency of air breathing pulse detonation engine is studied by Wu et al. [33]. They found that chocked convergent-divergent nozzle is required to improve the efficiency. Vutthivithayarak et al. [34] discussed the Humphrey and F. J. (Fickett-Jacobs) cycles in PDE. These cycles are illustrated with hydrogen-air combustion for generic heat release.

#### **3. Review on chemical kinetics and entropy transport**

The two-step chemical kinetics model of detonation combustion has been studied by Fomin [35]. This kinetics model has been used for stoichiometric, lean and rich mixture for combustion. This model is also followed by Le Chatelier's principle and 2nd law of thermodynamics. The pulse detonation combustor has lower entropy change and self-pressure gain compared to isobaric combustion process for same operating conditions [36]. Mehdi Safari et al. [37] studied on entropy generation with species transport equation for detonation combustion by large eddy simulation. Detonation initiation in hydrogen-air depends on mixture sensivity and geometrical parameters. Qi et al. [38] investigated the thermodynamics characteristics of methane-air detonation in pressure gain combustor. They compared the entropy change in detonation combustion process with gas turbine cycle. They found that cycle efficiency enhance rate up to 11.89%. Lu et al. [39] studied on DDT in a channels with obstacles using chemical diffusive model (CDM) integrated with reactive

#### *A State of Art Review on Thermodynamics Performance Analysis in Pulse Detonation Combustor DOI: http://dx.doi.org/10.5772/intechopen.103005*

Navier stokes equation. They found that CDM reduces the ignition time of detonation wave. Wu et al. [40] studied on atomization of liquid fuel detonation combustion. They found that nozzle can effectively atomize fuel-air mixtures under high pressure condition. Maciel and Marques [41] studied on hydrogen fuelled single cycle pulse detonation engine in Ansys Fluent. When OH\* kinetics added to the reaction set, they found cellular structure of detonation wave front in reaction zone. Ivanov et al. [42] studied on hydrogen-oxygen flame acceleration and transition from DDT in a channel using reactive Navier-Stokes equations. They found that steady detonation wave front is form in wider detonation channels of 10 mm and closed to C-J detonation propagation speed. Srihari et al. [43] studied on stoichiometric ethylene-air mixture of detonation combustion with one-step overall reaction model. They found that chemical reaction models have capable to predict the detonation wave velocity with reasonable accuracy.

#### **4. Review on energy and exergy analysis**

Ma et al. [44] studied on temporal variation of activation energy release rate of iso-octane vapor-air mixture in an obstacle-filled detonation tube. Their result shows that the activation energy influences the flame propagation parameters and deflagration-to-detonation transition process. Hutchins and Metghalchi [45] studied on exergy analysis of pulse detonation engine. They found that during deflagration to detonation transition period exergy loss is more. Bellini and Lu [46] studied on exergy analysis of fuel-air mixture at high frequency source within the detonation chamber. They found that combustion product accelerates inside the combustor in presence of Shchelkin spiral. The exergy analysis of pulse detonation power device designed for power production using gaseous fuel methane (CH4) and propane (C3H8) is analysis by Bellini and Lu [47]. The exergetic efficiency was analyzed for different cycle frequency corresponding to detonation tube length. Rouboa et al. [48] studied on exergy loss of hydrogen-air detonation during shock. They also observed exergy destruction increase with augmentation of hydrogen concentration in reacting mixture. Petela [49] studied on exergy analysis of gaseous fuel-air detonation. They observed that exergy gives a quantitative theoretical useful work that is obtained from different energy form combustion process and it is a function of system and environment. Som and Datta [50] and Som and Sharma [51] studied on theoretical model of energy and exergy balance in a spray combustion process. They found that exergy destruction in this combustion process can be reduced through proper control of chemical reactions.

#### **5. Results from CFD simulation and calorimetric analysis**

The numerical investigations have been done in Ansys fluent platform. The **Figure 1** shows that the ejector effects on unsteady detonation combustion wave phenomena in pulse detonation combustor. The time dependent detonation wave contour plots clearly shows that 0.033 seconds is required to reach the fully developed detonation wave [52]. They also found that ejector plays the vital role for vortex formation of reacting mixture in PDE combustor. They also observed that leading vortex rings are found in shrouded ejector taper angle of +4°. The **Figure 2**

#### **Figure 1.**

*Effect of shrouded ejectors on vortex ring formation of detonation wave [53].*

#### **Figure 2.**

*The mass fraction contour analysis of NOx pollutant number of (a) hydrogen-air and (b) kerosene-air combustion [53].*

*A State of Art Review on Thermodynamics Performance Analysis in Pulse Detonation Combustor DOI: http://dx.doi.org/10.5772/intechopen.103005*

#### **Figure 3.**

*Deflagration and detonation zone define by C-J velocity for exergy analysis [54].*

shows the mass fraction contour analysis of NOx pollutant number of hydrogen-air and kerosene-air detonation [53]. Lesser the fuel mass fraction higher the exergetic efficiency was found in pulse detonation combustor. The **Figure 3** shows the deflagration and detonation control volume for exergy analysis [54]. Alam et al. [55, 56] numerically studied on hydrogen-air detonation in pulse detonation combustor. Later on they also studied detonation combustion using alternative fuels, i.e. octane (C8H18), hexane (C6H14), pentane (C5H12)-air combustion in PDE combustor. They observed combustion efficiency of pentane-air mixture is higher than that of other fuels. Alam et al. [57] studied the combustion wave propagation in obstructed detonation tube. Their simulation results were carried out for stoichiometric mixture of kerosene-air and butane-air mixture at atmospheric conditions. They found that mixing of butane-air combustion process is better than kerosene-air mixture. Furthermore, the stoichiometric ethane-air (C2H6-air) and ethylene-air (C2H4-air) fuel mixture at atmospheric pressure conditions has been studied by Alam et al. [58]. The effect of blockage ratio of 0.4, 0.5, 0.6 and 0.7 in channel for detonation wave acceleration are shown in **Figure 4**. The contour plot analysis shows the shock wave initiation and propagation time period in detonation tube is reduced by smaller blockage ratio of 0.5 [59]. Tripathi et al. [60] computationally studied on effect of obstacle on flame propagation velocity. Alam et al. [61] studies on flame acceleration in pulse detonation engine with changing the obstacle clearance. They found that combustor pressure is reduced as increase the obstacle clearance. P. Debnath and Pandey [62] studied on deflagration to detonation transition in PDE combustor with Schelkin spiral effect inside the detonation tube. They found that Schelkin spiral accelerate the flame propagation. Alam et al. [63] numerically studied on flame propagation in obstructed pulse detonation combustor with hydrogen-air mixture. They found that performance is increase up to 4.46% and this value increase for = 1.3. Debnath and Pandey [64] studied on effect of different nozzle on flame acceleration and they found that divergent nozzle has more effect on flame acceleration. The **Figure 5** shows the comparison of thrust power for PDE combustor with several nozzle. Chourasia et al. [65] studied on progress and motivation of research in pulse detonation combustor. Xudong Zhang et al. [66] studied on critical mode

#### *Applications of Calorimetry*

**Figure 4.** *Effect of blockage ratio on detonation wave propagation [59].*

#### **Figure 5.**

*The propulsive thrust variation for PDE combustor with C-D nozzle, C-nozzle, D-nozzle and without nozzle at different Mach number [64].*

of gaseous methane-air detonation propagation in an annular tube based on reactive Navier Stokes Equations. They found that trajectories of triple point of the shock wave cell structures are petal pattern. Wang et al. [67] studied on effect of oxygen concentration on propane-air detonation in pulse detonation engine with straight nozzle, convergent nozzle, and convergent-divergent (CD) nozzle. Their results indicate that for the PDE with straight nozzle requires the shortest possible time for reacting gas burnt with high-temperature in detonation tube. Jishnu Chandran and Salih [68] studied on development of a benchmark solution in compressible liquid flows for shock tube problems. The compressibility effects in liquid water have been studied using the high-accuracy modified NASG equation of state. Arjun Singh et al. [69] studied on thermodynamic parameters for the formation of activation energy

*A State of Art Review on Thermodynamics Performance Analysis in Pulse Detonation Combustor DOI: http://dx.doi.org/10.5772/intechopen.103005*

and self-acceleration for thermal explosion from critical temperature. They found that the thermal stability has been significantly reduced in presence of hydroxylterminated poly butadiene. Dong et al. [70] studied on correlations among detonation velocity, thermal stability, heat of combustion and decomposition kinetics of nitric esters. They found that oxygen coefficient plays positive role on decomposition of heat release efficiency of detonation combustion.

#### **6. Concluding remarks**

The above literature survey represents that there is more research is needed in pulse detonation combustor for shortest possible pulse time of deflagration to detonation transition. The future proposed research can be analyzed by changing the design of PDE combustor and operating conditions. The series of numerical simulations and optimization can be performed desire research objectives of pulse detonation engine. From the CFD and calorimetric analysis the smaller blockage ratio of 0.5 is found better to reduce detonation wave run up distance. The ejector enhance the shortest possible time of 0.033 s, which is required for fully developed detonation wave. More possible pulse time can be reduced by ejector geometry modification. Lesser the hydrogen fuel mass fraction of 0.25 higher the exergetic efficiency of 67.55% is obtained from detonation combustion process. Once the computational model is validated, further simulation can be carried out with accuracy. There are several detonation tube geometry is steal in debate for acoustics atomization and evaporative characteristics of liquid fuel detonation wave.

#### **Author details**

Pinku Debnath1 \* and Krishna Murari Pandey2

1 Department of Mechanical Engineering, National Institute of Technology Agartala, Agartala, Tripura, India

2 Department of Mechanical Engineering, National Institute of Technology Silchar, Silchar, Assam, India

\*Address all correspondence to: er.pinkunits@yahoo.com

© 2022 The Author(s). Licensee IntechOpen. This chapter is distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/3.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

### **References**

[1] Hoffmann H. Reaction-propulsion produced by intermittent detonative combustion. German Research Institute for Gliding, Report ATI-52365. 1940

[2] Nicholls JA, Wilkinson HR, Morrison RB. Intermittent detonation as a thrust-producing mechanism. Journal of Jet Propulsion. 2012;**27**(5):534-541

[3] Krzycki LJ. "Performance characteristics of an intermittent detonation device," NAVWEPS, Report 7655, Naval ordnance test station. China Lake, California: 1962

[4] Lynch ED, Edelman R, Palaniswamy S. Computational fluid dynamic analysis of the pulse detonation wave engine concept. In: 32nd Aerospace Sciences Meeting & Exhibit; January 10-13; Reno, NV. 2012

[5] Bussing T, Bratkovich TE, Hinkey JB. Practical implementation of pulse detonation engines, AIAA-1997- 2748. In: 33rd AIAA/ASME/SAE/ASEE Joint Propulsion Conference and Exhibit; July 6-9; Seattle, WA. 2012

[6] Lu FK. Pulse Detonation Propulsion Systems Introduction. Texas: The University of Texas Arlington; 2007

[7] Bussing T, Pappas G. An introduction to pulse detonation engines. In: 32nd Aerospace Sciences Meeting and Exhibit; Reno, NV. 2012

[8] Eidelman S, Grossmann W, Lottati I. Review of propulsion applications and numerical simulations of the pulsed detonation engine concept. Journal of Propulsion and Power. 1991;**7**(6):857-865

[9] Li JL, Fan W, Qiu H, Yan C, Wang YQ. Preliminary study of a pulse detonation

wave engine. Aerospace Science and Technology. 2010;**14**:161-167

[10] Panicker PK. The development and testing of pulsed detonation engine ground demonstrators [Ph.D. thesis]. The University of Texas at Arlington; 2008

[11] Oates G. Aerothermodynamics of Gas Turbine and Rocket Propulsion. 3rd ed. American Institute of Aeronautics & Astronautics (AIAA). 1997

[12] Munipalli R, Shankar V, Wilson DR, Lu FK. Preliminary design of a pulse detonation based combined cycle engine. In: 15th International Symposium on Air Breathing Engines; 2-7 September; Bangalore, India. 2001

[13] Kuo KK. Principles of Combustion. 2nd ed. John Wiley & Sons; 2005. ISBN 10: 9780471046899

[14] Nekkanti K. Analysis of thrust development in a pulse detonation engine [M.S. thesis]. The University of Texas at Arlington. Libraries Research Commons; 2010

[15] Pandey KM, Debnath P. Reviews on recent advances in pulse detonation engines. Journal of Combustion. 2016;**2016**:1-16

[16] Debnath P, Pandey KM. Performance investigation on single phase pulse detonation engine using computational fluid dynamics. In: Proceedings of the ASME 2013 International Mechanical Engineering Congress & Exposition, IMECE2013; November 15-21; San Diego, CA. 2013

[17] Pandey KM, Kumar J. CFD analysis of pulse detonation engine: A review. Journal of Material Science and Mechanical Engineering (JMSME). 2016;**3**(2):111-116

#### *A State of Art Review on Thermodynamics Performance Analysis in Pulse Detonation Combustor DOI: http://dx.doi.org/10.5772/intechopen.103005*

[18] Kailasanath K. Application of detonations to propulsion: A review. In: 37th AIAA Aerospace Sciences Meeting and Exhibit; January 11-14; Reno, NV. 1999

[19] Rankine WJM. On the thermodynamic theory of finite longitudinal disturbance. Philosophical Transactions of the Royal Society of London. 1870;**160**:277-288

[20] Lee JHS. Dynamic parameters of gaseous detonations. Annual Review of Fluid Mechanics. 1984;**16**:311-336

[21] Fickett W, Davis WC. Detonation Theory and Experiment. India: Dover Publications Inc; 1979. ISBN 10: 0486414566

[22] Taylor G. The dynamics of the combustion products behind plane and spherical detonation fronts in explosives. Proceedings of the Royal Society of London A. 1950;**200**(1061):235-247

[23] Wintenberger E. Application of steady and unsteady detonation waves to propulsion [Ph.D. dissertation]. Pasadena, CA: California Institute of Technology; 2004

[24] Kindracki J, Wolanski P, Gut Z. Experimental research on the rotating detonation in gaseous fuels-oxygen mixture. Shock Waves. 2011;**21**(2):75-84

[25] Hugoniot PH. Memoire sur la propagation du movement dans les corps et plus specialement dans les gaz parfaits premiere partie. Journal de l'Ecole Polytechnique. 1887;**57**:3-97

[26] Lee JHS. The Detonation Phenomenon. Cambridge University Press; 2014

[27] Fickett W, Davis WC. Detonation: Theory and Experiment. Dover Publications; 2010

[28] Oran ES, Gamezo VN. Origins of the deflagration-to-detonation transition in gas-phase combustion. Combustion and Flame. 2007;**148**:4-47

[29] Mattingly JD. Elements of Propulsion: Gas Turbines and Rockets. Reston, VA: AIAA; 2006

[30] Alhussan K, Assad M, Penazkov O. Analysis of the actual thermodynamic cycle of the detonation engine. Applied Thermal Engineering. 2016;**107**:339-344

[31] Povinelli LA, Yungster S. Thermodynamic Cycle and CFD Analyses for Hydrogen Fueled Air-breathing Pulsed Detonation Engines. American Institute of Aeronautics and Astronautics. National Aeronautics and Space Administration, Glenn Research Center; 2002

[32] Alam N, Sharma KK, Pandey KM. Thermodynamic performance of pulse detonation engine: A Technical Report". In: Proceedings of International Conference on Sustainable Computing in Science, Technology and Management (SUSCOM); February 26-28, 2019; Amity University Rajasthan, Jaipur, India. 2019. DOI: 10.2139/ssrn.3355295

[33] Wu Y, Ma F, Yang V. System performance and thermodynamics cycle analysis of air breathing pulse detonation engine. Journal of Propulsion and Power. 2003;**18**(4):556-567

[34] Vutthivithayarak R, Braun EM, Lu FK. On thermodynamics cycle for detonation engine. In: 28th International Symposium on shock waves. Berlin, Heidelberg: Springer; 2012

[35] Fomin PA. Reduced chemical kinetic model of detonation combustion of one and multi fuel gaseous mixtures with air. AIP Conference Proceedings. 2018;**1939**:020012

[36] Wintenberger E, Sheperd JE. Thermodynamic cycle of analysis for propagating detonations. Journal of Propulsion and Power. 2006;**22**:694-698

[37] Mehdi S, Reza M, Sheikhi H, et al. Entropy transport equation in large Eddy simulation for exergy analysis of turbulent combustion systems. Entropy. 2010;**12**:434-444. DOI: 10.3390/e1230434

[38] Qi L, Wang Z, Zhao N, Dai Y, Zheng H, Meng Q. Investigation of the pressure gain characteristics and cycle performance in gas turbines based on Interstage bleeding rotating detonation combustion. Entropy. 2019;**21**:265

[39] Lu X, Kaplan CR, Oran ES. Predictions of flame acceleration, transition to detonation and detonation propagation using the chemicaldiffusive model. Combustion and Flame. 2022;**235**:111705

[40] Wu Y, Han Q, Yang G. Effects of an acoustic atomizer upon liquid-fueled detonation initiations in a detonation tube. Experimental Thermal and Fluid Science. 2019;**109**:109863

[41] Maciel EC, Marques CST. 2-D simulation with OH\* kinetics of a single-cycle pulse detonation engine. Journal of Applied Fluid Mechanics. 2019;**12**(4):1249-1263

[42] Ivanov MF, Kiverin AD, Liberman MA. Flame acceleration and DDT of hydrogen-oxygen gaseous mixtures in channels with no-slip walls. International Journal of Hydrogen Energy. 2011;**36**(13):7714-7727

[43] Srihari P, Mallesh MA, Sai Krishna G, Charyulu BVN, Reddy DN. Numerical study of pulse detonation engine with one-step overall reaction model. Defence Science Journal. 2015;**65**(4):265-271

[44] Hu M, Xia Z, Gao W, Zhuo C, Wang D. Numerical simulation of the deflagration to detonation transition of iso-octane vapor in an obstaclefilled tube. International Journal of Spray and Combustion Dynamics. 2018;**10**(3):244-259

[45] Hutchins TE, Metghalchi M. Energy and exergy analysis of the pulse detonation engine. Journal of Engineering for Gas Turbines and Power. 2003;**125**(4): 1075-1080

[46] Bellini R, Lu FK. Exergy analysis of a pulse detonation power device. In: Proceedings of the 10th Brazilian Congress of Thermal sciences and Engineering-ENCIT 2004, Braz. Soc. of Mechanical Sciences and Engineering-ABCM; November 29–December 03; Rio de Janerio, Brazil. 2004

[47] Bellini R, Frank KL. Exergy analysis of a hybrid pulse detonation power device. Journal of Propulsion and Power. 2010;**26**(4):875-878. DOI: 10.2514/1.4414

[48] Rouboa A, Silva V, Couto N. Exergy analysis in hydrogen-air detonation. Journal of Applied Mathematics. 2012;**2012**:16

[49] Petela R. Application of exergy analysis to the hydrodynamic theory of detonation in gases. Fuel Processing Technology. 2000;**67**:131-145

[50] Som SK, Datta A. Thermodynamic irreversibilities and exergy balance in combustion processes. Progress in Energy and Combustion Science. 2008;**34**:351-376

[51] Som SK, Sharma NY. Energy and exergy balance in the process of spray combustion in a gas turbine combustor. ASME Journal of Heat Transfer. 2002;**124**:828-836

*A State of Art Review on Thermodynamics Performance Analysis in Pulse Detonation Combustor DOI: http://dx.doi.org/10.5772/intechopen.103005*

[52] Debnath P, Pandey KM. Numerical investigation of detonation combustion wave in pulse detonation combustor with ejector. Journal of Applied Fluid Mechanics. 2017;**10**(2):725-733

[53] Debnath P, Pandey KM. Numerical analysis of detonation combustion wave in pulse detonation combustor with modified ejector with gaseous and liquid fuel mixture. Journal of Thermal Analysis and Calorimetry. 2021;**145**:3243-3254

[54] Debnath P, Pandey KM. Exergetic efficiency analysis of hydrogenair detonation in pulse detonation combustor using computational fluid dynamics. International Journal of Spray and Combustion Dynamics. 2016;**9**(1):44-54

[55] Alam N, Pandey KM, Sharma KK. Combustion characteristics of hydrogenair mixture in pulse detonation engines. Journal of Mechanical Science and Technology. 2019;**33**(5):1-7

[56] Alam N, Sharma KK, Pandey KM. Effects of various composition of fuel-air mixture on performance of pulse detonation engine. Combustion, Explosion, and Shock Waves. 2019;**55**:708-717

[57] Alam N, Pandey KM, Sharma KK. Numerical investigation of combustion wave propagation in obstructed channel of pulse detonation engine using kerosene and butane fuels. Journal of Applied Fluid Mechanics. 2019;**12**(3):883-890

[58] Alam N, Sharma KK, Pandey KM. Numerical investigation of combustion phenomena in pulse detonation engine with different fuels. AIP Conference Proceedings. 2018;**1966**:020015

[59] Debnath P, Pandey KM. Effect of blockage ratio on detonation flame

acceleration in pulse detonation combustor using CFD. Applied Mechanics and Materials. 2014;**656**:64- 71. DOI: 10.4028/www.scientific.net/ AMM. 656.64

[60] Tripathi S, Pandey KM, Randive P. Computational study on effect of obstacles in pulse detonation engine. International Journal of Engineering & Technology. 2018;**7**:113-117

[61] Alam N, Sharma KK, Pandey KM. Numerical investigation of flame propagation in pulse detonation engine with variation of obstacle clearance. Journal of Thermal Analysis and Calorimetry. 2021;**140**:2485-2495

[62] Debnath P, Pandey KM. Computational study of deflagration to detonation transition in pulse detonation engine using Shchelkin spiral. Applied Mechanics and Materials. 2015;**772**:136-140

[63] Alam N, Sharma KK,

Pandey KM. Numerical investigation of fame propagation and performance of obstructed pulse detonation engine with variation of hydrogen and air. Journal of the Brazilian Society of Mechanical Sciences and Engineering. 2019;**41**:502

[64] Debnath P, Pandey KM. Numerical investigation of detonation combustion wave propagation in pulse detonation combustor with nozzle. Advances in Aircraft and Spacecraft Science. 2020;**7**(3):187-202

[65] Chourasia A, Verma KA, Pandey KM. Review of computational work in pulse detonation engines. International Journal of Innovative Technology and Exploring Engineering. 2019;**8**:398-401

[66] Zhang X, Gui M, Pan Z, Zhang H. Numerical investigation on the rotating detonation critical mode for a methaneair mixture in an annular tube using

#### *Applications of Calorimetry*

reactive Navier-Stokes equations. Journal of Thermal Analysis and Calorimetry. 2021;**144**:2285-2293

[67] Wang Z, Wei L, Qin W, Liang Z, Zhang K. Oxygen concentration distribution in a pulse detonation engine with nozzle-ejector combinational structures. Proceeding of IMechE Part G: Journal of Aerospace Engineering. 2021;**235**:1-10

[68] Jishnu Chandran R, Salih A. Development of a benchmark solution in compressible liquid flows: Analytical solution to the water shock tube problem. Journal of Thermal Analysis and Calorimetry. 2021. [Published online: 08 June 2021]

[69] Singh A, Soni PK, Sarkar C, Mukherjee N. Thermal reactivity of aluminized polymer-bonded explosives based on non-isothermal thermogravimetry and calorimetry measurements. Journal of Thermal Analysis and Calorimetry. 2019;**136**:1021-1035

[70] Dong J, Yan Q-L, Liu P-J, He W, Qi X-F, Zeman S. The correlations among detonation velocity, heat of combustion, thermal stability and decomposition kinetics of nitric esters. Journal of Thermal Analysis and Calorimetry. 2018;**131**:1391-1403

#### **Chapter 3**

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

*Viorel Chihaia, Valentin Alexiev and Hasan S. AlMatrouk*

#### **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.

**Keywords:** 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–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.

#### *Assessment of the Heat Capacity by Thermodynamic Approach Based on Density Functional… DOI: http://dx.doi.org/10.5772/intechopen.104083*

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 *φ<sup>i</sup> r* !� � <sup>¼</sup> <sup>P</sup>*<sup>ϖ</sup> <sup>μ</sup>*¼<sup>1</sup>*cμiχμ <sup>r</sup>* !� � is developed in analytic or numeric basis sets *χμ <sup>r</sup>* !� �**,** 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μ<sup>i</sup>* 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 *φ<sup>i</sup> r* !� � ! *<sup>φ</sup>nk* ! *r* !� � <sup>¼</sup> *eik* ! *r* ! *u nk* ! *r* !� �, where *<sup>k</sup>* ! 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 <sup>∈</sup>*n 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 *<sup>n</sup>*, the energy levels <sup>∈</sup> *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 denser*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 nk* ! *r* !� � <sup>¼</sup> <sup>P</sup><sup>∞</sup> *G* ! ¼1 *c G* ! *i n*, *k* � �! *eiG* ! *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 *<sup>φ</sup>i k* ! *r* !� � <sup>¼</sup> *ei k* ! *r* ! *u nk* ! *r* !� � <sup>¼</sup>

P<sup>∞</sup> *G* ! ¼1 *c G* ! *i n*, *k* � �! *e i k*! þ*G* ! � �*<sup>r</sup>* ! . 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 normconserving 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]**.**

*Assessment of the Heat Capacity by Thermodynamic Approach Based on Density Functional… DOI: http://dx.doi.org/10.5772/intechopen.104083*

#### **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 correlationexchange 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

#### *Applications of Calorimetry*

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*Assessment of the Heat Capacity by Thermodynamic Approach Based on Density Functional… DOI: http://dx.doi.org/10.5772/intechopen.104083*

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 Å�<sup>1</sup> , 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 *V*<sup>0</sup> and the corresponding energy *U*<sup>0</sup> can be determined by a full optimization of the lattice parameters and atom fractional coordinates. Furthermore, the volume dependency of the unit cell *U V*ð Þ 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 *V0*, the bulk modulus *<sup>B</sup>* ¼ �*<sup>V</sup> <sup>∂</sup><sup>P</sup> ∂V* � � *<sup>T</sup>* and its derivative *<sup>B</sup>*<sup>0</sup> <sup>¼</sup> *<sup>∂</sup><sup>K</sup> ∂P* � � *<sup>T</sup>* by

$$U(V) = U\_0 + \frac{B\_0 V}{B'} \left[ \frac{(V\_0/V)^{B'}}{B'-1} + 1 \right] - \frac{V\_0 B\_0}{B'-1} \tag{1}$$

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

*Assessment of the Heat Capacity by Thermodynamic Approach Based on Density Functional… DOI: http://dx.doi.org/10.5772/intechopen.104083*

**Figure 1**. The obtained values of the parameters *U0*, *B*<sup>0</sup> and *B*<sup>0</sup> 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]).

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

$$P(V) = -\frac{dE}{dV} = \frac{B\_0}{B'} \left[ \left(\frac{V\_0}{V}\right)^{B'} - 1 \right] \tag{2}$$

#### **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.*


#### **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.*

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

$$U = U\_o + U\_1 + U\_2 + U\_3 + \dots \tag{3}$$

where:

*U*0- is the static energy of the system in the equilibrium geometry,

*<sup>U</sup>*<sup>1</sup> <sup>¼</sup> <sup>P</sup> *mIα*Φ*<sup>α</sup> mIu<sup>α</sup>*,*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\_{m \text{la}} \sum\_{l \parallel \beta} \Phi\_{m \text{I}, l\text{J}}^{a\beta} u\_{a, m \text{I}} u\_{\beta, l\text{J}} \tag{4}$$

$$U\_3 = \frac{1}{6} \sum\_{m|a} \sum\_{l|\beta} \sum\_{nK\gamma} \Phi\_{mI,l}^{a\beta\gamma} u\_{a,mI} u\_{\beta,l} u\_{\gamma,nK} \tag{5}$$

where Φ*αβ mI*,*lJ* and <sup>Φ</sup>*αβγ 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 Hamiltonian*H*HA <sup>¼</sup> <sup>1</sup> 2 P *mI<sup>α</sup>mIu*\_ <sup>2</sup> *<sup>α</sup>*,*mI* þ *U*2, which includes the kinetic energy of the atom *Assessment of the Heat Capacity by Thermodynamic Approach Based on Density Functional… DOI: http://dx.doi.org/10.5772/intechopen.104083*

*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

$$\sum\_{\beta,J} D^{a\beta}\_{I\dot{\jmath}} \left( \overrightarrow{q} \right) \varepsilon\_{\beta J} \left( \overrightarrow{q}, \overrightarrow{j} \right) = \alpha^2\_{\overrightarrow{q}, \dot{\jmath}} \varepsilon\_{a, l} \left( \overrightarrow{q}, \overrightarrow{j} \right) \tag{6}$$

of the dynamical matrix

$$D\_{ll}^{a\beta} \left( \overrightarrow{q} \right) = \frac{1}{\left( m\_{I} m\_{I} \right)^{\natural\_{\ddagger}}} \sum\_{l} \Phi\_{0l,ll}^{a\beta} e^{i\overrightarrow{q} \left( \overrightarrow{r} \, y - \overrightarrow{r}\_{\ddagger} \right)} \tag{7}$$

The solutions *ω<sup>q</sup>* !,*<sup>j</sup>* and *eα*,*<sup>I</sup> q* !, *j* � � 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 *Dαβ IJ 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α*,*<sup>I</sup> q* !, *j* � � 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* !,*<sup>j</sup>* ¼ exp ℏ*ω<sup>q</sup>* !,*j =kBT* � � � <sup>1</sup> � ��<sup>1</sup> . 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 � N*a*, where N*<sup>a</sup>* is the total number of atoms in the unit cell of the crystal.

The relation between *ω<sup>q</sup>* !,*<sup>j</sup>* and *q* ! for each mode *<sup>j</sup>*, namely *<sup>ω</sup><sup>j</sup>* <sup>¼</sup> *<sup>ω</sup><sup>j</sup> <sup>q</sup>* !� �, is called phonon dispersion. The number of the phonon modes with the frequency between *ω* and *<sup>ω</sup>* <sup>þ</sup> *Δω* gives the phonon density of states *<sup>g</sup>*ð Þ¼ *<sup>ω</sup>* <sup>1</sup>*=N*<sup>P</sup> *q* !,*j δ ω* � *ω<sup>q</sup>* !,*j* � �, where *<sup>N</sup>* is the number of unit cells in the crystal. The normalization factor 1*=*ð Þ 3 � N*<sup>a</sup>* is introduced in order to reduce to 1 the integral of *<sup>g</sup>*ð Þ *<sup>ω</sup>* over frequency <sup>Ð</sup> *g*ð Þ *ω dω* ¼ 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 *Z N*ð Þ¼ ,*V*, *<sup>T</sup>* <sup>P</sup> *i* e �*Ei*ð Þ *<sup>N</sup>*,*<sup>V</sup> kBT* , 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

$$\mathbf{U} = k\_B T^2 \left(\frac{\partial \ln Z}{\partial T}\right)\_{N,V} - \text{average potential energy} \tag{8}$$

$$\mathcal{S} = k\_B \ln Z + U/T-\text{entropy} \tag{9}$$

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

$$F = U - T \cdot \mathbf{S} - \text{Helmholtz free energy} \tag{11}$$

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* þ *P* � *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 P*ð Þ¼ , *T F* þ *p* � *V* ¼ *H* � *T* � *S*.

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


*Assessment of the Heat Capacity by Thermodynamic Approach Based on Density Functional… DOI: http://dx.doi.org/10.5772/intechopen.104083*

• the isobaric heat capacity *CP* <sup>¼</sup> *<sup>∂</sup><sup>H</sup> ∂T* � � *<sup>P</sup>* <sup>¼</sup> *CV* � *<sup>T</sup> <sup>∂</sup><sup>V</sup> ∂T* � �<sup>2</sup> *P ∂P ∂V* � � *<sup>T</sup>*, and the volumetric thermal expansion as *<sup>α</sup><sup>V</sup>* <sup>¼</sup> <sup>1</sup> *V ∂V ∂T* � � *<sup>P</sup>*. 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:

$$F(V, T) = U\_0(V) + F\_{el}(V, T) + F\_{phonon}(V, T) \tag{12}$$

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 *<sup>f</sup>*ð Þ¼ <sup>∈</sup>, *<sup>T</sup>* exp <sup>∈</sup> �*EF kBT* � � <sup>þ</sup> <sup>1</sup> � ��<sup>1</sup> , where the *EF* is the Fermi level, the last occupied electronic level, with the value given by the normalization condition Ð *f*ð Þ ∈ , *T d*∈ ¼ 1. The energy due to the electronic excitation is *Uel*ð Þ¼ *<sup>V</sup>*, *<sup>T</sup>* <sup>Ð</sup> *n V*ð Þ , <sup>∈</sup> *<sup>f</sup>* <sup>∈</sup>*d*<sup>∈</sup> � <sup>Ð</sup> *n V*ð Þ , ∈ ∈*d*∈, where *n V*ð Þ , ∈ is the electronic density of states computed for the volume *V*. The corresponding electronic entropy is *Sel*ð Þ¼� *V gkB* Ð ½ � *f V*ð Þ , *T* ln *f V*ð , *T*Þ þ ð Þ 1 � *f V*ð Þ , *T* ln 1ð Þ � *f V*ð Þ , *T* .

where *g* is equal to 1 for collinear spin polarized or 2 for non-spin-polarized systems. The electronic Helmholtz free energy is *Fel*ð Þ¼ *V*, *T Uel*ð Þ� *V*, *T TSel*ð Þ *V*, *T* . The electronic excitation term is calculated based on the Mermin theorem [54]. *Fel*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 selfconsistent 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

$$Z\_{vib} = \prod\_{\vec{q}, \vec{j}} \left( \sum\_{n=0}^{\infty} e^{-\left(n + \frac{\hbar}{2}\right) \frac{\vec{q} \cdot \vec{j}}{k\_B T}} \right) \text{ and the vibrational Helmholtz energy in the harmonic oscillator } \vec{q}$$

*approximation* is depending on temperature only by the phonon contribution calculated for the equilibrium volume *Veq*

$$F\_{\rm HA}(V\_{eq}, T) = U\_0(V) + \frac{1}{2} \sum\_{\vec{q}, j} \hbar o\_{\vec{q}, j} + k\_B T \sum\_{\vec{q}, j} \log \left( 1 - e^{-\frac{\hbar \nu\_{\vec{q}, j}}{k\_B T}} \right) \tag{13}$$

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

$$\mathbf{C}\_{V} = -T \left( \frac{\partial^2 F}{\partial T^2} \right)\_V = k\_B N\_d \sum g(\alpha) \left( \frac{\hbar \alpha}{k\_B T} \right)^2 \frac{\exp \left( \hbar \alpha\_{\overrightarrow{q}, \overrightarrow{j}} / k\_B T \right)}{\exp \left( \hbar \alpha\_{\overrightarrow{q}, \overrightarrow{j}} / k\_B T \right) - 1} \tag{14}$$

#### **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 *ω<sup>q</sup>* !,*j* ð Þ *V* and the Helmholtz free energy becomes

$$F\_{QHA}(V,T) = U\_0(V) + \frac{1}{2} \sum\_{\vec{q},j} \hbar o\_{\vec{q},j}(V) + k\_B T \sum\_{\vec{q},j} \log \left[ 1 - \exp\left( -\frac{\hbar o\_{\vec{q},j}(V)}{k\_B T} \right) \right] \tag{15}$$

or by the grouping of the terms [61].

$$F\_{QHA}(V,T) = U\_{cold}(V) + F\_{th}(V) \tag{16}$$

where *Ucold*ð Þ¼ *<sup>V</sup> <sup>U</sup>*0ð Þþ *<sup>V</sup>* <sup>1</sup> 2 P *q* !,*j* ℏ*ω<sup>q</sup>* !,*j* ð Þ *V* is the cold potential energy (T = 0 K), and *Fth*ð Þ¼ *V kBT* P *q* !,*<sup>j</sup>* log 1 � exp � <sup>ℏ</sup>*<sup>ω</sup> <sup>q</sup>* !,*j* ð Þ *V kBT* � � � � is the thermal component of the Helmholtz free energy given by the phonons.

*Assessment of the Heat Capacity by Thermodynamic Approach Based on Density Functional… DOI: http://dx.doi.org/10.5772/intechopen.104083*

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* ffi *Ucold* þ *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 *g <sup>q</sup>* !ð Þ¼ *<sup>ω</sup>*, *<sup>T</sup>* <sup>Ð</sup> *e<sup>i</sup>ω<sup>t</sup>* P*Na <sup>I</sup>*¼<sup>1</sup>*ei q*!*<sup>R</sup>* ! *<sup>I</sup> <sup>v</sup> <sup>I</sup>*ð Þ*t v* ! h i *<sup>I</sup>*ð Þ <sup>0</sup> *<sup>T</sup> v* ! *<sup>I</sup>*ð Þ 0 *v* ! h i *<sup>I</sup>*ð Þ <sup>0</sup> *<sup>T</sup> dt*, where *v* ! *<sup>I</sup>*ð Þ*t* is

velocity of each atom *I* at time *t*, and *R* ! *<sup>I</sup>* 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 temperaturedependent 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(V,T) = -\left(\frac{\partial F}{\partial V}\right)\_T = P\_{\text{stat}}(V, \mathbf{0} \,\mathbf{K}) + P\_{\text{phonon}}(V, T) + P\_{\text{act}}(V, T) \tag{17}$$

where *P*0ð Þ *V*,0K is the static pressure computed by Eq. (2), *Pphon*ð Þ *V*, *T* is phonon contribution to the pressure, and *Pae*ð Þ *V*, *T* is the anharmonic and electronic thermal pressure. The thermal pressure is calculated as the derivative of the thermal free energy

$$P\_{phonon} = -\left(\frac{\partial F\_{phonon}}{\partial V}\right)\_T = -\sum\_{\vec{q}j} \left(\frac{1}{2} + \frac{\hbar o\_{\vec{q}j}}{e^{\frac{\hbar o\_{\vec{q}j}}{k\_B T}} - 1}\right) \hbar \frac{d o\_{\vec{q}j}}{dV} = -\sum\_{\vec{q}j} \chi\_{\vec{q}j} \left(\frac{1}{2} + \frac{\hbar o\_{\vec{q}j}}{e^{\frac{\hbar j}{k\_B T}} - 1}\right) \tag{18}$$

where *γ <sup>q</sup>* !,*<sup>j</sup>* ¼ � *<sup>V</sup> ω q* !,*j <sup>∂</sup>* ln*<sup>ω</sup> <sup>q</sup>* !,*j <sup>∂</sup><sup>V</sup>* ¼ � *<sup>V</sup>* 2*ω*<sup>2</sup> *q* !,*j eq* !,*j ∂D q*! ð Þ *∂V* � � � � � � � � *eq* !,*j* � � is the so-called mode Grüneisen parameter, which characterizes the volume dependences of the frequency of the mode

*q* !, *j* � �. Similarly to Eq. (2), the pressure *P V*ð Þ¼� *dF dV* � � *<sup>T</sup>* 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ð Þþ *V PV*.

#### **6.2 Debye method**

The Debye model [68] considers a simple form for the DOS of the vibration modes *<sup>g</sup>*ð Þ¼ *<sup>ω</sup> <sup>C</sup>* � *<sup>ω</sup>*<sup>2</sup> � <sup>Θ</sup>ð Þ *<sup>ω</sup>* � *<sup>ω</sup><sup>D</sup>* , where <sup>Θ</sup> is the Heaviside step function, and *<sup>C</sup>* <sup>¼</sup> <sup>9</sup>*Na=ω*<sup>3</sup> *<sup>D</sup>* is a constant determined from the condition *C*Ð *g*ð Þ *ω dω* ¼ 3*Na*. 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

$$F\_{phonon}^{Delay\epsilon}(T) = N\_d k\_B T \left(\frac{9T\_D}{8T} + 3\ln\left(1 - e^{-\frac{T\_D}{T}}\right) - D\left(\frac{T\_D}{T}\right)\right) \tag{19}$$

where the function *D TD T* � � <sup>¼</sup> <sup>3</sup> ð Þ *TD=<sup>T</sup>* <sup>3</sup> Ð *TD=T* 0 *x*3 *ex*�<sup>1</sup> *dx* is the Debye integral and the parameter *TD* <sup>¼</sup> <sup>ℏ</sup> *kB* <sup>6</sup>*π*<sup>2</sup>*V*<sup>1</sup>*=*<sup>2</sup>*Na* � �<sup>1</sup>*=*3 *f*ð Þ *σ* ffiffiffi *Bs M* q is the Debye Temperature, which can be interpreted as the temperature at which each mode below the highest-frequency mode *ω<sup>D</sup>* 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},\overrightarrow{q}} \chi\_{\overrightarrow{q},\overrightarrow{q}} \mathbf{C}\_{V} \left(\overrightarrow{q}, \overrightarrow{j}\right)}{\sum\_{\overrightarrow{q},\overrightarrow{q}} \mathbf{C}\_{V} \left(\overrightarrow{q}, \overrightarrow{j}\right)} \tag{20}$$

is the total Grüneisen parameter, where

$$C\_V\left(\overrightarrow{q},j\right) = k\_B \left(\frac{\hbar o\_{\overrightarrow{q},j}}{k\_B T}\right)^2 \left[\exp\left(\frac{\hbar o\_{\overrightarrow{q},j}}{k\_B T}\right) - 1\right]^{-2} \exp\left(\frac{\hbar o\_{\overrightarrow{q},j}}{k\_B T}\right) \tag{21}$$

is the contribution of each vibration mode *q* !, *j* � � to the isochoric heat capacity

$$\mathbf{C}\_{V} = \sum\_{\vec{q}, \vec{j}} \mathbf{C}\_{V} \left( \vec{q}, \vec{j} \right) \tag{22}$$

The constant volume (isochoric) heat capacity was obtained from the quasiharmonic phonon frequencies calculated at each fixed volume. The isochoric heat *Assessment of the Heat Capacity by Thermodynamic Approach Based on Density Functional… DOI: http://dx.doi.org/10.5772/intechopen.104083*

capacity at each temperature is then obtained interpolating at the temperaturedependent volume values obtained at each temperature from the minimization of the free energy.

The constant pressure (isobaric) heat capacity is obtained as

$$\mathbf{C}\_{P} = \mathbf{C}\_{V} + T\mathbf{V}\boldsymbol{\beta}^{2}\mathbf{B}\_{T} \tag{23}$$

where *BT*ð Þ¼ *<sup>V</sup>*, *<sup>T</sup>* <sup>1</sup> *V <sup>∂</sup>*2*F V*ð Þ , *<sup>T</sup> ∂V*<sup>2</sup> *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�<sup>1</sup> K�<sup>1</sup> , 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.

The heat capacity in the Debye model is *CDebye <sup>V</sup>* <sup>¼</sup> <sup>3</sup>*NakB* <sup>4</sup>*<sup>D</sup> TD T* � <sup>3</sup>*TD=<sup>T</sup> <sup>e</sup>TD=T*�<sup>1</sup> . 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 3*NakB* at high temperatures.

#### **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 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=<sup>P</sup>* <sup>¼</sup> *<sup>E</sup>*<sup>2</sup> h i�h i *<sup>E</sup>* <sup>2</sup> *kBT*<sup>2</sup> , 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.

#### **Acknowledgements**

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.

*Assessment of the Heat Capacity by Thermodynamic Approach Based on Density Functional… DOI: http://dx.doi.org/10.5772/intechopen.104083*

#### **Author details**

Viorel Chihaia<sup>1</sup> \*, Valentin Alexiev<sup>2</sup> and Hasan S. AlMatrouk<sup>3</sup>

1 Institute of Physical Chemistry Ilie Murgulescu, Romanian Academy, Bucharest, Romania

2 Institute of Catalysis, Bulgarian Academy of Sciences, Sofia, Bulgaria

3 Kuwait Institute for Scientific Research, Safat, Kuwait City, Kuwait

\*Address all correspondence to: vchihaia@icf.ro

© 2022 The Author(s). Licensee IntechOpen. This chapter is distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/3.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

### **References**

[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] 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] 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] 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] 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] 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] AlMatrouk HS, Al-Ajmi F, Do NS, Chihaia V, Alexiev V. The pressuretemperature 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] 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-consistentfield ab initio calculations. The Journal of Chemical Physics. 2014;**141**:104108- 104109. DOI: 10.1063/1.4895113

[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] Folland NO. Finite-sum approximations to cubic Brillouin-zone integrals. Physical Review B. 1980;**22**: 3669-3677. DOI: 10.1103/ PhysRevB.22.3669

[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] Methfessel M, Paxton AT. Highprecision sampling for Brillouin-zone integration in metals. Physical Review B. 1989;**40**:3616-3621. DOI: 10.1103/ PhysRevB.40.3616

[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] Schwerdtfeger P. The pseudopotential approximation in electronic structure theory.

*Assessment of the Heat Capacity by Thermodynamic Approach Based on Density Functional… DOI: http://dx.doi.org/10.5772/intechopen.104083*

ChemPhysChem. 2011;**12**:3143-3155. DOI: 10.1002/cphc.201100387

[15] Kresse G, Hafner J. Normconserving 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] 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] 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] 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] 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] 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] 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] 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] 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] 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] 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] 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] 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] 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] 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] Runge E, Gross EKU. Densityfunctional theory for time-dependent systems. Physical Review Letters. 1984; **52**:997-1000. DOI: 10.1103/ PhysRevLett.52.997

[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] 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] 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] Reining L. The GW approximation: content, successes and limitations. WIREs Computational Molecular Science. 2018;**8**:e1344-e1326. DOI: 10.1002/wcms.1344

[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] Dral PO, Wu X, Spörkel L, Koslowski A, Thiel W. Semiempirical Quantum-chemical orthogonalizationcorrected methods: Benchmarks for ground-state properties. Journal of Chemical Theory and Computation. 2016;**12**:1097-1120. DOI: 10.1021/acs. jctc.5b01047

[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] 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] Dove MT. Introduction to Lattice Dynamics. Cambridge: Cambridge University Press; 1993. pp. 179-194. ISBN: 9780521392938

[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] 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

*Assessment of the Heat Capacity by Thermodynamic Approach Based on Density Functional… DOI: http://dx.doi.org/10.5772/intechopen.104083*

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] 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] 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. www.quantum-espresso.org [Accessed on 21 November, 2021]

[44] Available from: https://dalcorso.gith ub.io/thermo\_pw/ [Accessed on 21 November, 2021]

[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/PhysRevLe tt.78.1396

[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] 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] 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] 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] 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] 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] 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] 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] Mermin ND. Thermal properties of the inhomogeneous electron gas. Physics Review. 1965;**137**:A1441-A1443. DOI: 10.1103/PhysRev.137.A1441

[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] Ashcroft NW, Mermin ND. Solid State Physics. Philadelphia: Saunders College Publishing; 1976. ISBN: 0-03- 083993-9, 9780030839931

[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] 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] Born M, Huang K. Dynamical Theory of Crystal Lattices. New York: Oxford University Press; 1988. ISBN: 978-0198503699

[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] Zhang H, Shang SL, Wang Y, Saengdeejing A, Chen LQ, Liu ZK. Firstprinciples 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] 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] 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] 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] 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] 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] Wang CZ, Chan CT, Ho KM. Tightbinding 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] Debye P. Zur Theorie der spezifischen Wärmen. Annals of Physics. 1912;**39**:789-839. DOI: 10.1002/ andp.19123441404

[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] Deus P, Schneider HA, Voland U. Estimation of the Debye temperature of *Assessment of the Heat Capacity by Thermodynamic Approach Based on Density Functional… DOI: http://dx.doi.org/10.5772/intechopen.104083*

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] Wolf U, Bohmhammel K, Wolf G. Supports open access. Thermochim Acta. 1998;**310**:37-42. DOI: 10.1016/ S0040-6031(97)00382-1

[72] NIST-JANAF. Thermochemical Tables. Last Update to Data Content. College Park, Maryland, United States DOI: Available from: https://janaf.nist.g ov: American Institute of Physics; 1998 [Accessed on 21 November, 2021]

[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] 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] 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

Section 2
