The Madrid-2019 force field for electrolytes in water using TIP4P/2005 and scaled charges: Extension to the ions F-, Br-, I-, Rb+, and Cs.

In this work, an extension of the Madrid-2019 force field is presented. We have added the cations Rb+ and Cs+ and the anions F-, Br-, and I-. These ions were the remaining alkaline and halogen ions, not previously considered in the Madrid-2019 force field. The force field, denoted as Madrid-2019-Extended, does not include polarizability and uses the TIP4P/2005 model of water and scaled charges for the ions. A charge of ±0.85e is assigned to monovalent ions. The force field developed provides an accurate description of aqueous solution densities over a wide range of concentrations up to the solubility limit of each salt studied. Good predictions of viscosity and diffusion coefficients are obtained for concentrations below 2 m. Structural properties obtained with this force field are also in reasonable agreement with the experiment. The number of contact ion pairs has been controlled to be low so as to avoid precipitation of the system at concentrations close to the experimental solubility limit. A comprehensive comparison of the performance for aqueous solutions of alkaline halides of force fields of electrolytes using scaled and integer charges is now possible. This comparison will help in the future to learn about the benefits and limitations of the use of scaled charges to describe electrolyte solutions.

A force field of Li+, Na+, K+, Mg2+, Ca2+, Cl-, and SO4 2- in aqueous solution based on the TIP4P/2005 water model and scaled charges for the ions.

In this work, a force field for several ions in water is proposed. In particular, we consider the cations Li+, Na+, K+, Mg2+, and Ca2+ and the anions Cl- and SO4 2-. These ions were selected as they appear in the composition of seawater, and they are also found in biological systems. The force field proposed (denoted as Madrid-2019) is nonpolarizable, and both water molecules and sulfate anions are rigid. For water, we use the TIP4P/2005 model. The main idea behind this work is to further explore the possibility of using scaled charges for describing ionic solutions. Monovalent and divalent ions are modeled using charges of 0.85 and 1.7, respectively (in electron units). The model allows a very accurate description of the densities of the solutions up to high concentrations. It also gives good predictions of viscosities up to 3 m concentrations. Calculated structural properties are also in reasonable agreement with the experiment. We have checked that no crystallization occurred in the simulations at concentrations similar to the solubility limit. A test for ternary mixtures shows that the force field provides excellent performance at an affordable computer cost. In summary, the use of scaled charges, which could be regarded as an effective and simple way of accounting for polarization (at least to a certain extend), improves the overall description of ionic systems in water. However, for purely ionic systems, scaled charges will not adequately describe neither the solid nor the melt.

A potential model for sodium chloride solutions based on the TIP4P/2005 water model.

Despite considerable efforts over more than two decades, our knowledge of the interactions in electrolyte solutions is not yet satisfactory. Not even one of the most simple and important aqueous solutions, NaCl(aq), escapes this assertion. A requisite for the development of a force field for any water solution is the availability of a good model for water. Despite the fact that TIP4P/2005 seems to fulfill the requirement, little work has been devoted to build a force field based on TIP4P/2005. In this work, we try to fill this gap for NaCl(aq). After unsuccessful attempts to produce accurate predictions for a wide range of properties using unity ionic charges, we decided to follow recent suggestions indicating that the charges should be scaled in the ionic solution. In this way, we have been able to develop a satisfactory non-polarizable force field for NaCl(aq). We evaluate a number of thermodynamic properties of the solution (equation of state, maximum in density, enthalpies of solution, activity coefficients, radial distribution functions, solubility, surface tension, diffusion coefficients, and viscosity). Overall the results for the solution are very good. An important achievement of our model is that it also accounts for the dynamical properties of the solution, a test for which the force fields so far proposed failed. The same is true for the solubility and for the maximum in density where the model describes the experimental results almost quantitatively. The price to pay is that the model is not so good at describing NaCl in the solid phase, although the results for several properties (density and melting temperature) are still acceptable. We conclude that the scaling of the charges improves the overall description of NaCl aqueous solutions when the polarization is not included.

