Reactive Grand-Canonical Monte Carlo Simulations for Modeling Hydration of MgCl2.

Thermochemical heat-storage applications, based on the reversible endo-/exothermic hydration reaction of salts, are intensively investigated to search for compact heat-storage devices. To achieve a truly valuable storage system, progressively complex salts are investigated. For these salts, the equilibrium temperature and pressure conditions are not always easy to predict. However, these conditions are crucial for the design of thermochemical heat-storage systems. A biased grand-canonical Monte Carlo (GCMC) tool is developed, enabling the study of equilibrium conditions at the molecular level. The GCMC algorithm is combined with reactive force field molecular dynamics (ReaxFF), which allows bond formation within the simulation. The Weeks-Chandler-Andersen (WCA) potential is used to scan multiple trial positions for the GCMC algorithm at a small cost. The most promising trial positions can be selected for recomputation with the more expensive ReaxFF. The developed WCA-ReaxFF-GCMC tool was used to study the hydration of MgCl2·nH2O. The simulation results show a good agreement with experimental and thermodynamic equilibriums for multiple hydration levels. The hydration shows that water, present at the surface of crystalline salt, deforms the surface layers and promotes further hydration of these deformed layers. Additionally, the WCA-ReaxFF-GCMC algorithm can be used to study other, non-TCM-related, reactive sorption processes.

Introduction

Energy-storage systems are a vital link in the sustainable energy infrastructure. Thermochemical energy storage can fulfill this essential link regarding thermal applications and is therefore widely studied.[1−5] Thermochemical energy storage is compelling because it can realize a relatively high thermal energy-storage density and no thermal losses occur during storage.[2] The concept relies on a reversible endo-/exothermic chemical reaction between a sorbate (A) and a sorbent (B)

When there is a surplus of (thermal) energy, heat is used to separate material AB—via an endothermic reaction—into components A and B. These components are stored separately, and no energy dissipates during storage. In times of a lack in (thermal) energy, components A and B are combined and react—via an exothermic reaction—back to AB, releasing heat. The operation of the thermochemical energy storage is dictated by the given equilibrium conditions of the reaction as described by thermodynamics.[6,7]

Ideal thermochemical materials (TCMs) for such storage require high reaction energies between the sorbate and sorbent, stability over many storage cycles, and a reaction equilibrium around desired operating temperatures and partial pressures of the sorbate. Salt hydrates are promising TCMs[8−10] for applications in the built environment because of the high sorption energy of water vapor and dehydration conditions, which are reachable by built environment climate systems. Their endo-/exothermic reversible chemical reaction of (de)hydration is described by

Most common salts in their pure form suffer from drawbacks regarding their application as TCM.[8−10] For example, undesired irreversible side reactions which degrade the storage capacity;[11,12] slow kinetics, and low thermal conductivities,[13] which decrease thermal in- and output power; a metastable zone around the equilibrium conditions;[6,7] or the occurrence of undesired melting or deliquescence of the salts, which forms blocking agglomerates in the storage system. Deliquescence is the phenomenon of a water-soluble substance that absorbs so much water vapor from the atmosphere that the substance will dissolve in its own absorbed water. These drawbacks lower the stability and cyclability of the TCM and, thereby, the storage system.

Recently, much research focuses on enhancing pure salts to overcome the previously mentioned drawbacks, for example, a mixture of multiple salts,[14] double salts,[15,16] encapsulation,[17] composites,[18] impregnation of salts in porous materials,[19] or doping of pure salts.[16,20,21] For most common salts, the equilibrium conditions can be estimated from thermodynamic tables.[22,23] However, due to the increased complexity of interacting elements that influence the equilibrium conditions or missing thermodynamic tables for new enhanced salts, the reaction equilibrium description becomes increasingly complex compared to pure salts. As a result, the equilibrium cannot always be described by simple thermodynamic rules without assumptions or educated guesses for unknown model parameters.

