Computation of Electrical Conductivities of Aqueous Electrolyte Solutions: Two Surfaces, One Property.

In this work, we computed electrical conductivities under ambient conditions of aqueous NaCl and KCl solutions by using the Einstein-Helfand equation. Common force fields (charge q = ±1 e) do not reproduce the experimental values of electrical conductivities, viscosities, and diffusion coefficients. Recently, we proposed the idea of using different charges to describe the potential energy surface (PES) and the dipole moment surface (DMS). In this work, we implement this concept. The equilibrium trajectories required to evaluate electrical conductivities (within linear response theory) were obtained by using scaled charges (with the value q = ±0.75 e) to describe the PES. The potential parameters were those of the Madrid-Transport force field, which accurately describe viscosities and diffusion coefficients of these ionic solutions. However, integer charges were used to compute the conductivities (thus describing the DMS). The basic idea is that although the scaled charge describes the ion-water interaction better, the integer charge reflects the value of the charge that is transported due to the electric field. The agreement obtained with experiments is excellent, as for the first time electrical conductivities (and the other transport properties) of NaCl and KCl electrolyte solutions are described with high accuracy for the whole concentration range up to their solubility limit. Finally, we propose an easy way to obtain a rough estimate of the actual electrical conductivity of the potential model under consideration using the approximate Nernst-Einstein equation, which neglects correlations between different ions.

"In Silico" Seawater.

Many important processes affecting the earth's climate are determined by the physical properties of seawater. In addition, desalination of seawater is a significant source of drinking water for the human population living in coastal areas. Since the physical properties of seawater governing these processes depend on the molecular interactions among its components, a deeper knowledge of seawater at the molecular level would contribute to a better understanding of these phenomena. However, in strong contrast with the situation in other areas such as biomolecules or materials science, molecular simulation studies reporting the physical properties of seawater are currently lacking. This is probably due to the usual perception of the seawater composition being too complex to approach. This point of view ignores the fact that physical properties of seawater are dependent on a single parameter representing the composition, namely the salinity. This is because the relative proportions of any two major constituents of seasalt are always the same. Another obstacle to performing molecular simulations of seawater could have been the unavailability of a satisfactory force field representing the interactions between water molecules and dissolved substances. However, this drawback has recently been overcome with the proposal of the Madrid-2019 force field. In this work we show for the first time that molecular simulation of seawater is feasible. We have performed molecular dynamics simulations of a system, the composition of which is close to the average composition of standard seawater and with the molecular interactions given by the Madrid-2019 force field. In this way we are able to provide quantitative or semiquantitative predictions for a number of relevant physical properties of seawater for temperatures and salinities from the oceanographic range to those relevant to desalination processes. The computed magnitudes include static (density), dynamical (viscosity and diffusion coefficients), structural (ionic hydration, ion-ion distribution functions), and interfacial (surface tension) properties.

Two-structure thermodynamics for the TIP4P/2005 model of water covering supercooled and deeply stretched regions.

One of the most promising frameworks for understanding the anomalies of cold and supercooled water postulates the existence of two competing, interconvertible local structures. If the non-ideality in the Gibbs energy of mixing overcomes the ideal entropy of mixing of these two structures, a liquid-liquid phase transition, terminated at a liquid-liquid critical point, is predicted. Various versions of the "two-structure equation of state" (TSEOS) based on this concept have shown remarkable agreement with both experimental data for metastable, deeply supercooled water and simulations of molecular water models. However, existing TSEOSs were not designed to describe the negative pressure region and do not account for the stability limit of the liquid state with respect to the vapor. While experimental data on supercooled water at negative pressures may shed additional light on the source of the anomalies of water, such data are very limited. To fill this gap, we have analyzed simulation results for TIP4P/2005, one of the most accurate classical water models available. We have used recently published simulation data, and performed additional simulations, over a broad range of positive and negative pressures, from ambient temperature to deeply supercooled conditions. We show that, by explicitly incorporating the liquid-vapor spinodal into a TSEOS, we are able to match the simulation data for TIP4P/2005 with remarkable accuracy. In particular, this equation of state quantitatively reproduces the lines of extrema in density, isothermal compressibility, and isobaric heat capacity. Contrary to an explanation of the thermodynamic anomalies of water based on a "retracing spinodal," the liquid-vapor spinodal in the present TSEOS continues monotonically to lower pressures upon cooling, influencing but not giving rise to density extrema and other thermodynamic anomalies.

Molecular mechanism for cavitation in water under tension.

