Reactive force field development for magnesium chloride hydrates and its application for seasonal heat storage

• A submitted manuscript is the version of the article upon submission and before peer-review. There can be important differences between the submitted version and the official published version of record. People interested in the research are advised to contact the author for the final version of the publication, or visit the DOI to the publisher's website. • The final author version and the galley proof are versions of the publication after peer review. • The final published version features the final layout of the paper including the volume, issue and page numbers.


Introduction
Thermal energy storage is necessary for the successful implementation of solar thermal energy. 1 Solar thermal energy can be stored in three forms: sensible, latent, and thermochemical. 2 Thermal storage is a proven concept but it is volume restricted. The energy storage density is higher in thermochemical form than in sensible or latent form. Therefore thermochemical materials (TCMs) may offer a compact seasonal heat storage solution. 3 Salt hydrates are one class of TCMs. 1 These materials convert thermal energy (physical form of energy) into chemical form by decomposing into lower hydrate form or anhydrous form (charging cycle). Anhydrous and lower hydrates recombine with water vapor (present in moisture) to form hydrates again, thereby releasing energy (discharging cycle). The operating principle of TCM can be expressed with the following reaction: In addition to higher energy storage capacity, TCMs can store thermal energy with almost no heat loss, even for longer periods, making TCM a unique seasonal heat storage material. 4,5 Selection of an appropriate salt hydrate is the key aspect in designing a thermochemical heat storage system. 6 Halide based salt hydrates (CaX 2 Á6H 2 O, MgX 2 Á6H 2 O, and SrX 2 Á6H 2 O, X = Cl/Br), sulphate based salt hydrates (LiSO 4 ÁH 2 O, MgSO 4 Á7H 2 O), and nitrate based salt hydrates (LiNO 3 Á3H 2 O, Zn(NO 3 ) 2 Á3H 2 O) are the major class of promising reversible TCMs. 6 These hydrates have theoretical crystal storage capacities between 2-3 GJ m À3 . 6 Magnesium chloride hydrates (MgCl 2 ÁnH 2 O, n = 0, 1, 2, 4, and 6) have been investigated extensively as an appropriate candidates for solar thermal storage. [7][8][9][10][11][12][13] MgCl 2 Á6H 2 O (Bischofite) has been successfully tested for 1000 cycles for application in solar based cooking. 11 Thermal dehydration of Bischofite is an important process in production of pure Mg, Mg alloys, MgO, and in heat storage applications. [14][15][16][17] The dehydration kinetics for MgCl 2 hydrates is faster than for other materials like MgSO 4 hydrates, 8,18 but a problem with this material is the presence of the hydrolysis along with the dehydration. Hydrolysis is an irreversible side reaction in the dehydration process, which produces HCl gas and affects the durability of the heat storage system. It was experimentally observed that the hydrolysis reaction could be avoided by use of external HCl pressure 14 or by adding other halides. 19 However, there is no fundamental understanding of the conditions for the hydrolysis.
Next to the dehydration, the HCl production is also a key aspect in designing the MgCl 2 based heat storage system. Experimentally, thermal dehydration and hydrolysis of MgCl 2 hydrates have been studied by thermogravimetric analysis (TGA), 17 derivative thermogravimetry (DTG), 20 TGA/Differential scanning calorimetry (DSC) 10 and thermochemical analysis. 14 Galwey et al. 21 has experimentally observed the hydrolysis reaction in MgCl 2 Á2H 2 O in the temperature range of 623-703 K. 21 Kirsh et al. 20 has not observed hydrolysis below 473 K for MgCl 2 hydrates, but found 43.8% of chlorine released as HCl at 673 K. However, Huang et al. 17 has reported the onset temperature of HCl formation at 440 K in MgCl 2 Á2H 2 O. Kipouros et al. 14 has reported hydrolysis at 417 K for MgCl 2 Á2H 2 O and at 413 K for MgCl 2 ÁH 2 O for an equilibrium HCl pressure of 5 Â 10 À3 atm. The onset temperature of HCl formation is thus ambiguous from experiments. There is no theoretical study to determine the hydrolysis kinetics and the onset temperature of HCl formation.
Atomic and molecular level simulations turned to be useful tools to gain insight in the hydrolysis coupled dehydration reactions. Weck et al. 12 has computed thermal properties from thermodynamical calculations of MgCl 2 ÁnH 2 O using Density Function Theory (DFT) under harmonic assumption. Wang et al. 16 has considered MgCl 2 ÁH 2 O to be the potential candidate for hydrolysis and proposed a reaction path for thermolysis of MgCl 2 ÁH 2 O from DFT calculations. Smeets et al. 8 has obtained the safety limit of hydrolysis reactions for all hydrates of MgCl 2 from DFT calculations using equilibrium thermodynamic principles under the ideal poly-atomic gas assumption. The atomic scale calculations can simulate systems upto a system size of 100-1000 atoms and a time scale upto 200 picosecond over a week of high performance computing. 22 To the best of our knowledge, no kinetic study was done on molecular level to investigate the hydrolysis and dehydration reactions of MgCl 2 ÁH 2 O and MgCl 2 Á2H 2 O from first principle atomic level calculations, integrating the quantum level results to molecular level modeling under reactive force field (ReaxFF) formalism.
For an efficient TCM storage cycle, heat and mass transfer in the storage volume are essential. 2 Experimentally it has been observed that water transport through the solid MgCl 2 hydrate affects the dehydration rate and may become the rate limiting step under specific reaction conditions. 23 To understand the complex reaction (hydration and dehydration) coupled with water transport, it is desirable to investigate the diffusive transport of water molecules, which are formed by the dehydration reaction. These water molecules have to diffuse to the surface of the solid salt hydrate. The diffusion path may vary from several atomic layers up to a few micrometer. 24 The water transport through the solid salt hydrate, the temperature, the diffusion path length, the external H 2 O pressure and the crystal defects influence considerably the dehydration rate. 25 To study the effect of these aspects on the complex dehydration reaction on molecular level, a force field is required that is able to capture both the chemical reactivity and the mass transport of H 2 O in MgCl 2 ÁH 2 O and MgCl 2 Á2H 2 O crystals.
Chemical reactivity has been incorporated using bond order dependent ReaxFF. 26,27 According to previous studies, 28 In the present study, we parameterize the ReaxFF to model the dehydration kinetics of MgCl 2 ÁH 2 O and MgCl 2 Á2H 2 O, as these hydrates are prone to hydrolysis reactions. Proton transfer is an important step in the hydrolysis reaction. The energy barrier for proton transfer in MgCl 2 Á2H 2 O is obtained from the ReaxFF simulations and validated against results from DFT. The extent of the reaction depends on the temperature of the solid and the vapor transport inside the solid. The new ReaxFF is used to study the effect of temperature on the dehydration in the temperature range of 300-500 K. This temperature range is the typical operating range of TCM based heat storage systems. The onset temperature of HCl formation is investigated in the same temperature range. The dependence of the diffusion coefficient (DC) on temperature is analyzed using the Green-Kubo (GK) method.