Nucleation free-energy barriers with Hybrid Monte-Carlo/Umbrella Sampling.

The aim of this work is to evaluate nucleation free-energy barriers using molecular dynamics (MD). More specifically, we use a combination of Hybrid Monte Carlo (HMC) and an Umbrella Sampling scheme, and compute the crystallisation barrier of NaCl from its melt. Firstly the convergence and performance of HMC for different time-steps and the number of MD steps within a HMC cycle are assessed. The calculated potential energies and densities converge regardless of the chosen time-step. However the acceptance ratio of the Metropolis step within the HMC scheme strongly depends on the time-step and affects the performance. It is shown that the acceptance ratio is close to 100% for time-steps of the order of those commonly used in molecular dynamics runs. We then explore the results obtained with a "non-Metropolised" version of HMC where the MD trajectories are always accepted (omitting the Metropolis criteria) and conclude that they are satisfactory for time-steps below 5 fs. Next, HMC is combined with Umbrella Sampling (HMC/US) to compute the nucleation free-energy for both the standard and the "non-Metropolised" HMC (using a small time-step) and in both cases find excellent agreement with the reported values. To conclude, we explore approximations to the HMC/US technique implementing HMC with isothermal-isobaric MD trajectories. The computed nucleation free-energy curve is coincident, within the statistical error, with previous calculations.

Determining the phase diagram of water from direct coexistence simulations: the phase diagram of the TIP4P/2005 model revisited.

Direct coexistence simulations between the fluid and solid phases are performed for several ices. For ices Ih and VII it has already been shown that the methodology is successful and the melting point is in agreement with that obtained from free energy calculations. In this work the methodology is applied to ices II, III, V, and VI. The lengths of the direct coexistence runs for the high pressure polymorphs are not too long and last less than 20 ns for all ices except for ice II where longer runs (of about 150 ns) are needed. For ices II, V, and VI the results obtained are completely consistent with those obtained from free energy calculations. However, for ice III it is found that the melting point from direct coexistence simulations is higher than that obtained from free energy calculations, the difference being greater than the statistical error. Since ice III presents partial proton orientational disorder, the departure is attributed to differences in the partial proton order in the water model with respect to that found in the experiment. The phase diagram of the TIP4P/2005 model is recalculated using the melting points obtained from direct coexistence simulations. The new phase diagram is similar to the previous one except for the coexistence lines where ice III is involved. The range of stability of ice III on the p-T plot of the phase diagram increases significantly. It is seen that the model qualitatively describes the phase diagram of water. In this work it is shown that the complete phase diagram of water including ices Ih, II, III, V, VI, VII, and the fluid phase can be obtained from direct coexistence simulations without the need of free energy calculations.

Homogeneous ice nucleation at moderate supercooling from molecular simulation.

Among all of the freezing transitions, that of water into ice is probably the most relevant to biology, physics, geology, or atmospheric science. In this work, we investigate homogeneous ice nucleation by means of computer simulations. We evaluate the size of the critical cluster and the nucleation rate for temperatures ranging between 15 and 35 K below melting. We use the TIP4P/2005 and the TIP4P/ice water models. Both give similar results when compared at the same temperature difference with the model's melting temperature. The size of the critical cluster varies from ∼8000 molecules (radius = 4 nm) at 15 K below melting to ∼600 molecules (radius = 1.7 nm) at 35 K below melting. We use Classical Nucleation Theory (CNT) to estimate the ice-water interfacial free energy and the nucleation free-energy barrier. We obtain an interfacial free energy of 29(3) mN/m from an extrapolation of our results to the melting temperature. This value is in good agreement both with experimental measurements and with previous estimates from computer simulations of TIP4P-like models. Moreover, we obtain estimates of the nucleation rate from simulations of the critical cluster at the barrier top. The values we get for both models agree within statistical error with experimental measurements. At temperatures higher than 20 K below melting, we get nucleation rates slower than the appearance of a critical cluster in all water of the hydrosphere during the age of the universe. Therefore, our simulations predict that water freezing above this temperature must necessarily be heterogeneous.