Molecular modeling could act as a solution; it can predict many material properties for these new complex TCMs. In this regard, quantum mechanical (QM) methods, such as density functional theory (DFT), were used before to investigate chemical bonding between a new combination of sorbate and sorbent regarding heat-storage systems[16] or to predict equilibrium conditions for gas-phase systems.[24] However, due to the high computational cost of QM simulations, it is limited to a relatively small number of atoms over a short time period. On a larger scale, force field-based molecular dynamics (MD) is a powerful method, which is computationally much cheaper and therefore applicable to much larger systems over a much longer time. However, when standard classical force fields are used, the ability to model bond breaking and formation like in DFT is lost. Thereby, the study of TCMs including reactions is not possible. This gap between quantum mechanical methods and MD, in terms of modeling bond formation and the required computational cost, is bridged by reactive force fields (ReaxFF).[25−27] ReaxFF is able to model bond breaking and formation but only at a slightly higher computational cost than nonreactive MD with classical force fields. As a result, ReaxFF has been used before to study dynamic properties such as hydrolysis, diffusion, and dehydration of TCMs.[12,13,28] Despite the lower computational cost, ReaxFF is still not used to study hydration, deliquescence, or equilibrium conditions of TCMs. These phenomena are related to rare events and/or high energy barriers, and the MD time scale is too short to overcome these barriers.

Molecular Monte Carlo (MC) methods are closely related to MD. However, where MD uses time integration methods to sample the phase space, MC is a stochastic approach to relate statistical properties to a mechanical property.[29] Because MC does not depend on the dynamical method of time integration, these high energy barriers and rare events can be circumvented by smart MC algorithms.

Herein, we introduce an advanced grand-canonical MC (GCMC) model which is able to predict the chemical equilibrium at a given temperature and pressure between a sorbate and sorbent.[30] Moreover, it can reveal changes in the reaction equilibrium of salts upon chemical or physical enhancement. With the GCMC algorithm, molecules within a simulated system are exchanged with an infinite large reservoir at a constant chemical potential. These exchanges are accepted according to acceptance rules that enclose the phase space of the system. Consequently, for a large number of molecule exchanges, the system’s chemical potential will equilibrate with the reservoir’s chemical potential. In this way, a relation is established between the reservoir with an imposed pressure and temperature and the number of molecules in the system. This makes GCMC practically preferable to study sorption over time-dependent MD simulations that are hindered by a limited computationally available time regarding processes with high energy barriers like diffusion in confined regions.

MC in combination with reactive systems has been applied before, for example, reaction ensemble MC.[31] However, for such a system, the reaction product should be known. This is not the case for ReaxFF MD, where the reaction will follow from the dynamics, given the required chemical environment. The combination between GCMC and ReaxFF has been introduced by Senftle et al.[32−34] It has been used in catalytic studies of oxidation,[32,34,35] hydrogenation,[32] and carbonation.[32] Islam et al.[36] used the same method for battery applications and studied lithium insertion in α-sulfur. Jung et al.[37] developed a ReaxFF grand-canonical MD (GCMD) combination, which performed GCMC trial moves at predefined intervals in ReaxFF-MD. These authors showed its applicability to gas-phase water formation from oxygen and hydrogen molecules on platinum catalysts. All the previous studies either focus on the gas-phase abovementioned surfaces[34,37] or on bulk system interactions with monatomic sorbents.[32,33,35,36] When considering the gas phase, insertions are still practically feasible due to a large amount of available voids. However, in dense bulk materials, GCMC insertion of molecules becomes cumbersome due to the high possibility of overlapping atoms when random insertions are performed. Insertions that result in overlapping atoms will be rejected due to the corresponding high energy, and these low acceptance probabilities make the system hard to reach equilibrium. To avoid overlaps between atoms, advanced GCMC algorithms are successfully developed[30,38−41] for dense systems, where a bias is used to prevent overlaps. In the existing ReaxFF-GCMC combination, molecular insertions have been increased using a forward bias that performs energy minimization after each insertion. Accordingly, this bias has to be counterbalanced to avoid unrealistic overloading. Senftle et al.[32] compensated the bias by reducing the accessible volume.

