Molecular Force Field Development for Aqueous Electrolytes: 2. Polarizable Models Incorporating Crystalline Chemical Potential, and their Accurate Simulations of Halite, Hydrohalite, Aqueous Solutions of NaCl, and Solubility

The current state-of-the-art force ﬁelds (FFs) for Na + a Cl − ions are not capable of simultaneously predicting the thermodynamic properties of the aqueous solution and the crystalline phase. This is primarily due to an oversimpliﬁcation of the interaction models used but partially also due to the insuﬃcient parametrization of the FFs. We have devised a straightforward and simple parametrization procedure for determining the ion-ion interaction parameters in complex molecular models of NaCl electrolytes which involves ﬁtting the density, lattice energy and chemical potential of crystalline NaCl at ambient conditions. Starting from the AH/BK3 and MAH/BK3 FFs, the parametrization approach is employed to develop a complex and accurate polarizable molecular model for the NaCl electrolyte by parametrizing the ion-ion interactions. The performance of the reﬁned polarizable NaCl FF is assessed by evaluating the dif-ferent thermodynamic and mechanical properties of the crystal, density of crystalline and molten NaCl along with the melting temperature, properties of aqueous solutions, and the structure and stability of hydrohalite. The simulation results conﬁrm the superiority of the reﬁned FF in comparison with the existing state-of-the-art FFs to accurately predict a wide range of system properties in diﬀerent NaCl phases, including NaCl aqueous solubility. The reﬁned FF may ﬁnd applications in the accurate simulations of NaCl electrolytes including inhomogeneous environment, phase equilibria and interfaces, and metastable states. Finally, the parametrization strategy is robust and general, and can be used to devise molecular models for other electrolytes.


