{'fraction_non_alphanumeric': 0.028959027999817293, 'fraction_numerical': 0.0136573333942356, 'mean_word_length': 4.871279163314561, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Transferability of polarizable models for ion-water electrostatic interaction

Studies of ion-water systems at condensed phase and at interfaces have pointed out that molecular and ionic polarization plays an important role for many phenomena ranging from hydrogen bond dynamics to water interfaces' structure. Classical and ab initio Molecular Dynamics simulations reveal that induced dipole moments at interfaces (e.g. air-water and water-protein) are usually high, hinting that polarizable models to be implemented in classical force fields should be very accurate in reproducing the electrostatic properties of the system. In this paper the electrostatic properties of three classical polarizable models for ion-water interaction are compared with ab initio results both at gas and condensed phase. For Li+- water and Cl--water dimers the reproducibility of total dipole moments obtained with high level quantum chemical calculations is studied; for the same ions in liquid water, Car-Parrinello Molecular Dynamics simulations are used to compute the time evolution of ionic and molecular dipole moments, which are compared with the classical models. The PD2-H2O model developed by the author and coworkers [Masia et al. J. Chem. Phys. 2004, 121, 7362] together with the gaussian intermolecular damping for ion-water interaction [Masia et al. J. Chem. Phys. 2005, 123, 164505] showed to be the fittest in reproducing the ab initio results from gas to condensed phase, allowing for force field transferability.


Introduction
It is widely accepted that the implementation of polarizable force fields in classical Molecular Dynamics (MD) simulations is important for reproducing with high reliability the features of inhomogeneous systems [1,2,3,4,5,6,7]. Even though dipolar interactions decay faster than coulomb interactions they are known to be responsible of many macroscopic properties of interfaces [8,9,10,11]. In particular, in water-air, water-protein and water-framework (in confined geometries) interfaces, the induced dipole moments vary over a wide range in different regions of the system and are responsible of many phenomena such as nanobubbles' coalescence [12,13]; in the bulk, ionic and molecular polarization influences solvation shell dynamics [14,15,16,17,18,19,20]. For these reasons it is required that the polarizable force fields implemented in MD simulations are able to faithfully reproduce both the energetic and electrostatic properties without diverging for high electric fields. Most of the studies on polarizable models focus only on the energetic properties; as usual in force field development, parameters are optimized to fit the ab initio energies of a large sample of configurations. Only recently it has been observed that, although common force fields fit with high accuracy the energetic properties, they perform poorly in reproducing the ab initio induced dipole moments [21]. This perspective guided the author in pursuing a method to develop polarizable models which could faithfully reproduce the electrostatic properties of simple systems such as ion-molecule (either water or carbon tetrachloride) dimers [21,22,23,24,25]. In particular the focus of previous papers was on the validation of polarizable models usually implemented in MD simulations: classical calculations for simple systems at gas phase were compared to high level quantum chemical calculations. The rationale behind this approach is that the electrostatic properties (arising from the interactions among point charges and point dipoles) should be in agreement with the ab initio calculations; this is already the standard way for developing intermolecular no-polarizable force fields. Point charges are obtained by a fit of the ab initio electrostatic potential with the additional constraint of reproducing the permanent dipole moment of the molecules [26,27,28,29].
In the case of polarizable models the approach proposed is to find a set of functions and parameters (vide infra) which fit the ab initio induced dipole moments for a certain number of sample configurations for dimers. Once a good polarizable model is found, the following step for the development of the complete force field should be, then, to parameterize short range interactions, what could be done with the common methods of force field parameterization (e.g. energy fitting of many sample configurations). Before that, the issue of the model portability should be addressed; in fact this feature is required for a successful implementation of a force field in the study of a wide range of systems from gas to condensed phase. Recently, many authors [30,31,23,24,25,32,33] have pointed out that one of the largest limitations of standard molecular-mechanics polarizable force fields lies in that the transferability of gas phase derived potentials to condensed phase is hindered by the absence of polarization exchangecoupling in classical models (since it is due to short range electron cloud repulsion, it is also known as Pauli effect [33]). The portability of water-ion polarizable models is the focus of this contribution. In previous works, the author and his coworkers showed that, among the most popular methods implemented in MD algorithms to include polarization, the induced Point Dipole method performs better than fluctuacting charges and Drude models. Moreover, in the case of ion-molecule electrostatic interaction, it has been shown that an intermolecular damping function should be added to account for short range electron repulsion which cause the effective polarizability to be lower. This effect is independent of the species polarizability and it is rather related to the Van der Waals radius; in fact it was found to be important both for ions with high (anions) and low polarizabilities (cations) [23,24,25]. The outline of this article is as follows: in section 2 we discuss the computational details and the methods used; results and conclusion are given respectively in sections 3 and 4.

Ab Initio Calculations
To explore the portability of the polarizable force field, the results of ab initio calculations for gas phase and condensed phase were compared with classical calculations. Reference ab initio data were obtained using different approaches depending on the size of the system.
In the case of ion-water dimers five different configurations were chosen; for each of them the ion-oxygen distance was varied within an interval of ∼ 5Å. The shortest distance was selected to be be equal to the distance where the potential energy was ≈ 25 kJ mol −1 (≈ 10k B T at standard temperature); this criterium is very strict since some points are very likely to be inaccessible in condensed phase [21,22,23]. This fact should considered when assessing the significance of the results at very short distances.
All quantum chemical calculations were performed with the commercial package Gaussian 03 [34]. The B3LYP/aug-cc-pvtz model chemistry was used for all atoms. The density functional used was choosen because it is known to perform well, with small errors for the computed dipole moments and polarizabilities [35]. In a previous work we estimated that for such systems B3LYP slightly overestimate CCSD calculations, the largest difference being ∼ 4% [24]. The total dipole moment of the dimer has been defined taking the position of the ion as the origin of the reference system. In this way, in the limit of zero polarizability for the ion the total dipole moment of the dimer is equal to that of the molecule. The pure quantum chemical approach for condensed phase system would not be feasible, being the size of the system too large; in fact, since the configuration space to sample would increase exponentially with the number of molecules, performing ab initio calculations for each point would require huge computational resources. To overcome this problem, the Car-Parrinello (CP) [36,37] scheme for propagating the wavefunctions and the configurations has been used. This allows to sample a great number of configurations that could occur in typical MD simulations. The system under study was formed by a single ion surrounded by 64 water molecules at ambient conditions. The BLYP/PW (Plane Waves) model chemistry [38,39] was used with Troullier-Martins pseudopotentials [40] with a 70 Ry cutoff. The time step was set to 4 a.u. and the fictitious mass for the orbital was chosen to be 400 a.m.u. Periodic boundary conditions were applied. Production runs of 30 ps followed equilibration runs of 10 ps. Every five time steps, the Wannier centers' [41,42,43] coordinates r j (t) were computed (see figure 1); in this way, given the ionic coordinates R i (t), the dipole moments of molecule (or ion) I were obtained as: Q i and q j being respectively the charge of the i th ion and of the j th Wannier center. Contrary to the case of ion-water dimers, the frame of reference chosen for charged species was the molecular (ionic) center of mass. The advantage of this approach as compared to the one used for dimers is that the contribution of each species to the ab initio dipole moment can be known. On the other hand the approach used for dimers gives high quality results (which are needed for a fine tuning of parameters) but are limited to the total dipole moment of the dimer.

Classical Calculations
For all the configurations sampled with ab initio methods, the molecular (and ionic) dipole moments were calculated with the Point Dipoles (PD) method [44,7]. Within this method, a value of polarizability α j is assigned to site j; it is not required that such sites are localized on the atoms, neither that they bear physically meaningful values of α j . The only requirement is that the polarizability tensorα I (which could be easily obtained from the site polarizabilities and their coordinates) of the whole molecule (ion) is comparable to the experimental one. The  ions were assigned their polarizability at gas phase, while water molecules were described with three different models: (i) PDM-H2O: one single polarizable site "M" is located at a distance d = 0.215Å from the oxygen along the permanent dipole moment vector; (ii) PD1-H2O: two site polarizabilities are assigned to the hydrogen atoms and one polarizable site "M" is located at a distance d = 0.22Å from the oxygen along the permanent dipole moment vector; (iii) PD2-H2O: two site polarizabilities are assigned to the hydrogen atoms and one polarizable site "M" is located at a distance d = 0.0606Åfrom the oxygen along the permanent dipole moment vector.
Independently on the model, all sites bear the same charge; more details are given in table 1.
A detailed description of polarization methods was given in Refs. [21,23]. Here the main features of the PD method are outlined: the total electric field ( E) acting on each polarizable site is obtained by the field E q produced by the external partial charges and by the field E µ due to the intramolecular and external induced dipole moments: where T q i and T µ ij denote respectively the charge and the dipole field tensor, which for a charge and a point dipole located at the origin can be written as: The induced point dipole on site i is obtained from the total field according to and is computed iteratively until a given threshold of convergence for the induced dipole is reached. It has been shown that for a realistic description of the short range interactions the electric field created by a fixed charge and/or that created by a point dipole should be damped by functions f q (r) and f µ (r); eqs. 5 and 6 are substituted by the following: T where, in the limit of point charges and/or point dipoles we recover the usual expressions, being f q (r) = f µ (r) = 1. On the other hand, when charges and dipoles are thought to be spatially extended, the form of the damping depends on the charge distribution assumed; in ref. [23] the following functional forms were proposed: (i) linear damping: (ii) exponential damping: (iii) gaussian damping: In any of the above possibilites a and b play the role of the characteristic distance of maximum approach. Their value depends on the functional form chosen, and is dependent on each different pair.
Since the linear damping shows a pathological behaviour at r = a or r = b which might be problematic in Molecular Dynamics simulations, in this paper only the exponential and gaussian damping functions are explored in combination with the polarizable models introduced above. All the parameters used are shown in table 1 where, to the water models introduced we add the prefix EXP or GAU depending on the kind of damping function used. The damping parameters for lithium-water dimers are reported in ref. [23]. In the case of chloride-water dimers, the a and b parameters of eqs. 12, 13, 14 and 15 were optimized using the simplex algorithm [45]; the root mean square difference (RMSD) between ab initio and classical dipole moments for all configurations sampled was minimized. In this sense, those parameters represent the best choice which reproduce the ab initio results.

Dimers
Lithium and chloride ions were chosen because of their electronic characteristics: Li + has a high charge density and a low polarizability (α = 0.02875Å 3 ) [46], while Cl − has a lower charge density (being its Van der Waals radius bigger than Lithium) and a high polarizability (α Cl = 5.4822Å 3 ) [47]. The two ions feature opposite qualities which are useful for testing the classical models against ab initio calculations. The total dipole moment of the dimers is sketched in figures 2 and 3. The configurations showed (respectively for the lithium-water and chloride-water dimers), being energetically favoured both at gas and at condensed phase, are the most representative and they weight more than others for condensed phase simulations. It can be appreciated that the models without damping overestimate the ab initio curve at short distances. As explained above, this is due to the so called "Pauli effects" in the multi-body interactions [33]. The only way to reproduce the ab initio curve would be to vary the site polarizabilities both for water and for the ion. In MD simulations of halide-water at condensed phase, this is indeed the typical way to circumvent the "polarization catastrophe", i.e. the rise to infinity of the induced dipoles [50,48,49,14]. Nevertheless, as shown in refs. [24,25] the recourse to a lower polarizability causes the total dipole moment of the dimer to be underestimated at intermediate distances. Moreover, by changing site polarizabilities to hinder the model divergence, force field transferability would be affected, since it can not be assured that in different conditions it works out well.
The development of PDX-H2O models was done under the constraint that the gas phase molecular (ionic) polarizability should be reproduced; some model is even able to reproduce the water polarizability tensor. Accounting for short range electronic repulsion does not require to tune site polarizabilities; as can be seen in figures 2 and 3, in fact, the performance of all models is greatly increased by implementing the simple damping functions introduced in the previous section. In these cases the dipole moment at long and intermediate distances is still preserved, and the short range part is reproduced with small errors. Some differences between all models can be appreciated at short distance. For lithium-water (figure 2) it can be seen that the ab initio curve is fitted very well by all models in both configurations, and a better performance is obtained by the XXX-PD2-H2O models. Both the exponential and the gaussian dampings give comparable results.
On the other hand, for chloride-water, the performance of all models depends on the configuration considered. In the back configuration (panel a) the EXP-PD1-H2O and the GAU-PDM-H2O (not shown) models show the smallest displacement with respect to the other models; while in the cis configuration, the GAU-PD2-H2O and the EXP-PD1-H2O models are the best.
For both ions, though, it should be noticed that, overall, the differences among all models are not big and that they occur at (so short) distances where, as observed in section 2.1, it is very unlikely to access in MD simulations at standard conditions. Moreover those small differences cancel when other configurations (in addition to the most probable ones shown in figures 2 and 3) are considered. It can be concluded that all of the models proposed are good enough to reproduce the ab initio dipole moment of the dimers, for both ion-water dimers.

Condensed Phase
For condensed phase simulations it would not make sense to compare the ab initio total dipole moment of the system with the classical ones. In fact the contributions of many different molecular orientations would involve big cancelations and there would exist many parameter sets which faithfully reproduce the ab initio results. Moreover no insight could be obtained on the contributions of every single species in solution. By calculating the Wannier Centers, ionic and molecular dipole moments are easily computed (see section 2.1). Figure 4 shows the time evolution of the dipole moment as computed from the GAU-PD2-H2O classical model (left panel) and from the Wannier Centers (right panel) for the ion and for one water molecule in the chloride-water system. It can be noticed that for both species the classical and ab initio results are almost superposable.  Further insight in the goodness of each model could be obtained from dipole-dipole correlation plots: in the case of first shell molecules in the lithium-water system, the comparison of PDX-H2O models with different damping functions is shown in figure 5. The PD2-H2O model perform much better than PD1-H2O and PDM-H2O (which, being almost superposable, can not be shown in the same figure); underestimation of ab initio results is common to other classical models recently published where the shift in the classical induced dipole moments rises up to 50% [51,52].
The same conclusion is obtained when the dipole moment distribution functions P (µ) are considered (see figure 6). Either the EXP-PD2-H2O and the GAU-PD2-H2O models give good results both for first shell and bulk water molecules. The corresponding curves for Li + are peaked at 0.01 Debye both for ab initio and classical models. From the analysis of P (µ) it can be observed a small displacement of the PD2-H2O models with respect to ab initio calculations for first shell molecules which escaped from the study of the dipole-dipole correlation plot. Such small shift (around 0.3 Debye, less than 10%) could be due to the fact that water-water interactions are not damped. Li + solvation shell has in fact a tetrahedral geometry where water molecules occupy each vertex at an average distance of ≈ 2Å and are rigidly oriented towards the ions [53]. The fact that first shell molecules' dipole moment is so highly influenced highlights that no-polarizable models would miss important electrostatic contributions and that accounting for polarizable interactions is needed also for ions with low (almost zero) polarizability. For what concerns bulk water molecules, it seems that they are not polarized by the ion, since their average dipole moment agrees with what was found in previous calculations [54] of pure liquid water.
A similar (but smaller) shift in the point dipole distribution functions of first shell molecules is found for chloride-water (see figure 7). Also in this case both exponential and gaussian damping functions are almost equally good in reproducing water induced dipole moments.
On the other hand, if we look at the induced dipole moment on the chloride ion, it can be observed that the exponential damping fails, since it gives a much broader distribution peaked    at higher values of µ; only by implementing the gaussian damping it is possible to reproduce with high precision the ab initio distriution. This is in contrast to what was found for the dimer, where both kind of damping functions could be used to model the Cl − -water system.

Concluding Remarks
The use of intermolecular damping functions for modeling ion-water polarizable interaction is of great importance for accounting for short range electron cloud repulsion. This is evident already for the dimer at gas phase; in this case the PD2-H2O water model with the two types of damping functions studied, showed to be slightly better then the PDX-H2O models explored in this and previous contributions. For dimers all models are interchangeable, since all of them fit the ab initio curve with small errors; condensed phase systems are characterized by manybody interactions which account for the different performance of all models: for the systems studied, only the PD2-H2O model is able to catch all features of ion and water polarization, giving good results both for lithium-water and chloride-water systems. In fact, PDM-H2O and PD1-H2O highly underestimate water polarization such that in the dipole-dipole correlation plot all points are below the bisetrix. Such feature could not be observed in the calculations at gas phase because only the total dipole moment, rather than the molecular (ionic) one was computed; in that case it is highly probable that the ionic dipole moment overestimation could cancel out with the error in the water dipole moment. The good performance of PD2-H2O is probably due to its ability to describe the anisotropic polarizability of the water molecule [32].
For what concerns the parameterization of the damping functions, it has been shown that only the GAU-PD2-H2O model is able to reproduce all ab initio results both for ion and water polarization. This model is the only candidate for further force field development given its high portability from gas to condensed phase for the ions considered here.