Despite its relevance in biology and engineering, the molecular mechanism driving cavitation in water remains unknown. Using computer simulations, we investigate the structure and dynamics of vapor bubbles emerging from metastable water at negative pressures. We find that in the early stages of cavitation, bubbles are irregularly shaped and become more spherical as they grow. Nevertheless, the free energy of bubble formation can be perfectly reproduced in the framework of classical nucleation theory (CNT) if the curvature dependence of the surface tension is taken into account. Comparison of the observed bubble dynamics to the predictions of the macroscopic Rayleigh-Plesset (RP) equation, augmented with thermal fluctuations, demonstrates that the growth of nanoscale bubbles is governed by viscous forces. Combining the dynamical prefactor determined from the RP equation with CNT based on the Kramers formalism yields an analytical expression for the cavitation rate that reproduces the simulation results very well over a wide range of pressures. Furthermore, our theoretical predictions are in excellent agreement with cavitation rates obtained from inclusion experiments. This suggests that homogeneous nucleation is observed in inclusions, whereas only heterogeneous nucleation on impurities or defects occurs in other experiments.

A comprehensive scenario of the thermodynamic anomalies of water using the TIP4P/2005 model.

The striking behavior of water has deserved it to be referred to as an "anomalous" liquid. The water anomalies are greatly amplified in metastable (supercooled and/or stretched) regions. This makes difficult a complete experimental description since, beyond certain limits, the metastable phase necessarily transforms into the stable one. Theoretical interpretation of the water anomalies could then be based on simulation results of well validated water models. But the analysis of the simulations has not yet reached a consensus. In particular, one of the most popular theoretical scenarios-involving the existence of a liquid-liquid critical point (LLCP)-is disputed by several authors. In this work, we propose to use a number of exact thermodynamic relations which may shed light on this issue. Interestingly, these relations may be tested in a region of the phase diagram which is outside the LLCP thus avoiding the problems associated to the coexistence region. The central property connected to other water anomalies is the locus of temperatures at which the density along isobars attain a maximum (TMD line) or a minimum (TmD). We have performed computer simulations to evaluate the TMD and TmD for a successful water model, namely, TIP4P/2005. We have also evaluated the vapor-liquid (VL) spinodal in the region of large negative pressures. The shape of these curves and their connection to the extrema of some response functions, in particular the isothermal compressibility and heat capacity at constant pressure, provides very useful information which may help to elucidate the validity of the theoretical proposals. In this way, we are able to present for the first time a comprehensive scenario of the thermodynamic water anomalies for TIP4P/2005 and their relation to the vapor-liquid spinodal. The overall picture shows a remarkable similarity with the corresponding one for the ST2 water model, for which the existence of a LLCP has been demonstrated in recent years. It also provides a hint as to where the long-sought for extrema in response functions might become accessible to experiments.

Bubble nucleation in simple and molecular liquids via the largest spherical cavity method.

In this work, we propose a methodology to compute bubble nucleation free energy barriers using trajectories generated via molecular dynamics simulations. We follow the bubble nucleation process by means of a local order parameter, defined by the volume of the largest spherical cavity (LSC) formed in the nucleating trajectories. This order parameter simplifies considerably the monitoring of the nucleation events, as compared with the previous approaches which require ad hoc criteria to classify the atoms and molecules as liquid or vapor. The combination of the LSC and the mean first passage time technique can then be used to obtain the free energy curves. Upon computation of the cavity distribution function the nucleation rate and free-energy barrier can then be computed. We test our method against recent computations of bubble nucleation in simple liquids and water at negative pressures. We obtain free-energy barriers in good agreement with the previous works. The LSC method provides a versatile and computationally efficient route to estimate the volume of critical bubbles the nucleation rate and to compute bubble nucleation free-energies in both simple and molecular liquids.

Detecting vapour bubbles in simulations of metastable water.

The investigation of cavitation in metastable liquids with molecular simulations requires an appropriate definition of the volume of the vapour bubble forming within the metastable liquid phase. Commonly used approaches for bubble detection exhibit two significant flaws: first, when applied to water they often identify the voids within the hydrogen bond network as bubbles thus masking the signature of emerging bubbles and, second, they lack thermodynamic consistency. Here, we present two grid-based methods, the M-method and the V-method, to detect bubbles in metastable water specifically designed to address these shortcomings. The M-method incorporates information about neighbouring grid cells to distinguish between liquid- and vapour-like cells, which allows for a very sensitive detection of small bubbles and high spatial resolution of the detected bubbles. The V-method is calibrated such that its estimates for the bubble volume correspond to the average change in system volume and are thus thermodynamically consistent. Both methods are computationally inexpensive such that they can be used in molecular dynamics and Monte Carlo simulations of cavitation. We illustrate them by computing the free energy barrier and the size of the critical bubble for cavitation in water at negative pressure.