INTRODUCTION
Sodium chloride and other alkali-halide electrolytes are important in many common, industrial, and biological processes, from simple cooking to advanced power storage or nuclear power technologies. Understanding such processes involves understanding interactions between the Na + and Cl − ions and also between the ions and the molecules of water. These interactions cannot be probed directly by experiments due to their microscopic nature. Computer modeling is thus used to develop molecular interaction models, force fields (FFs), that can be used to simulate the microscopic behavior of both liquid and solid phases, and to predict their thermodynamic, structural and dynamical properties. The accuracy of such computer simulations inevitably depends on the accuracy of the FFs employed, and it reflects our understanding of electrolytes, their solutions, and their hydrates at the microscopic level.
There is a plethora of FFs for Na + and Cl − ions in literature, e.g., see our reviews, 1,2 and new models continue to be developed. [3][4][5][6][7][8] It is thus surprising that there is no FF yielding accurate predictions of thermodynamic properties of both the NaCl aqueous solutions and NaCl crystal (halite). It is common that an electrolyte FF completely fails in predicting some of the key properties. For example, the complex polarizable AH/BK3 FF predicts six times lower aqueous solubility 9 as compared to experiment, the recent Madrid FF misses the NaCl crystal lattice energy and also the chemical potential both in the solution and solid phases by ∼ 182 kJ mol −1 , 5 the NaCl/ε FF 3 predicts incorrect concentration dependence of activity coefficent and very low solubility, 10 and 12 of 13 electrolyte FFs assessed in our previous study 11 fail to predict of the density of NaCl crystal. One can attribute these problems to oversimplified or incorrect functional forms of interactions employed in the molecular models, e.g., using the Lennard-Jones potential or simple combining rules for short-ranged repulsive forces, neglecting polarizability, considering point distributions of charges, using inaccurate water FFs, and employing charge scaling without applying the electronic continuum corrections. In addition, these problems can also be related to insufficient parametrization of the FFs, often targeting a limited set of properties and excluding properties that are important for the correct overall behavior of the molecular models; e.g., targeting only properties of infinitely diluted solutions, neglecting energetic and dynamical properties, and phase behavior, and targeting non-reliable experimental data.
In our previous study, 12 we attempted to develop a new set of parameters for a simple nonpolarizable model based on Lennard-Jones and point charge interactions accompanied by Lorentz-Berthelot rules for NaCl aqueous solutions at ambient conditions. We showed that the model was not able to simultaneously fit the dependence of the chemical potential and density on NaCl concentration over the entire experimentally accessible concentration range, and we thus concluded that more complex models are needed when the concentration dependence of several properties are of interest. The inadequacy of simple nonpolarizable models for electrolyte solutions is also supported by other studies. 11 For example, it is well known that all nonpolarizable FFs yield too slow dynamics in concentrated solutions unless their charges are scaled down. However, the decrease of charges results in an enormous increase of potential energy. 5 Furthermore, the simulation studies of electrolytes typically deal with spatially heterogeneous systems, e.g., electrolytes in confined systems such as mineral pores, [13][14][15] at solid surfaces with the surface polarizability such as graphene, 16 where molecular polarizability plays very important roles compared to bulk liquids. Hence, polarizable models should be favored in principle.
The AH/BK3 FF 17 turned out to be a promising candidate for an advanced model of aqueous electrolyte solutions because it employs complex interaction potentials while also being tractable for a wide range of applications. It considers polarizability, Gaussian distributions of charges, and the Buckingham intermolecular potential. Moreover, the AH/BK3 FF relies on the polarizable BK3 water model 18 which predicts the majority of the properties of pure water extremely well, e.g., density, self-diffusivity, viscosity, compressibility, permittivity, and chemical potential. 9 The BK3 FF has its origin in the GCP model of Cummings and coworkers. 19 The AH/BK3 FF for NaCl yields correct concentration dependence of density and activity at ambient conditions, and its chemical potential in NaCl aqueous solutions is lower by only ∼ 3.35 kJ mol −1 compared to experiment. However, the AH/BK3 FF has a few drawbacks as it exhibits very strong Na + -Cl − ion pairing in aqueous solutions and too low chemical potential of NaCl crystal by ∼ 15 kJ mol −1 , resulting in unrealistically low solubility, 9 and its temperature dependence of activity coefficients is incorrect. 20 The deficiencies of the AH/BK3 FF are related to its parametrization as its authors did not target any properties of concentrated solutions, and they considered the properties of crystalline NaCl only loosely; e.g., missing the density of NaCl crystal by ∼ 45 kg m −3 and its lattice energy by ∼ 9.3 kJ mol −1 . 17 Later, the parameters of the ion-ion interactions were refitted by Kolafa, yielding the MAH/BK3 FF 4 for NaCl. The MAH/BK3 FF successfully predicts several mechanical properties of the NaCl crystalline phase at ambient conditions, lattice energy, and additional properties at the melting point. However, Kolafa's explicit interface simulations indicate that the prediction of the NaCl aqueous solubility by the MAH/BK3 FF is too low as compared to experiment by ∼ 40%. Our simulations confirmed this low aqueous solubility and showed this is due to a low value of the chemical potential of NaCl crystal by the MAH/BK3 FF when compared to experiment by ∼ 4.4 kJ mol −1 . It is puzzling that even the complex polarizable NaCl MAH/BK3 FF, targeting many properties of crystalline phase and including experimental lattice energy, yields too low a value of the chemical potential and aqueous solubility when compared with experiments.
Hence, the aims of this work are the following: (i) to investigate, if and how the NaCl MAH/BK3 FF can be refined to yield correct predictions of the aqueous solubility and chemical potential of the crystal while maintaining accurate predictions of other system properties; (ii) to determine new Na + and Cl − FF parameters; (iii) to test the refined NaCl FFs under different conditions, including NaCl aqueous solutions, crystalline NaCl (halite), NaCl·2H 2 O low temperature crystalline phase (hydrohalite), and NaCl melt; and (iv) to provide insight into the microscopic behavior of NaCl systems and suggestions for possible improvements of FFs for modeling of other alkali-halide electrolytes.
The remainder of the paper is organized as follows: In section 2, we propose a method for fitting the FF parameters of Na + and Cl − ions, targeting crystalline chemical potential and pressure at specified density and temperature. In section 3, we describe the simulation methodology employed. In section 4, we present our refined NaCl FF and assess its performance at various state conditions and phases. Section 5 provides our conclusions, including suggestions for the possible improvement or development of FFs for other electrolytes.