To study the hydration of MgCl2·nH2O, non-monatomic gas—H2O—molecules need to be inserted. To counterbalance energy minimization for an inserted H2O molecule, via an assumed reduced volume, would be increasingly ambiguous for a molecule compared to atomic insertion. In this work, a novel ReaxFF–GCMC combination is developed with an alternative biasing scheme to increase its insertion efficiency. To avoid excessive energy calculations with ReaxFF for unrealistic overlapping insertions, first, k trial insertions are performed with a computationally cheap short-range Weeks–Chandler–Andersen (WCA)[42] interaction potential. From these cheaply calculated insertion trials, realistic insertions without hardcore overlaps can be selected for recalculation with ReaxFF. Thus, this bias favors insertions that are more likely to be accepted, and as a result, equilibrium can be reached for dense systems. Because the number of trial positions (k) and their corresponding energies are exactly known, this bias is exactly known and can be counterbalanced in the acceptance rules of the GCMC algorithm.

The article is organized as follows. In Section , the GCMC algorithms are validated and accordingly used to study and discuss the hydration of MgCl2·nH2O. In Section , conclusions are drawn concerning the algorithm and its future potential use to study equilibrium conditions. In Section , the basic GCMC algorithm is explained. Furthermore, the biased WCA–ReaxFF–GCMC algorithms are explained, which improve the insertion of molecules, together with the WCA and ReaxFF potentials.

Results and Discussion

In Sections and 2.2, the biased GCMC algorithms are tested and validated. In Sections and 2.4, the WCA–ReaxFF–GCMC algorithm is used to study MgCl2·nH2O hydration. Images to visualize MgCl2·nH2O hydration are created with iRASPA.[43]

WCA–GCMC Validation

To test and validate the WCA–ReaxFF combination with the GCMC algorithm, the hydration of a MgCl2 crystal with an artificial cavity was studied for different numbers of k trial positions.

In Figure , the resulting loading of H2O is shown. These systems were simulated at a given temperature of 300 K and a vapor pressure of 12 mbar. The simulation with k = 1 is considered as a reference system, where it represents the conventional (unbiased) GCMC algorithm with k = 1. The computational times required for the other simulations were normalized to this k = 1 reference. This reference system took 5.6 h for 30,000 MC moves, initially including 288 Mg/Cl atoms, with single-core MC moves (insertion/deletion) energy calculations and 8-core ReaxFF–MD (translation) calculations, on a Haswell (Intel Xeon Processor E5-2690 v3) node. From Figure , the gain by the WCA bias is clear. The k = 1 reference is loading much slower than the other structures and does not load more than 28 H2O water molecules over the entire simulation, where the other systems (k = 5, k = 10, k = 20, and k = 40) already reach a similar loading at ca. 5–7% of the computational time. Furthermore, with multiple k trial positions, the final loading is much higher. The k = 40 test even reaches a loading of 44 water molecules, at only 121% of the time of the k = 1 simulation. It is doubtful if the reference k = 1 case would ever reach such a loading within acceptable computational time, where it was not successful for the last ∼40% of the time to insert another water molecule. In the Supporting Information, we show for a simple system, a rarefied system where hardcore overlap is hardly present for insertion, that the number of k-trial positions does not change the final equilibrium. The insets of Figure represent the systems at the indicated times. After approximately 5% of the time, the first successful insertion of a H2O molecule outside the initial cavity occurred. Before this point, hydration only happened in the cavity. This indicates that the initial bulk MgCl2 layers are too densely packed to accept a H2O molecule, and the created H2O water layer in the cavity is needed to deform the densely packed layers and create a disordered region of salt and water. This phenomenon is more extensively discussed in Section .

Figure 1

Loading of MgCl2 cavities, with k trial positions, normalized to the total required simulation time for k = 1. The insets are the structures at the indicated times.

Gaussian Preference Validation

To increase the loading of H2O molecules near the MgCl2·nH2O spherical clusters, a Gaussian selection preference is used. This preference gives a higher selection probability for inserted H2O molecules near the center of the box compared the outer vacuum region. In our systems, this center region is important because we use MgCl2·nH2O spherical clusters that are positioned at the center of the simulation box. This algorithm was first validated with an empty box at 400 K at 1 atm. These results are shown in the Supporting Information. It shows that more insertion trial moves were selected near the center. Furthermore, the Gaussian preference selection bias for the center of the box, does not change the total loading of the box, where it correctly predicts a H2O vapor density close to the NIST ref (44) value.