Anomalies in bulk supercooled water at negative pressure.

Water anomalies still defy explanation. In the supercooled liquid, many quantities, for example heat capacity and isothermal compressibility κT, show a large increase. The question arises if these quantities diverge, or if they go through a maximum. The answer is key to our understanding of water anomalies. However, it has remained elusive in experiments because crystallization always occurred before any extremum is reached. Here we report measurements of the sound velocity of water in a scarcely explored region of the phase diagram, where water is both supercooled and at negative pressure. We find several anomalies: maxima in the adiabatic compressibility and nonmonotonic density dependence of the sound velocity, in contrast with a standard extrapolation of the equation of state. This is reminiscent of the behavior of supercritical fluids. To support this interpretation, we have performed simulations with the 2005 revision of the transferable interaction potential with four points. Simulations and experiments are in near-quantitative agreement, suggesting the existence of a line of maxima in κT (LMκT). This LMκT could either be the thermodynamic consequence of the line of density maxima of water [Sastry S, Debenedetti PG, Sciortino F, Stanley HE (1996) Phys Rev E 53:6144-6154], or emanate from a critical point terminating a liquid-liquid transition [Sciortino F, Poole PH, Essmann U, Stanley HE (1997) Phys Rev E 55:727-737]. At positive pressure, the LMκT has escaped observation because it lies in the "no man's land" beyond the homogeneous crystallization line. We propose that the LMκT emerges from the no man's land at negative pressure.

Homogeneous bubble nucleation in water at negative pressure: a Voronoi polyhedra analysis.

We investigate vapor bubble nucleation in metastable TIP4P/2005 water at negative pressure via the Mean First Passage Time (MFPT) technique using the volume of the largest bubble as a local order parameter. We identify the bubbles in the system by means of a Voronoi-based analysis of the molecular dynamics trajectories. By comparing the features of the tessellation of liquid water at ambient conditions to those of the same system with an empty cavity we are able to discriminate vapor (or interfacial) molecules from the bulk ones. This information is used to follow the time evolution of the largest bubble until the system cavitates at 280 K above and below the spinodal line. At the pressure above the spinodal line, the MFPT curve shows the expected shape for a moderately metastable liquid from which we estimate the bubble nucleation rate and the size of the critical cluster. The nucleation rate estimated using Classical Nucleation Theory turns out to be about 8 order of magnitude lower than the one we compute by means of MFPT. The behavior at the pressure below the spinodal line, where the liquid is thermodynamically unstable, is remarkably different, the MFPT curve being a monotonous function without any inflection point.

A flexible model for water based on TIP4P/2005.

A new flexible water model, TIP4P/2005f, is developed. The idea was to add intramolecular degrees of freedom to the successful rigid model TIP4P/2005 in order to try to improve the predictions for some properties, and to enable the calculation of new ones. The new model incorporates flexibility by means of a Morse potential for the bond stretching and a harmonic term for the angle bending. The parameters have been fitted to account for the peaks of the infrared spectrum of liquid water and to produce an averaged geometry close to that of TIP4P/2005. As for the intermolecular interactions, only a small change in the σ parameter of the Lennard-Jones potential has been introduced. The overall predictions are very close to those of TIP4P/2005. This ensures that the new model may be used with the same confidence as its predecessor in studies where a flexible model is advisable.

Simulating water with rigid non-polarizable models: a general perspective.

Over the last forty years many computer simulations of water have been performed using rigid non-polarizable models. Since these models describe water interactions in an approximate way it is evident that they cannot reproduce all of the properties of water. By now many properties for these kinds of models have been determined and it seems useful to compile some of these results and provide a critical view of the successes and failures. In this paper a test is proposed in which 17 properties of water, from the vapour and liquid to the solid phases, are taken into account to evaluate the performance of a water model. A certain number of points between zero (bad agreement) and ten (good agreement) are given for the predictions of each model and property. We applied the test to five rigid non-polarizable models, TIP3P, TIP5P, TIP4P, SPC/E and TIP4P/2005, obtaining an average score of 2.7, 3.7, 4.7, 5.1, and 7.2 respectively. Thus although no model reproduces all properties, some models perform better than others. It is clear that there are limitations for rigid non-polarizable models. Neglecting polarizability prevents an accurate description of virial coefficients, vapour pressures, critical pressure and dielectric constant. Neglecting nuclear quantum effects prevents an accurate description of the structure, the properties of water below 120 K and the heat capacity. It is likely that for rigid non-polarizable models it may not be possible to increase the score in the test proposed here beyond 7.6. To get closer to experiment, incorporating polarization and nuclear quantum effects is absolutely required even though a substantial increase in computer time should be expected. The test proposed here, being quantitative and selecting properties from all phases of water can be useful in the future to identify progress in the modelling of water.