FORCE-FIELDS TARGETING CHEMICAL POTEN-TIAL AND DENSITY OF HALITE
We argue that, besides different properties of the solution phase, development of any electrolyte FF should always include the chemical potential of crystalline electrolyte, µ cr , within the set of target properties. Particularly, excluding µ cr from the target properties often results in an inaccurate prediction of aqueous solubility, which limits the usability of the FF and may lead to precipitation at low concentrations or spurious behavior in general. We are not aware that considering µ cr as the target property has been applied in developing polarizable FFs, which is not surprising as its implementation poses a serious challenge.
Calculation of chemical potentials by molecular simulations is a challenging task, always involving tens of separate simulations at different integration stages, in systems of different sizes, and extrapolations to infinite systems. Developing a complex polarizable FF targeting µ cr thus relies on a fast converging process that needs to avoid any iterative procedure involving a long sequence of iteration steps; e.g, Newton's method used by Kolafa comprised of 25 iteration steps. 4 Such iterative approaches would be very time consuming if µ cr (typically also its partial derivatives with respect to all the fitted FF parameters are required) was calculated within each of the numerous steps. Moreover, such iterative approaches do not provide further insight into the microscopic behavior of crystalline NaCl. Hence, we have devised a fast converging procedure based on several physical assumptions, which allows us to straightforwardly and easily reparametrize the MAH/BK3 FF to yield correct values of µ cr and in addition, also pressure P cr at the given density of crystalline NaCl, ρ cr , and temperature, T . The fitting procedure can be further extended to target also other properties or properties at different thermodynamic conditions.

Force-Field Functional Form
We consider the same functional form of the molecular interactions as the AH/BK3, MAH/BK3, and BK3 FFs, which is briefly summarized below. The AH/BK3 and MAH/BK3 FFs model the nonpolar interactions between ions by the Buckingham (EXP6) potential where subscripts i and j refer to the particles. The electrostatic interaction between ions is modeled by two Gaussian charge distributions per each ion. One of the Gaussian charges is fixed to the center of the ion. The other Gaussian charge is massless and is attached to the center of the ion by a virtual harmonic spring, modeling polarizability by the Drude chargeon-spring approach, where the spring constants are obtained from the particle polarizability, α ia . The interactions between Gaussian charges of different particles converge to the Coulomb potential at long distances but become weaker at short distances due to the overlap of the Gaussians depending on their distribution widths, σ ia : where subscripts a and b refer to the interaction sites of the particles, q ia are the values of charges, ε 0 is the permittivity in vacuum, and the summation over a, b runs over the two charge sites.

Unchanged Force-Field Parameters
To achieve the desired values of target properties without altering other properties, we aim to minimally modify the MAH/BK3 FF, keeping most of the interaction parameters of the MAH/BK3 and AH/BK3 FFs unchanged. We preserve the total charge magnitudes as there is no reason to scale charges in polarizable models. We also did not modify polarizability, α ia , as the values of α ia used in the AH/BK3 FF are plausible. We further kept unchanged all other electrostatic parameters, q at ia , q Dr ia , and σ ia because their roles are rather technical, and they affect the behavior of crystalline NaCl only marginally (e.g., Kiss and Baranyai stated that the size of the charge-on-spring is arbitrary, apart from very extreme choices. 18 ), or their changes can be compensated by changes in the short-ranged repulsive exponential interactions. Finally, we did not change the values of C as we argue that the London's attractions at long distances are negligible compared to strong electrostatic interactions, and they can be compensated by changing the repulsive exponential interactions at short distances. We note that these are just assumptions whose validity will be judged by the parametrization procedure and performance of the refined FF.