In Figure , the Gaussian preference is tested for a MgCl2 spherical cluster in the center of the simulation box. It is shown that when the Gaussian distribution is used to increase loading in the center of the box, near the MgCl2 spherical cluster, a marginal gain is achieved when the amount of water molecules is still low and as a consequence, the MgCl2 spherical cluster is relatively small. However, if the amplitude of the Gaussian distribution a is set high, it has a contrary effect. This is also the case when the Gaussian preference method is used without the WCA–GCMC algorithm since the possibility of atomic overlap is much higher in the center of the box. Due to the small gain with relatively small clusters, which could already provide relevant information regarding the hydration of MgCl2, the remaining calculations are performed including this Gaussian bias.

Figure 2

Gaussian preference test for a MgCl2 spherical cluster at 300 K with a (kcal/mol) as the height of the Gaussian distribution.

Hydration of MgCl2·6H2O Clusters

The motivation of this work was to develop a tool that can provide insights into the hydration mechanism of complex (new) materials for thermochemical heat-storage applications and predict equilibrium conditions at a given pressure and temperature. Therefore, the deliquescence equilibrium and behavior are studied by means of a spherical MgCl2·6H2O cluster of approximately 50 Å. In terms of hygroscopic salts as TCMs, this deliquescence could seriously affect the performance of the storage system. The experimental deliquescence equilibrium[45] line for MgCl2 is given by the black solid line in Figure . Under conditions below this deliquescence line, solid MgCl2·nH2O crystals occur, with n = 6 as the highest hydrated close to the deliquescence equilibrium. Under conditions above this deliquescence line, MgCl2 in an aqueous solution occurs.

Figure 4

MgCl2·6H2O deliquescence equilibrium, black solid line from experimental reference,[45] and symbols of GCMC prediction (red▲ = hydration, ● stable, and blue▼ = dehydration) as shown in Figure and in the Supporting Information.

In the GCMC simulation, the MgCl2·6H2O spherical cluster is placed at the center of a vacuum box and hydrated at vapor pressures of 6, 12, and 50 mbar, which are pressures around design conditions for thermochemical heat-storage systems with domestic applications.[8] In Figure , (de)hydration trends from the WCA–ReaxFF–GCMC algorithm are shown for the 12 mbar systems. Upward triangles red▲ represent increasing trends (hydration), and downward triangles blue▼ represent decreasing trends (dehydration). Black ● are given for stable trends (nondehydrating and nonhydrating). The results for the 6 and 50 mbar vapor pressures are given in the Supporting Information. These obtained trend symbols are plotted in Figure and show a close match with the experimental deliquescence equilibrium,[45] where the MgCl2·6H2O spherical cluster hydrates above the experimental equilibrium line and dehydrates under conditions below this line. The WCA–ReaxFF–GCMC results, given in Figure , predict a slightly higher equilibrium temperature (∼10 K) at lower vapor pressure. However, one must note that for the GCMC systems, steps of 10 K are used and it would be ambiguous to determine hydration or dehydration with smaller temperature steps (Figure ). Furthermore, a higher equilibrium temperature is expected for microparticles, since they have a higher solubility than the bulk material.[46] This is described by the Ostwald–Freundlich equation and caused by the relatively large factor of the surface energy for microparticles compared to the bulk material. As a result, a shift of the equilibrium curve to higher temperatures/lower vapor pressure will be present.

Figure 3

WCA–ReaxFF–GCMC results for a MgCl2·6H2O spherical cluster at p = 12 mbar, with k = 20 trial positions and a = 5 kcal/mol, where n on the y-axis represents the hydration level (MgCl2·nH2O) and the number of MC moves is given on the x-axis.

Hydration of MgCl2·nH2O Clusters