Note: Equation of state and compressibility of supercooled water: simulations and experiment.

The equation of state and the isothermal compressibility of supercooled water for pressures up to 3000 bar obtained from computer simulations of the TIP4P/2005 model are compared to recent experimental results. The agreement between the simulations and experimental results is quite good. This reinforces the idea that the Widom line and the liquid-liquid phase separation found in the simulations should also exist in real water.

Heat capacity of water: A signature of nuclear quantum effects.

In this note we present results for the heat capacity at constant pressure for the TIP4PQ/2005 model, as obtained from path-integral simulations. The model does a rather good job of describing both the heat capacity of ice I(h) and of liquid water. Classical simulations using the TIP4P/2005, TIP3P, TIP4P, TIP4P-Ew, simple point charge/extended, and TIP5P models are unable to reproduce the heat capacity of water. Given that classical simulations do not satisfy the third law of thermodynamics, one would expect such a failure at low temperatures. However, it seems that for water, nuclear quantum effects influence the heat capacities all the way up to room temperature. The failure of classical simulations to reproduce C(p) points to the necessity of incorporating nuclear quantum effects to describe this property accurately.

What ice can teach us about water interactions: a critical comparison of the performance of different water models.

The performance of several popular water models (TIP3P, TIP4P, TIP5P and TIP4P/2005) is analyzed. For that purpose the predictions for ten different properties of water are investigated, namely: 1. vapour-liquid equilibria (VLE) and critical temperature; 2. surface tension; 3. densities of the different solid structures of water (ices); 4. phase diagram; 5. melting-point properties; 6. maximum in the density of water at room pressure and thermal coefficients alpha and KT; 7. structure of liquid water and ice; 8. equation of state at high pressures; 9. self-diffusion coefficient; 10. dielectric constant. For each property, the performance of each model is analyzed in detail with a critical discussion of the possible reason of the success or failure of the model. A final judgement on the quality of these models is provided. TIP4P/2005 provides the best description of almost all properties of the list, the only exception being the dielectric constant. In second position, TIP5P and TIP4P yield a similar performance overall, and the last place with the poorest description of the water properties is provided by TIP3P. The ideas leading to the proposal and design of the TIP4P/2005 are also discussed in detail. TIP4P/2005 is probably close to the best description of water that can be achieved with a non-polarizable model described by a single Lennard-Jones (LJ) site and three charges.

Properties of ices at 0 K: a test of water models.

The properties of ices Ih, II, III, V, and VI at zero temperature and pressure are determined by computer simulation for several rigid water models (SPC/E, TIP5P, TIP4P/Ice, and TIP4P/2005). The energies of the different ices at zero temperature and pressure (relative to the ice II energy) are compared to the experimental results of Whalley [J. Chem. Phys. 81, 4087 (1984)]. TIP4P/Ice and TIP4P/2005 provide a qualitatively correct description of the relative energies of the ices at these conditions. In fact, only these two models provide the correct ordering in energies. For the SPC/E and TIP5P models, ice II is the most stable phase at zero temperature and pressure whereas for TIP4P/Ice and TIP4P/2005 ice Ih is the most stable polymorph. These results are in agreement with the relative stabilities found at higher temperatures. The solid-solid phase transitions at 0 K are determined. The predicted pressures are in good agreement with those obtained from free energy calculations.

Computer simulation of the ionic atmosphere around Z-DNA.

We describe a coarse-grained model for Z-DNA that mimics the DNA shape with a relatively small number of repulsive interaction sites. In addition, negative charges are placed at the phosphate positions. The ionic atmosphere around this grooved Z-DNA model is then investigated with Monte Carlo simulation. Cylindrically averaged concentration profiles as well as the spatial distribution of ions have been calculated. The results are compared to those for other DNA models differing in the repulsive core. This allows the examination of the effect of the DNA shape in the ionic distribution. It is seen that the penetrability of the ions to the DNA groove plays an important role in the ionic distribution. The results are also compared with those reported for B-DNA. In both conformers the ions are structured in alternating layers of positive and negative charge. In Z-DNA the layers are more or less concentric to the molecular axis. Besides, no coions enter into the single groove of this conformer. On the contrary, the alternating layers of B-DNA are also structured along the axial coordinate with some coions penetrating into the major groove. In both cases we have found five preferred locations of the counterions and two for the coions. The concentration of counterions reaches its absolute maximum at the narrow Z-DNA groove and at the minor groove of B-DNA, the value of the maximum being higher in the Z conformer.