Reparametrized Force-Field Parameters
We decided to modify the exponential repulsion interactions of the ion-ion interactions by changing the values of parameters A and B for the ions (a similar approach was adopted by Kolafa 4 ). We argue that the exponential repulsion interactions play the key role and they balance strong electrostatic attractions between unlike ions while simultaneously keeping realistic interactions between like ions. We further assume that the exponential repulsion between Na cations is negligible due to the small size of Na + and typical Na + -Na + distances in the NaCl aqueous solutions and crystals. On the other hand, the exponential repulsion between Cl anions may play a role as Cl − is relatively large with respect to Na + . Hence, we focus on the Na + -Cl − and Cl − -Cl − interactions and obtain the Na + -Na + A and B parameters from the Kong rules.

Relation between Short-Ranged Repulsion Potentials and Crystalline Properties
We argue that small changes in the exponential repulsion between ions in crystalline NaCl significantly affect the energetic properties and equilibrium density of the crystal phase.
Considering a canonical ensemble (given temperature, T ; volume V ; number of ion pairs in the crystal, N NaCl ), the lattice energy, U cr , and pressure, P cr , include contributions due to the exponential repulsion, which can be expressed using radial distribution functions, g ij (r), In Eqs. (6) and (7), the sums run over two types of ions, and are the exponential repulsion potential between the ions and its derivative. Considering the face-centered cubic lattice of the NaCl crystal, one can devise simple formulae which approximate the exponential contributions, U EXP cr and P EXP cr , as where the values r = 2.775 and 3.975Å correspond to the locations of first maxima in g NaCl (r) and g ClCl (r) as obtained from a molecular simulation using the MAH/BK3 model at an experimental ambient crystal density 2163 kg m −3 and temperature 298.15 K.

Fitting Crystalline Density and Lattice Energy
Eqs. (10) and (11) are employed in the parametrization procedure as follows. We start with the interaction parameters of the MAH/BK3 FF and simulate the lattice energy, U cr , and pressure, P cr , at ρ = 2163 kg m −3 and T = 298. 15 (10) and (11). Values of A NaNa and B NaNa are then obtained via the Kong rules, Eqs. (3) and (4). However, as we fit six parameters, we need six equations. Fitting two properties (U cr and P cr at given ρ and T ) yields two equations, the Kong rules are another two equations, and we thus need two additional conditions. Keeping the minimal change in the MAH/BK3 FF, we opt to fix the fractions of the energetic and pressure contributions from the Na + -Cl − and Cl − -Cl − interactions for the MAH/BK3 FF, i.e., If necessary, these two conditions can be modified. For example, one can consider different values of θ U and θ P or can substitute the conditions by other equations similar to Eqs. (10) and (11), corresponding to additional target values of P cr and U cr at different density and temperature and/or can employ other approximate equations related to different target properties. Our results show that such changes are not necessary unless very accurate predictions of the crystal or melt properties at high temperatures are required.
However due to the approximate nature of Eqs. (10) and (11), the fitting procedure does not directly give a perfect agreement between the simulation and target properties and it was necessary to repeat the parametrization procedure a few times. However, this is a rather easy task as the procedure is fast and quickly converging. We found a very good agreement between the simulation and target values of U cr and P cr within three iterations.
We parametrized our refined FF at ρ = 2163 kg m −3 and T = 298.15 K to the target value of P cr = (1 ± 10) bar. Due to the scattered values of the experimental lattice energy found in literature, we repeated the fitting procedure for several different target values of U cr , and we discuss the effects of the choice of the different values for experimental lattice energy below.