Thermochemical heat-storage systems for domestic heating are typically studied at a vapor pressure of 12 mbar, over a temperature range from 300 to 500 K.[8] In this sense, the MgCl2·nH2O clusters are studied under these conditions. Since the GCMC algorithm is computationally demanding, multiple starting structures are used at different relevant hydration levels (n = 0, 2, and 6), making the prediction of equilibrated hydration levels easier to estimate. These results are given in Figure . The corresponding estimated equilibrium hydration levels are compared with theoretical values, predicted by thermodynamics, in Figure . The computed equilibrium curves from thermodynamics, in which the left-hand side of chemical reaction eq is in equilibrium with the right-hand side, are described by the thermodynamic relation[7]in which ΔH0 and ΔS0 are the reaction enthalpy and entropy, respectively, per mole of water under standard conditions for the different components on the left-hand side and right-hand side of equilibrium 1. R is the universal gas constant, p0 is the standard pressure, and peq is the resulting vapor equilibrium pressure at a given temperature T. The exponential term in this equation makes the direct prediction of the equilibrium from entropy and enthalpy calculations (e.g., with DFT) hard, since small energy deviations can result in large equilibrium deviations. In Figure , the equilibrium curves are given for MgCl2·nH2O (n = 0, 1, 2, 4, and 6), in which ΔH0 and ΔS0 are taken from the NBS tables.[47]

Figure 5

Hydration of MgCl2·nH2O clusters at different temperatures (different colored lines) and different initial hydration levels (n = 0, 2, and 6) and a vapor pressure of 12 mbar. The insets visualize the highest hydrated cluster at 300 K with the initial MgCl2·6H2O cluster and the lowest hydrated structure at 500 K with the initial MgCl2 clusters.

Figure 6

Phase diagram of MgCl2·nH2O n = 0, 1, 2, 4, and 6, red lines are computed equilibrium lines by eq and NBS thermodynamic values[47] and the gray line is the experimental MgCl2·6H2O equilibrium line.[45] Blue crosses represent GCMC simulation conditions (T, p) of Figure and are given with the estimated water loading from the simulations. The black solid circles are the estimated MgCl2·6H2O deliquescence equilibrium from Figure . The dotted line is the water saturation equilibrium line. The values of the data points are given in the Supporting Information.

From Section , it was already shown that at 300 K and 12 mbar, MgCl2 would go to a higher hydration than the hexahydrate (n = 6). In Figure , this is observed for all simulated systems at different initial hydration settings (n = 6, 2, and 0). However, after 300,000 MC moves, the systems which started at a hydration level n of 0 and 2 are still far from equilibrium. At 340 K, an estimated hydration level between 2 and 4 appears; from the thermodynamic values, the tetrahydrate (n = 4) crystalline structure is expected. At 380 and 400 K, the GCMC simulations equilibrate approximately around a hydration level of 2, where the thermodynamic dihydrate (n = 2) can be found between 369 and 390 K, from eq , at the given vapor pressure. At 500 K, the GCMC algorithm predicts a hydration level of 1, where the thermodynamic rules indicate anhydrous MgCl2 (n = 0). The deviations between the results from the GCMC algorithm and the calculated equilibrium lines from thermodynamic values (eq ) can be explained by the use of disordered microparticles (23–50 Å) in the simulations, versus the used crystalline bulk thermodynamic values (H0, S0). This also indicates that it would require a much higher temperature to obtain a completely dry microparticle, compared to the bulk material.

From Section , it was shown that the water layer on top of the salt crystal breaks the crystal itself and forms a disordered structure. This effect reappears in the hydration of the MgCl2 clusters, where the added water creates a disordered MgCl2·nH2O cluster that allows further hydration within the created voids of the disordered structure. This is according to the hypothesized theory of a two-step hydration process by Sögütoglu et al.[6,7] in which the complete hydration process of a salt is described by two distinctive steps; step (1) water adsorption from the atmosphere to a wetting layer on the surface of the salt and accordingly dissolution of ions; step (2) nucleation into the crystal of the final hydrate. This second step, the nucleation of a crystal, is a rare event which occurs on time scales far beyond MD time scales.[48] Hence, this second step is not reached by or observed in the GCMC hydration modeling. Due to the absence of nucleation in crystalline structures, smooth transitions are expected from the GCMC simulations, compared to the distinctive zones from thermodynamic calculations. To visualize the hydration, two movies of the MgCl2 cluster, at 300 K and 12 mbar vapor pressure, are added to the Supporting Information. One movie shows the hydration of the particle from the first MC cycle until 250,000 cycles. The second movie shows a section at the center of the same particle in the yz plane, from the first MC cycle until 90,000 cycles. The initial wetting on the surface from the particle, the disordered development of the initial crystalline starting structures, and loading within the particle can be observed from these movies.