Vapor-liquid equilibria from the triple point up to the critical point for the new generation of TIP4P-like models: TIP4P/Ew, TIP4P/2005, and TIP4P/ice.

The vapor-liquid equilibria of three recently proposed water models have been computed using Gibbs-Duhem simulations. These models are TIP4P/Ew, TIP4P/2005, and TIP4P/ice and can be considered as modified versions of the TIP4P model. By design TIP4P reproduces the vaporization enthalpy of water at room temperature, whereas TIP4P/Ew and TIP4P/2005 match the temperature of maximum density and TIP4P/ice the melting temperature of water. Recently, the melting point for each of these models has been computed, making it possible for the first time to compute the complete vapor-liquid equilibria curve from the triple point to the critical point. From the coexistence results at high temperature, it is possible to estimate the critical properties of these models. None of them is capable of reproducing accurately the critical pressure or the vapor pressures and densities. Additionally, in the cases of TIP4P and TIP4P/ice the critical temperatures are too low and too high, respectively, compared to the experimental value. However, models accounting for the density maximum of water, such as TIP4P/Ew and TIP4P/2005 provide a better estimate of the critical temperature. In particular, TIP4P/2005 provides a critical temperature just 7 K below the experimental result as well as an extraordinarily good description of the liquid densities from the triple point to the critical point. All TIP4P-like models present a ratio of the triple point temperature to the critical point temperature of about 0.39, compared with the experimental value of 0.42. As is the case for any effective potential neglecting many body forces, TIP4P/2005 fails in describing simultaneously the vapor and the liquid phases of water. However, it can be considered as one of the best effective potentials of water for describing condensed phases, both liquid and solid. In fact, it provides a completely coherent view of the phase diagram of water including fluid-solid, solid-solid, and vapor-liquid equilibria.

Relation between the melting temperature and the temperature of maximum density for the most common models of water.

Water exhibits a maximum in density at normal pressure at 4 degrees above its melting point. The reproduction of this maximum is a stringent test for potential models used commonly in simulations of water. The relation between the melting temperature and the temperature of maximum density for these potential models is unknown mainly due to our ignorance about the melting temperature of these models. Recently we have determined the melting temperature of ice I(h) for several commonly used models of water (SPC, SPC/E, TIP3P, TIP4P, TIP4P/Ew, and TIP5P). In this work we locate the temperature of maximum density for these models. In this way the relative location of the temperature of maximum density with respect to the melting temperature is established. For SPC, SPC/E, TIP3P, TIP4P, and TIP4P/Ew the maximum in density occurs at about 21-37 K above the melting temperature. In all these models the negative charge is located either on the oxygen itself or on a point along the H-O-H bisector. For the TIP5P and TIP5P-E models the maximum in density occurs at about 11 K above the melting temperature. The location of the negative charge appears as a geometrical crucial factor to the relative position of the temperature of maximum density with respect to the melting temperature.

A potential model for the study of ices and amorphous water: TIP4P/Ice.

The ability of several water models to predict the properties of ices is discussed. The emphasis is put on the results for the densities and the coexistence curves between the different ice forms. It is concluded that none of the most commonly used rigid models is satisfactory. A new model specifically designed to cope with solid-phase properties is proposed. The parameters have been obtained by fitting the equation of state and selected points of the melting lines and of the coexistence lines involving different ice forms. The phase diagram is then calculated for the new potential. The predicted melting temperature of hexagonal ice (Ih) at 1 bar is 272.2 K. This excellent value does not imply a deterioration of the rest of the properties. In fact, the predictions for both the densities and the coexistence curves are better than for TIP4P, which previously yielded the best estimations of the ice properties.