Fitting Crystalline Chemical Potential via Lattice Energy
The molar chemical potential of NaCl in the crystal, µ cr , is the molar Gibbs free energy and it consists of U cr , the FF-independent kinetic energy, the entropic contribution −T S, and the (negligible) pressure contribution P V in one mole of NaCl. The structure of the NaCl crystal in a canonical ensemble only weakly depends on the FF parameters. Therefore, we assume that the same also applies to the entropy and the dependence of µ cr on U cr becomes approximately linear for different FF parametrizations. Fig. 1 displays the actual dependence of µ cr on U cr obtained by molecular simulations with different FF parametrizations, i.e., FFs fitted to ambient pressure and different values of U cr . Fig. 1 suggests that one can find two FF parametrizations targeting one P cr and two values of U cr at a given density and temperature, obtain their µ cr by simulations, and subsequently find a value of U cr , which gives the desired value of µ cr by a linear interpolation.
We implemented this approach with the target U cr equal to -780 and -786 kJ mol where P (±e) is the pressure of the expanded/compressed crystal and e = 0.01 was used.
The values of G were calculated from a simulation with the crystal lattice rotated around the z-axis by 45 deg with respect to the simulation cell, with the system size expanded exp(e)-times in the x-direction and exp(−e)-times in the y-direction, where P xx and P yy are components of the pressure tensor and e = 0.00707. The values of α were evaluated using with ∆T = 30 K. In all the cases, the pressure and components of the pressure tensor were simulated using the volume-perturbation method with an optimal ∆V selected from six different values to avoid any systematic errors due to numerical differentiation. The densities of crystalline and molten NaCl were obtained from simulations in an NPT ensemble with 500 ion pairs and cubic simulation cells. The equivalent isotropic displacement parameters in crystalline NaCl, U eq , were calculated using the equation where u is the displacement vector of an ion from its equilibrium position.
The chemical potential of crystalline NaCl, µ cr , was determined using the Frenkel-Ladd