Conclusions

The motivation of this work was to develop a molecular modeling tool that can be used to study the hydration of salt hydrates for thermochemical storage applications. In this sense, an efficient ReaxFF–GCMC model has been designed and implemented with a force field for MgCl2·nH2O.[13] To increase accepted water loading in dense salt hydrate structures, the ReaxFF–GCMC was combined with two forward bias methodologies, namely, scanning of k insertion trials with the computationally cheap WCA potential and Gaussian preference selection of trials near the salt hydrate cluster. The WCA forward bias significantly increased the speed of the ReaxFF–GCMC algorithm, where it reached a much higher loading in a substantially shorter time. The Gaussian preference selection near the salt hydrate cluster showed a minor improvement, for relatively small clusters.

The biased WCA–ReaxFF–GCMC algorithm was used to study the hydration of MgCl2·nH2O microparticles. The predicted results are in good agreement with the experimental deliquescence line of MgCl2·6H2O, for a vapor pressure range from 6 to 50 mbar. Furthermore, different hydration levels were found over the temperature range from 300 to 500 K, from the deliquescence phase to the monohydrate (n = 1), respectively. The trend of the different predicted hydration levels is in agreement with thermodynamically computed values. However, with a slight deviation which can be attributed to the use of microparticles compared to the thermodynamic values based on bulk crystalline hydrates. The studied hydrated structures are found to be in a disordered state, where nucleation to a given hydrated crystalline structure is beyond the time scale of MD.[48]

From the hydration process, observed from the WCA–ReaxFF–GCMC simulations, it was shown that H2O molecules first load on the surface of the salt. The presence of water at the surface creates a disordered MgCl2 region; consequently, H2O molecules are also loaded within the salt itself. This fluidized nucleation process is according to the two-step hydration process formulated by Sögütoglu et al.[6,7]

The given results from the biased ReaxFF–GCMC algorithm validate the potential of the method for future studies. For example, to use it for complex TCM (e.g., encapsulated, impregnated, double, doped, or mixed salts) for which thermodynamic values are not available in the literature. Additionally, not only the hydration of TCMs can be studied but also other sorption processes (e.g., oxidation, hydration, and carbonation) that involve larger sorbate gas molecules and dense sorbent material.

Methodology

GCMC

To study molecular properties, one could use statistical thermodynamic rules combined with an ensemble of atomic positions. The ensemble must contain all relevant states of a system and thereby resemble its entire phase space. MD can obtain such an ensemble, in which successive molecular states are sampled by trajectories over time via integration of Newton’s laws. Alternatively, one could sample relevant states using MC methods, in which the states are generated according to the imposed probability distribution. Contrary to MD, states in MC are not sampled over time. This is the strength of MC, where one is not limited by time and/or slow diffusion to pass high energy barriers. For practical efficiency, new states in the MC ensemble are generated from the previous states—generating a Markov Chain—by performing trial moves which are accepted or rejected based on the probability of finding the system in a certain state and the probability of attempting the trial move.[29,49]