Quantum mechanics calculations
To accurately predict the dehydration and hydrolysis reaction in the condensed phase of MgCl 2 hydrates, quantum mechanics (QM) calculations are performed in both periodic and non-periodic systems in order to generate a training set necessary to parameterize the new force field. DFT under the generalized gradient approximation (GGA) is used for the QM calculations. All the non-periodic MgCl 2 hydrate clusters are relaxed in the framework of DFT implemented in Amsterdam Density Functional (ADF) program. 34 The Perdew-Wang exchange and correlation functional (PW91) is used with a double polarized triple-z basis set. A similar GGA-DFT framework has been used to study the CaCl 2 , MgCl 2 , and MgSO 4 hydrates. 8 To simulate a reaction in solid phase of hydrates, the crystal structure of MgCl 2 , MgCl 2 ÁH 2 O, and MgCl 2 Á2H 2 O are compressed or stretched uniformly upto 20-30% of their relaxed volume to obtain the EOS. The periodic systems are relaxed in the framework of GGA-DFT using the VASP package. 36 The projector augmented wave (PAW) method 37 along with the Perdew Burke Ernzerhof (PBE) exchange-correlation functional 38 is used under the GGA approximation. The Monkhorst-Pack k-point mesh is used to integrate the Brillouin zone of MgCl 2 ( The unit cell of crystals are fully relaxed with break condition for ionic relaxation of 10 À6 eV and plane wave cut off of 500 eV using tetrahedron method with Blöchl corrections.

Reactive force field and parameterization
ReaxFF is a general bond-order-dependent many body empirical potential with a polarizable charge model that dynamically predicts a bond formation and dissociation during a reaction without connecting step function potentials. 39 ReaxFF formalism represents a continuous landscape of energy which is dependent on the bond order and the non-bonded distance. The bond order is obtained from the inter-atomic distance and energy contribution updated with time. In the ReaxFF formalism, the total energy of the system can be partitioned into the following terms: E total is the overall interaction energy of the system, E Coulomb is the electrostatic contribution, and E bond is the bond energy. E val and E tor are the three-body valence angle terms and the fourbody torsion terms respectively. E lp and E over/under are the energy contributions from lone-pair electrons, and the penalty energy coming from over coordination and under coordination, E conj represents the conjugation energy term, E H-bond represents the weak hydrogen bond and E others is introduced for other correcting types of species. E Coulomb and E vDW are the nonbonded interactions between all pairs irrespective of the connectivity in the system. The Coulomb charge is calculated dynamically from the electro-negativity equalization method. 39 Excessive electrostatic repulsions at very short interatomic distance are avoided by shielding parameters.
The ReaxFF parameters required to define the potential energy surface are optimized against the training set obtained from QM data. The initial O/H parameters are taken from the iron-oxyhydroxide force field. 40 The O/H parameters are kept fixed while bond parameters of Mg-O, Mg-Cl and angle parameters associated with these Mg and Cl are optimized. The optimization is done in such a way that optimized parameters are able to reproduce the energy and structural data present in the training set. The data in the training set used for comparison include bond dissociation energy, angle bending curve, reaction enthalpy and EOS for MgCl 2 , MgCl 2 ÁH 2 O, and MgCl 2 Á 2H 2 O. In an optimization process, ReaxFF parameters are varied such that the cumulative error (X) between the DFT data in the training set and their corresponding values obtained from ReaxFF is minimized. The cumulative error is defined as In this equation, X i,DFT is the DFT value in the training set for data point i, and X i,ReaxFF is the corresponding value obtained by ReaxFF. s i is a weight assigned to the ith data point to set its relative importance compared to other data points. In a singleparameter search algorithm, 32 a parabolic extrapolation method is used for parameterization. Further, we have used an MMC algorithm 33 to attain global minima of the cumulative error as ReaxFF parameters are correlated. In the MMC algorithm multiple parameters are changed with a random movement at every iteration and the proposed configuration is accepted with a probability given by E new is the error from the new set of parameter proposition and E old from the old set. T is an artificial temperature which gradually decreases after each iteration (simulated annealing) and k B is the Boltzmann constant. The MMC algorithm has shown a good agreement between QM and ReaxFF data for MgSO 4 hydrates systems. 33

Diffusion
The DC of H 2 O molecules through the MgCl 2 hydrate crystal is obtained from the GK formulation described as D is the diffusion coefficient of water molecules. The expression hv(t)Áv(0)i is known as the velocity autocorrelation function (VAF). 41 The DC of H 2 O molecules through the MgCl 2 hydrate crystals is obtained in a spherical control volume of radius R in such a way that the crystal slab resides completely in the sphere and the center of mass of this slab is chosen as the center of the sphere as shown in Fig. 9 3 Results and discussion

Force field parameterization results
The initial force field parameters are optimized against hydrated Cl À ion and MgCl 2 condensed phase data.  Fig. 3.
To predict the kinetics of solid MgCl 2 ÁH 2 O and MgCl 2 Á2H 2 O, the EOS is essential. The optimized force field, which is able to describe the gas phase reactions accurately, is further optimized with a condensed phase EOS of MgCl 2 hydrates obtained from DFT in VASP. 36 The parameterized force field is also able to predict the experimentally observed crystal structures of MgCl 2 , 42 Fig. 4-6.
The bulk modulus is obtained from fitting the Birch-Murnaghan equation of state with the EOS obtained from DFT and ReaxFF. The bulk modulus (B 0 ) and its first order derivative with respect to pressure (B 0 0 ) (rows 7 and 8 of

Proton transfer reaction
Hydrolysis is an irreversible reaction which is usually observed in MgCl 2 ÁH 2 O and MgCl 2 Á2H 2 O. Wang and Chen 16 have studied the thermolysis mechanism of MgCl 2 ÁH 2 O using a semi-empirical PM3 method. They observed that proton transfer is an important step in the hydrolysis reaction. To understand the reaction path for the proton transfer in MgCl 2 Á2H 2 O, the relative energy of a     Fig. 7(b). A barrier of 29.61 kcal mol À1 is obtained from DFT for proton transfer in MgCl 2 Á2H 2 O surrounded by one H 2 O molecule. This explains the lower hydrolysis in the higher hydrates (tetra and hexa) of MgCl 2 , as the neighboring H 2 O molecule increases the barrier for proton transfer and inhibits the hydrolysis. A barrier of 26.54 kcal mol À1 for proton transfer is obtained from ReaxFF, which is in agreement with the DFT result (29.61 kcal mol À1 ). The position of the transition state predicted from ReaxFF differs by 0.5 Å because surrounding H 2 O molecule may form many structures of similar energy (see ESI †).

ReaxFF-MD simulations of MgCl 2 hydrates
In the present study, we have investigated the kinetics of dehydration and hydrolysis reactions of solid MgCl 2 ÁH 2 O and MgCl 2 Á2H 2 O.      43 A similar peak location was reported from the ReaxFF parameters of Lindqvist Polyoxoanion in Bulk Water. 46 The first peak location for the Mg-O pair is observed at 2.0 Å, which represents the Mg-O atomic bond as shown in Fig. 10. The corresponding distance in the experimental unit cell of MgCl 2 Á2H 2 O is 2.0 Å. 44 The second peak location for Mg-O pair is observed at 4.4 Å, which represents the non-bonded Mg attractions with the H 2 O molecule of the adjacent MgCl 2 hydrate as shown in Fig. 10. The corresponding distance in the experimental unit cell of MgCl 2 Á2H 2 O is 4.0 Å 44 in this case.
The .08 in the temperature range of 300 to 500 K as shown in Fig. 11. A molecular dehydration rate could be estimated from  the gradient of a linear fit. The dehydration rate keeps on increasing with temperature in the range of 300 to 500 K as shown in Fig. 11. The dehydration rate at 500 K is 3.1 times faster when compared to the rate at 300 K.
The average number of H 2 O molecules leaving the solid slab of MgCl 2 Á2H 2 O (44.3 Â 1000 Â 21.9 Å) into the vacuum region (Y direction) is 51.17 AE 7.63 after 750 picosecond. The system is further equilibrated upto 1000 picosecond and H 2 O average increased by 4.61. The average density distribution of H 2 O molecules are obtained by diving the simulation box into the bins of 10 Å along the Y axis as shown in Fig. 9b. Blue lines represent the initial boundaries. The water molecules which escape from the solid slab accumulate on the surface at 300 K. There is a concentration gradient of H 2 O molecules present on the surface along the Y axis. Thereby, confirms the influence of water diffusion on the dehydration process. The average number of H 2 O molecules escaping from MgCl 2 Á2H 2 O increases from 54.93 AE 6.76 to 71.66 AE 7.27 in the temperature range of 300 to 450 K as shown in Fig. 12. Further, increment in temperature leads to reduction in the number of H 2 O molecules escaping from MgCl 2 Á2H 2 O to 50.73 AE 5.85 at 500 K as shown in Fig. 12. This reduction in number of water molecules escaping from solid slab of MgCl 2 Á2H 2 O could be explained from the fact that existing water molecules present in the vacuum region are pushing them back into the crystal. The effect of temperature on the dehydration of MgCl 2 ÁH 2 O and MgCl 2 Á2H 2 O is similar to experiments. 47 3.3.3 HCl formation. The HCl formation is investigated at various temperatures (300, 350, 400, 450, and 500 K). The first HCl molecule is observed from the 2D periodic slab of MgCl 2 ÁH 2 O at 350 K after 136.25 picosecond. To narrow down the onset temperature of HCl formation, MD simulations are carried out at the interval of 10 K between 300-350 K. The onset of HCl formation is observed at 340 K after 318.75 picosecond in MgCl 2 ÁH 2 O as shown in Fig. 13b. After dehydration/hydrolysis, the crystal slab elongates along the vacuum direction as some of MgCl 2 and H 2 O came out on the surface. The onset of HCl formation is experimentally 23 observed in the temperature range of 343-353 K, which is in good agreement with the present ReaxFF study. The hydrolysis  kinetics of the MgCl 2 ÁH 2 O is shown in Fig. 13c at different temperatures (350-500 K) upto 1000 picosecond. The number of HCl molecules keeps on increasing with temperature. The hydrolysis rate estimated from the linear fits keeps on increasing with temperature as shown in Fig. 13c 50 The DC follows the Arrhenius equation. The activation energy and pre-exponential factor are 5.7 Â 10 À4 m 2 s À1 and 42.6 kcal mol À1 . Thus, the strength of Mg-O    49 and experiments. 50 The DC increases with temperature and follows the Arrhenius law for MgCl 2 ÁH 2 O and MgCl 2 Á2H 2 O. The diffusivity of H 2 O through MgCl 2 ÁH 2 O at 500 K increases 64 times compared with diffusivity at 300 K. Similarly, the diffusivity of H 2 O through MgCl 2 Á2H 2 O at 500 K increased 57 times compared to the 300 K. The diffusivity trend suggests that by increasing the temperature upto 500 K, the H 2 O transport (mass diffusion) can be improved. These results demonstrate the ability of ReaxFF to explore the molecular reaction rate for dehydration/hydrolysis of MgCl 2 hydrates along with the water transport in the operational range of seasonal heat storage systems.