New Force-Field Parameters
We applied the fitting procedure described in section 2 to parametrize the ion-ion interactions and the resulting new interaction parameters are listed in Table 2   Crystalline NaCl at Ambient Conditions Table 3 lists several properties of crystalline NaCl obtained from simulations with our refined FF and compares them with our simulated values using other state-of-the-art models from literature, Madrid and MAH/BK3 FFs, and also with experiments. We see from Table 3 that the crystalline density, ρ cr , obtained at P = 1 bar and T = 298.15 K, for the refined and MAH/BK3 FFs agrees with the experimental value while the experimental ρ cr is slightly overestimated by the Madrid FF; by ∼ 8 kg m −3 which is rather small compared to many other FFs found in literature. 11 It is worth mentioning that experimental data are subject to uncertainties and possible lattice defects are neglected in simulations.
The chemical potential of crystalline NaCl, µ cr , reported in Table 3   The mechanical properties B, G and α obtained from our FF also agree rather well with experimental values, as seen in Table 3; differences between our and the experimental values are less than 5%. Our simulation results for B, G and α exhibit almost the same accuracy as those obtained using the MAH/BK3 model which was developed by directly targeting the mechanical properties of crystalline NaCl. The performance of the Madrid FF is rather poor with about 66%, 27%, and 23% deviations from experimental B, G, and α, respectively.
In addition for our refined FF, we evaluated equivalent isotropic displacement parameters, U eq , for both the ions. Our results, 0.0203Å 2 and 0.0178Å 2 respectively for Na + and Cl − ions, are in excellent agreement with experimental data, 38 0.0201(8)Å 2 and 0.0174(5)Å 2 , respectively.
Density of Crystalline and Molten NaCl, and Melting Temperature simulations result in the melting temperature of 1025(5) K, which is close to the experimental 1073.9 K. 34 We note that the explicit crystal-melt simulations may suffer from finite size and hysteresis effects. If desired, the high-temperature behavior of our FF can be improved, e.g., by considering high-temperature target properties instead of Eqs. (12) and (13). However, this would need to be performed with caution as the experimental high-temperature properties often suffer from large uncertainties.
Regarding comparison with other FFs in Fig. 2       First, we discuss comparison of results for the refined FF with experiments shown in Fig. 5. We see that µ cr simulation results (solid horizontal lines) are in an excellent agreement with experiment not only at ambient conditions (see Table 3  In Fig. 6d

CONCLUSIONS
We have shown that it is possible to straightforwardly and relatively easily reparametrize the ion-ion short range repulsion contributions of the AH/BK3 17 and MAH/BK3 4 NaCl FFs to yield experimental values of the crystal density, chemical potential, and lattice energy at ambient conditions. Our parametrization procedure employs several simple assumptions involving the radial distribution functions in the NaCl crystal and their approximate analytical relations to pressure and lattice energy, and another relation between the energy of the crystal and its chemical potential. These relations not only facilitated fast convergence of the iterative process within the parametrization procedure but also provided further insight into the microscopic behavior of the systems. The procedure can be generalized or further extended by additionally devised relations either based on statistical mechanics (e.g., exponential contributions to mechanical properties such as B, G and α), thermodynamics (e.g., targeting energy and pressure values at multiple state points to fit B and α indirectly), mechanics and geometry (e.g, substituting the constant r values in Eqs. (8) and (9) by geometrical rules to target crystals at different densities).
An example in which the procedure provides an important insight into crystalline electrolytes are the approximate relations for the exponential contributions to the lattice energy and crystal pressure, Eqs. (10) and (11). These relations show that simultaneous fitting of the density and energy of the crystal for simple molecular models based on Lennard-Jones (LJ) potential would fail. For the simple models, it is impossible to meaningfully change the derivative and value of the repulsive part of the LJ potential at particle separations corresponding to the nearest neighbors without unrealistically altering the London's attractions.
In contrast to the LJ potential, the Buckingham potential allows simultaneous fitting of both the density and energy of the crystal via a change of its interaction parameters A and B; cf.
The resulting NaCl FF can be viewed as a refinement of the MAH/BK3 FF of Kolafa 4 for the ion-ion interactions, which adopts the original AH/BK3 FF of Kiss and Baranyai 17 for the H 2 O-H 2 O and H 2 O-ion interactions. The refined FF provides excellent predictions of the system properties in aqueous solutions, anhydrous NaCl crystal and melt, and an accurate structure of hydrohalite. The performance of the refined FF is superior to the existing state-of-the-art FFs since no significant disagreement between the simulation and experimental results were found for a wide range of system properties, which does not apply to other FFs (The Madrid FF yields extremely high chemical potential and lattice energy, and incorrect mechanical properties of crystal; the AH/BK3 FF has very low solubility, low crystal density, very low crystal chemical potential, and very low lattice energy; the MAH/BK3 FF has low chemical potential of crystalline NaCl and low NaCl lattice energy).
Relatively small discrepancies between the simulation and experimental data were observed for the self-diffusivity in the aqueous solutions with errors comparable with those for other state-of-the-art FFs.
As regards general development or improvement of the molecular models for electrolytes, we provide several suggestions: • It is desirable to include the salt chemical potential in both the solution phase and crystal as a target property in the development of electrolyte FFs.
• The Buckingham potential is preferred to the LJ potential as the Buckingham potential allows targeting both the density and energy of the crystal while for the LJ potential, it leads to spurious London's attractions.
• Development of the electrolyte FFs requires critical assessment of existing experimental (target) data since they may be rather scattered or misleading, e.g., the mismatch in experimental values of NaCl lattice energy/enthalpy which affects the NaCl aqueous solubility by 40 %.
• Our results for water self-diffusivity are similar to those of the Madrid FF which uses the scaled charges. Although both the models improved prediction of the water selfdiffusivity when compared with nonpolarizable models they still exhibit faster decrease in the water self-diffusivity upon an increase of the salt concentration when compared with experiments. This indicates that the slower water self-diffusivity in the nonpolarizable models is not only caused by approximative treatment of polarizability but it may be also related to other features of the models, e.g., the TIP4P water geometry employed in both the BK3 and Madrid water models.
• Our simulations of hydrohalite suggest a way to parametrize electrolyte molecular models by targeting the properties of hydrates. The experimental microscopic properties of hydrates such as lattice constants or atomic displacement parameters are readily available, they are directly related to the interactions of ions with water molecules, and they can be obtained from molecular simulations relatively easily. This can find applications for example in the development of an accurate Li + FF where properties of different LiCl hydrates can be utilized. For