The GCMC algorithm is the most commonly used method to predict sorption phase equilibria with numerical modeling. The algorithm relies on three basic trial moves that are successively and randomly selected: translation of molecules, insertion of molecules, and deletion of molecules. The first one, translation (thermalization) of molecules, can be a rotation, translation, sampling the inner degrees of a molecule, or a combination of these. To improve efficient sampling of translational moves for complex systems, multiple advanced MC algorithms are proposed, for example, the AVBMC method.[50] In this work, due to practical ReaxFF implementations, we choose to perform a short MD simulation in the canonical ensemble. With this approach, the reaction will follow from the dynamics when the system is in the required chemical environment. The second one, the insertion trial move, is the insertion of a molecule within the system at a random location and random orientation. The inserted molecule is an “ideal” gas molecule taken from an infinite large reservoir.[29,38,49] The trial move is accepted with probabilityin which β is the reciprocal of the thermodynamic temperature 1/kBT with kB as the Boltzmann constant and T as the absolute temperature and f is the gas-phase fugacity computed from the pressure by the Peng–Robinson equation of state.[51] The critical temperature of 647.3 K, a critical pressure of 221.2 bar, and an acentric factor of 0.344 were used for water.[52] However, under the given conditions encountered in this work, fugacity and pressure are nearly identical as the pressure is low, V is the volume of the simulation box, Eig is the intramolecular energy of the isolated molecule, and ΔE is the change in energy of the system. The change in energy of the system is given by ΔE = Enew – Eold where Enew is the energy of the new state of the system and Eold is the energy of the old state of the system. The third trial move of the GCMC algorithm is the deletion of a randomly selected molecule from the system. This trial move is accepted according toin which ΔE = Enew – Eold, and with Eig as the intramolecular energy of the molecule in the conformation as in the system. Visualization of the GCMC algorithm is presented in Figure , in which the blue area is the newly developed biased algorithm as explained in the following sections.

Figure 7

Visualization of the GCMC algorithm, including the newly developed biased algorithm given in blue. The conventional basic GCMC moves are (1) translations, (2) insertions, and (3) deletions.

Biased ReaxFF–GCMC

In the GCMC algorithm, new trial configurations are created based on old states and the corresponding acceptance rules (eqs and 4). These rules are created based on the detailed balance condition[29] and it is key that the probability of creating the trial move from the old configuration to the new configuration is equal to the reverse way. Consequently, the insertion of a new molecule should be done at a random location and with a random orientation. This works well for inserting and deleting relatively small molecules in systems that contain many voids (e.g., MOFs[53,54] and zeolites[55]). However, this will be difficult for dense systems with very little suitable locations to achieve successful insertions. When the used interaction potential to compute the new energy is an expensive calculation, this could lead to very long computational times before equilibrium is reached. Hence, it makes sense to bias the position of inserted molecules toward more feasible positions.

WCA–ReaxFF–GCMC

We developed an advanced ReaxFF–GCMC method to avoid many unnecessary calculations with the ReaxFF formalism. In the case of an insertion trial move, our advanced ReaxFF–GCMC method first generates k trial positions that are evaluated with a computationally cheap WCA interaction potential (EWCA). The details of WCA interaction potential are explained in Section . Accordingly, from these k trials, one trial position is selected proportional to its normalized Rosenbluth factor[29,38,56]in which P is the probability of selecting trial i. Ebias, is the energy given by the computationally cheap WCA interaction potential (Ebias, = EWCA,) for trial position i in the new configuration, which is normalized with the total Rosenbluth weight

The energy of the selected trial position (k) is recalculated with the more expensive ReaxFF formalism. Due to the applied forward bias in selecting promising trials, the acceptance probability changes towhere ERxFF is the energy difference between the new and old configurations by ReaxFF. For the deletion of a molecule from the system, the Rosenbluth weight of the old configuration has to be computed, which is done by performing k – 1 random insertion trial positions, and the kth position is the old configuration itselfin which Ebias, is the energy of the kth position, the selected molecule for deletion, computed by the cheap WCA interaction potential. The modified acceptance rule is given by

In the Supporting Information, it is shown that eqs and 9 obey the detailed balance condition.

WCA–ReaxFF–GCMC, with Center Preference

To study the hydration of MgCl2, an anhydrous MgCl2 spherical cluster is placed in the center of a larger simulation box. This allows hydration on all the different surfaces of the anhydrous cluster. Since the cluster is placed in an empty simulation box, the possibility of sampling the empty space around the cluster is large compared to the dense region in and around the MgCl2 cluster. To increase selection of promising trial positions near the MgCl2 cluster, an extra biasing potential is applied in the form of a Gaussian energy distribution with its center at the center of the cluster, as illustrated in Figure .

Figure 8