The melting temperature of the most common models of water.

The melting temperature of ice I(h) for several commonly used models of water (SPC, SPC/E,TIP3P,TIP4P, TIP4P/Ew, and TIP5P) is obtained from computer simulations at p = 1 bar. Since the melting temperature of ice I(h) for the TIP4P model is now known [E. Sanz, C. Vega, J. L. F. Abascal, and L. G. MacDowell, Phys. Rev. Lett. 92, 255701 (2004)], it is possible to use the Gibbs-Duhem methodology [D. Kofke, J. Chem. Phys. 98, 4149 (1993)] to evaluate the melting temperature of ice I(h) for other potential models of water. We have found that the melting temperatures of ice I(h) for SPC, SPC/E, TIP3P, TIP4P, TIP4P/Ew, and TIP5P models are T = 190 K, 215 K, 146 K, 232 K, 245 K, and 274 K, respectively. The relative stability of ice I(h) with respect to ice II for these models has also been considered. It turns out that for SPC, SPC/E, TIP3P, and TIP5P the stable phase at the normal melting point is ice II (so that ice I(h) is not a thermodynamically stable phase for these models). For TIP4P and TIP4P/Ew, ice I(h) is the stable solid phase at the standard melting point. The location of the negative charge along the H-O-H bisector appears as a critical factor in the determination of the relative stability between the I(h) and II ice forms. The methodology proposed in this paper can be used to investigate the effect upon a coexistence line due to a change in the potential parameters.

A general purpose model for the condensed phases of water: TIP4P/2005.

A potential model intended to be a general purpose model for the condensed phases of water is presented. TIP4P/2005 is a rigid four site model which consists of three fixed point charges and one Lennard-Jones center. The parametrization has been based on a fit of the temperature of maximum density (indirectly estimated from the melting point of hexagonal ice), the stability of several ice polymorphs and other commonly used target quantities. The calculated properties include a variety of thermodynamic properties of the liquid and solid phases, the phase diagram involving condensed phases, properties at melting and vaporization, dielectric constant, pair distribution function, and self-diffusion coefficient. These properties cover a temperature range from 123 to 573 K and pressures up to 40,000 bar. The model gives an impressive performance for this variety of properties and thermodynamic conditions. For example, it gives excellent predictions for the densities at 1 bar with a maximum density at 278 K and an averaged difference with experiment of 7 x 10(-4) g/cm3.

Phase diagram of water from computer simulation.

The phase diagram of water as obtained from computer simulations is presented for the first time for two of the most popular models of water, TIP4P and SPC/E. This Letter shows that the prediction of the phase diagram is an extremely stringent test for any water potential function, and that it may be useful in developing improved potentials. The TIP4P model provides a qualitatively correct description of the phase diagram, unlike the SPC/E model which fails in this purpose. New behavior not yet observed experimentally is predicted by the simulations: the existence of metastable reentrant behavior in the melting curves of the low density ices (I,III,V) such that it could be possible to transform them into amorphous phases by adequate changes in pressure.

Tracing the phase diagram of the four-site water potential (TIP4P).

We present here the phase diagram for one of the most popular water models, the four-point transferable intermolecular potential (TIP4P) model. We show that TIP4P model, does indeed provide a qualitatively correct description of the phase diagram of water. The melting line of the five-point transferable intermolecular potential (TIP5P) at low pressures is also presented.

Order-disorder transition in the solid phase of a charged hard sphere model.

We investigate the solid phases of the restricted primitive model (RPM). Monte Carlo simulations show the existence of an order-disorder transition from a substitutionally disordered face centered cubic lattice (fcc) to a new ordered fcc structure which is proposed as the ground state of the RPM at the close packing density. Our results suggest that the new phase might turn out in a new triple point in the RPM phase diagram involving three solid phases: CsCl, fcc ordered and fcc disordered structures. The order-disorder transition is also studied using the cell theory. The theory shows good agreement with the simulation results and suggests that the transition is weakly first order.