Widom line and the liquid-liquid critical point for the TIP4P/2005 water model.

The Widom line and the liquid-liquid critical point of water in the deeply supercooled region are investigated via computer simulation of the TIP4P/2005 model. The Widom line has been calculated as the locus of compressibility maxima. It is quite close to the experimental homogeneous nucleation line and, in the region studied, it is almost parallel to the curve of temperatures of maximum density at fixed pressure. The critical temperature is determined by examining which isotherm has a region with flat slope. An interpolation in the Widom line gives the rest of the critical parameters. The computed critical parameters are T(c)=193 K, p(c)=1350 bar, and ρ(c)=1.012 g/cm(3). Given the performance of the model for the anomalous properties of water and for the properties of ice phases, the calculated critical parameters are probably close to those of real water.

The shear viscosity of rigid water models.

In this work, the shear viscosity at ambient conditions of several water models (SPC/E, TIP4P, TIP5P, and TIP4P/2005) is evaluated using the Green-Kubo formalism. The performance of TIP4P/2005 is excellent, that of SPC/E and TIP5P is more or less acceptable, whereas TIP4P and especially TIP3P give a poor agreement with experiment. Further calculations have been carried out for TIP4P/2005 to provide a wider assessment of its performance. In accordance with experimental data, TIP4P/2005 predicts a minimum in the shear viscosity for the 273 K isotherm, a shift in the minimum toward lower pressures at 298 K, and its disappearance at 373 K.

Triple points and coexistence properties of the dense phases of water calculated using computer simulation.

In recent years our group has performed a thorough study of the dense region of the phase diagram of water. In this paper we report the numerical results for the triple points and for the change in volume, Deltav, and enthalpy, DeltaH, along the coexistence lines involving liquid water and/or ices ih (hexagonal), II, III, V and VI for several simple water models. The predictions for Deltav using the TIP4P/2005 model are in excellent agreement with the experimental values. As to DeltaH for the same model, the computer simulation results are also satisfactory although there are small but significant differences between simulation and experiment.

Dipole-quadrupole force ratios determine the ability of potential models to describe the phase diagram of water.

We show that, in order to yield an acceptable phase diagram for water, the key ingredient is the relative strength of the dipolar and quadrupolar forces. We also demonstrate that for three-point-charge water models the relative strength of the dipolar to quadrupolar forces is determined by the position of the negative charge. This provides a theoretical background--supported by the experiment--to our knowledge of the water interactions which at present is essentially dominated by empirical models. In addition, it can be quite useful in the development of improved water models with a more physical basis.

The melting point of hexagonal ice (Ih) is strongly dependent on the quadrupole of the water models.

It is shown that quadrupolar interactions play a determinant role in the melting temperatures of common water models and that there is a simple relationship between the strength of the quadrupolar forces and the position of the negative charge; our conclusion is that acceptable predictions for the melting temperature can only be obtained when the negatively charged site is shifted from the oxygen atom towards the hydrogens.

The melting point of ice Ih for common water models calculated from direct coexistence of the solid-liquid interface.

In this work we present an implementation for the calculation of the melting point of ice I(h) from direct coexistence of the solid-liquid interface. We use molecular dynamics simulations of boxes containing liquid water and ice in contact. The implementation is based on the analysis of the evolution of the total energy along NpT simulations at different temperatures. We report the calculation of the melting point of ice I(h) at 1 bar for seven water models: SPC/E, TIP4P, TIP4P-Ew, TIP4P/ice, TIP4P/2005, TIP5P, and TIP5P-E. The results for the melting temperature from the direct coexistence simulations of this work are in agreement (within the statistical uncertainty) with those obtained previously by us from free energy calculations. By taking into account the results of this work and those of our free energy calculations, recommended values of the melting point of ice I(h) at 1 bar for the above mentioned water models are provided.