Example of preference selection by a Gaussian energy distribution. The x-axis is the relative x coordinate within L box size, the energy from the Gaussian distribution is in blue on the left y-axis, and in red on the right y-axis is the proportional selection probability from the Gaussian bias energy.

This added Gaussian energy distribution is described byin which a is the height of the Gaussian distribution (in [kcal/mol], the same energy units as EReaxFF and EWCA), c is the width of the distribution, x, y, and z are the coordinates of the trial position, and x, y, and z are the center of the cluster. The energy term EGauss is added to the energy of the WCA interaction potential (Ebias = EWCA + EGauss). This modification results in the fact that if a successful trial is found in the center of the cluster, the probability of selecting this one for insertion is higher than a successful trial position further away from the center of the simulation box.

Force Fields

As explained in the previous section, k trial insertions are performed for each GCMC trial move with a relatively cheap WCA potential followed by the more expensive ReaxFF calculation.

ReaxFF

ReaxFF MD[25−27] enables simulations including reactions. The ReaxFF interaction potential is a summation of different energy contributions

The terms EvdW and ECoul are the noncovalent bonded van der Waals and Coulomb terms, respectively. The Ebond term accounts for the covalently bonded atoms. The terms Eval, Etors, Epen, Eunder, Eover, and Econj describe the valence and torsion contributions, “penalty” energies, under- and over-coordination, and conjugated systems, respectively. Eothers can include other terms for specific systems, such as H bonds or extra dispersion interactions.[57] The bond order (BO) between atoms is described by a summation of empirical relations for the BOσ, BOπ, and BOππ bond, which depend on the distance r between the atoms i and j.[25]in which r0σ, r0π, and r0ππ are the bond radii for the σ, π, and ππ bond, respectively. The pbo values are fitted parameters to experimental or first-principles results. Each BO term has a maximum value of 1, and when all BOs contribute, BO could add up to 3. The empirical BO approach of ReaxFF enables us to model bond breaking and formation without expensive quantum mechanical calculation.

To study the hydration of MgCl2 hydrates, we used the ReaxFF force field developed by Pathak et al.[12,13] in combination with the long-range corrected H2O ReaxFF force field.[57] This force field was not explicitly trained to recreate the phase diagram of the salt; however, a transferable force field for TCM application should be able to recreate it. Furthermore, the MgCl2 force field has proved itself useful for multiple hydration levels MgCl·nH2O (n = 0, 1, 2, 4, and 6), and the H2O force field is able to accurately capture both condensed and vapor phase under saturation conditions.[57]

WCA Potential

Prior to the ReaxFF insertion, k computationally inexpensive WCA potential trial insertions are computed. This WCA[42] interaction potential is relatively cheap because it only considers short-range repulsive interactions. In this way, hard overlaps between atoms are avoided. The WCA interaction potential is described bywith ϵ and σ as characteristic energy and distance parameters, respectively.

In this work, we used the WCA parameters, as given in Table , combined with the Lorentz–Berthelot (LB) mixing rules between different elements. Careful determination of the WCA σ and ϵ parameters is key to use the WCA bias effectively. If these parameters are too large, promising small voids in the dense MgCl2 structure will never be selected. However, when they are too small, many unrealistic insertions with overlapping atoms can still be selected, and the algorithm becomes ineffective again. The parameter σ corresponds to the atomic/molecule size. For the H2O molecule, only the oxygen atom was considered in the WCA potential, with a diameter (σ) corresponding to the radial distribution function (RDF) O–O distance for liquid water obtained by the ReaxFF force field. Thereby, the oxygen in the WCA calculation resembles the diameter of a H2O molecule in liquid, and excessive WCA interaction calculations with twice as much H atoms are avoided. The magnesium and chlorine parameters were chosen such that the O–Mg and O–Cl repulsive part closely matches the repulsive distance of the ReaxFF force field. In the Supporting Information, the resulting O–O interaction potential is compared with the ReaxFF RDF, and the O–Mg and O–Cl interaction potentials are compared with the ReaxFF interaction energy.

Table 1

WCA Parameters

element	ϵ/[kcal/mol]	σ [Å]
O	0.1	2.3
Mg	0.1	0.8
Cl	0.1	3