E. Noya | CSIC (Consejo Superior de Investigaciones Científicas-Spanish National Research Council) (original) (raw)

Papers by E. Noya

Journal of Physical Chemistry C, 2007

A relatively simple equation of state is proposed for several forms of ice, whose parameters have... more A relatively simple equation of state is proposed for several forms of ice, whose parameters have been fitted to the results of extensive computer simulations using the TIP4P/Ice and TIP4P/2005 models of water. Comparison with available experimental data for ice Ih shows that both models reproduce the experimental density and isothermal compressibility to good accuracy over the entire range of thermodynamic stability except at low temperatures. The predictions for the thermal expansion coefficient are slightly worse but still reasonable. Results obtained with the TIP4P/2005 model are slightly better than those obtained with the TIP4P/Ice model. At temperatures below 150 K, the predictions of both models deviate significantly from experiment. As expected, at low temperatures, quantum effects become increasingly important, and classical simulations are unable to accurately describe the properties of ices. In fact, neither the heat capacity nor the thermal expansion coefficient go to zero at zero temperature (as they should be according to the third law of thermodynamics). Predicted compressibilites are however reliable even up to 0 K. Finally, the relative energies of the ices at 0 K have also been estimated and compared with the experiments. † Part of the "Keith E. Gubbins Festschrift".

The Journal of Physical Chemistry B, 2010

With the objective of establishing the importance of water flexibility in empirical models which ... more With the objective of establishing the importance of water flexibility in empirical models which explicitly include nuclear quantum effects, we have carried out path integral Monte Carlo simulations in water clusters with up to seven molecules. Two recently developed models have been used for comparison: the rigid TIP4PQ/ 2005 and the flexible q-TIP4P/F models, both inspired by the rigid TIP4P/2005 model. To obtain a starting configuration for our simulations, we have located the global minima for the rigid TIP4P/2005 and TIP4PQ/ 2005 models and for the flexible q-TIP4P/F model. All the structures are similar to those predicted by the rigid TIP4P potential showing that the charge distribution mainly determines the global minimum structure. For the flexible q-TIP4P/F model, we have studied the geometrical distortion upon isotopic substitution by studying tritiated water clusters. Our results show that tritiated water clusters exhibit an r OT distance shorter than the r OH distance in water clusters, not significant changes in the Φ HOH angle, and a lower average dipole moment than water clusters. We have also carried out classical simulations with the rigid TIP4PQ/2005 model showing that the rotational kinetic energy is greatly affected by quantum effects, but the translational kinetic energy is only slightly modified. The potential energy is also noticeably higher than in classical simulations. Finally, as a concluding remark, we have calculated the formation energies of water clusters using both models, finding that the formation energies predicted by the rigid TIP4PQ/2005 model are lower by roughly 0.6 kcal/mol than those of the flexible q-TIP4P/F model for clusters of moderate size, the origin of this difference coming mainly from the geometrical distortion of the water molecule in the clusters that causes an increase in the intramolecular potential energy.

The phase diagram of water has been calculated for the TIP4PQ/2005 model, an empirical rigid non-... more The phase diagram of water has been calculated for the TIP4PQ/2005 model, an empirical rigid non-polarisable model. The path integral Monte Carlo technique was used, permitting the incorporation of nuclear quantum effects. The coexistence lines were traced out using the Gibbs-Duhem integration method, once having calculated the free energies of the liquid and solid phases in the quantum limit, which were obtained via thermodynamic integration from the classical value by scaling the mass of the water molecule. The resulting phase diagram is qualitatively correct, being displaced to lower temperatures by 15-20K. It is found that the influence of nuclear quantum effects are correlated to the tetrahedral order parameter.

Physical Chemistry Chemical Physics, 2012

The melting point of ice I h , as well as the temperature of maximum density (TMD) in the liquid ... more The melting point of ice I h , as well as the temperature of maximum density (TMD) in the liquid phase, has been computed using the path integral Monte Carlo method. Two new models are introduced; TIP4PQ D2O and TIP4PQ T2O which are specifically designed to study D 2 O and T 2 O respectively. We have also used these models to study the "competing quantum effects" proposal of Habershon, Markland and Manolopoulos; the TIP4PQ/2005, TIP4PQ/2005 (D 2 O) and TIP4PQ/2005 (T 2 O) models are able to study the isotopic substitution of hydrogen for deuterium or tritium whilst constraining the geometry, while the TIP4PQ D2O and TIP4PQ T2O models, where the O-H bond lengths are progressively shortened, permit the study of the influence of geometry (and thus dipole moment) on the isotopic effects. For TIP4PQ D2O -TIP4PQ/2005 we found a melting point shift of 4.9 K (experimentally the value is 3.68K) and a TMD shift of 6K (experimentally 7.2K). For TIP4PQ T2O -TIP4PQ/2005 we found a melting point shift of 5.2 K (experimentally the value is 4.49K) and a TMD shift of 7K (experimentally 9.4K).

Physical Chemistry Chemical Physics, 2009

In this work the high pressure region of the phase diagram of water has been studied by computer ... more In this work the high pressure region of the phase diagram of water has been studied by computer simulation by using the TIP4P/2005 model of water. Free energy calculations were performed for ices VII and VIII and for the fluid phase to determine the melting curve of these ices. In addition molecular dynamics simulations were performed at high temperatures (440K) observing the spontaneous freezing of the liquid into a solid phase at pressures of about 80000 bar. The analysis of the structure obtained lead to the conclusion that a plastic crystal phase was formed.

Molecular Physics, 2011

In this work the path integral formulation for rigid rotors, proposed by Berne [Phys. Rev. Lett. ... more In this work the path integral formulation for rigid rotors, proposed by Berne [Phys. Rev. Lett. 77, 2638 (1996)], is described in detail. It is shown how this formulation can be used to perform Monte Carlo simulations of water. Our numerical results show that whereas some properties of water can be accurately reproduced using classical simulations with an empirical potential which, implicitly, includes quantum effects, other properties can only be described quantitatively when quantum effects are explicitly incorporated. In particular, quantum effects are extremely relevant when it comes to describing the equation of state of the ice phases at low temperatures, the structure of the ices at low temperatures, and the heat capacity of both liquid water and the ice phases. They also play a minor role in the relative stability of the ice phases. a eva.noya@iqfr.csic.es b cvega@quim.ucm.es ∞ J=0 J M =−J J K =−J

Molecular Physics, 2009

The so-called thermodynamic anomalies of water form an integral part of the peculiar behaviour of... more The so-called thermodynamic anomalies of water form an integral part of the peculiar behaviour of this both important and ubiquitous molecule. In this paper our aim is to establish whether the recently proposed TIP4P/2005 model is capable of reproducing a number of these anomalies. Using molecular dynamics simulations we investigate both the maximum in density and the minimum in the isothermal compressibility along a number of isobars. It is shown that the model correctly describes the decrease in the temperature of the density maximum with increasing pressure. At atmospheric pressure the model exhibits an additional minimum in density at a temperature of about 200K, in good agreement with recent experimental work on super-cooled confined water.

Journal of Physics: Condensed Matter, 2008

In this review we focus on the determination of phase diagrams by computer simulation with partic... more In this review we focus on the determination of phase diagrams by computer simulation with particular attention to the fluid-solid and solid-solid equilibria. The methodology to compute the free energy of solid phases will be discussed. In particular, the Einstein crystal and Einstein molecule methodologies are described in a comprehensive way. It is shown that both methodologies yield the same free energies and that free energies of solid phases present noticeable finite size effects. In fact this is the case for hard spheres in the solid phase. Finite size corrections can be introduced, although in an approximate way, to correct for the dependence of the free energy on the size of the system. The computation of free energies of solid phases can be extended to molecular fluids. The procedure to compute free energies of solid phases of water (ices) will be described in detail. The free energies of ices Ih, II, III, IV, V, VI, VII, VIII, IX, XI and XII will be presented for the SPC/E and TIP4P models of water. Initial coexistence points leading to the determination of the phase diagram of water for these two models will be provided. Other methods to estimate the melting point of a solid, as the direct fluid-solid coexistence or simulations of the free surface of the solid will be discussed. It will be shown that the melting points of ice Ih for several water models, obtained from free energy calculations, direct coexistence simulations and free surface simulations, agree within their statistical uncertainty. Phase diagram calculations can indeed help to improve potential models of molecular fluids. For instance, for water, the potential model TIP4P/2005 can be regarded as an improved version of TIP4P. Here we will review some recent work on the phase diagram of the simplest ionic model, the restricted primitive model. Although originally devised to describe ionic liquids, the model is becoming quite popular to describe the behaviour of charged colloids. Besides the possibility of obtaining fluid-solid equilibria for simple protein models will be discussed. In these primitive models, the protein is described by a spherical potential with certain anisotropic bonding sites (patchy sites).

The Journal of Chemical Physics, 2011

The expression for the quantum propagator for rigid tops, proposed by Müser and Berne [Phys. Rev.... more The expression for the quantum propagator for rigid tops, proposed by Müser and Berne [Phys. Rev. Lett. 77, 2638], has been extended to asymmetric tops. Path-integral Monte Carlo simulations are provided that show that the quantum propagator proposed in this work exactly reproduces the rotational energy of free asymmetric tops as evaluated from the partition function. This propagator can subsequently be used in path-integral simulations of condensed phases if a rigid molecular model is used. Noya, Vega, and McBride J. Chem. Phys. 134, 054117 (2011) ρ t,t+1 rot (β/P)

The Journal of Chemical Physics, 2007

The properties of ices I h , II, III, V, and VI at zero temperature and pressure are determined b... more The properties of ices I h , 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 I h 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.

The Journal of Chemical Physics, 2007

In this paper a new method to evaluate the free energy of solids is proposed. The method can be r... more In this paper a new method to evaluate the free energy of solids is proposed. The method can be regarded as a variant of the method proposed by Frenkel and Ladd ͓J. Chem. Phys. 81, 3188 ͑1984͔͒. The main equations of the method can be derived in a simple way. The method can be easily implemented within a Monte Carlo program. We have applied the method to determine the free energy of hard spheres in the solid phase for several system sizes. The obtained free energies agree within the numerical uncertainty with those obtained by Polson et al. ͓J. Chem. Phys. 112, 5339 ͑2000͔͒. The fluid-solid equilibria has been determined for several system sizes and compared to the values published previously by Wilding and Bruce ͓Phys. Rev. Lett. 85, 5138 ͑2000͔͒ using the phase switch methodology. It is shown that both the free energies and the coexistence pressures present a strong size dependence and that the results obtained from free energy calculations agree with those obtained using the phase switch method, which constitutes a cross-check of both methodologies. From the results of this work we estimate the coexistence pressure of the fluid-solid transition of hard spheres in the thermodynamic limit to be p * = 11.54͑4͒, which is slightly lower than the classical value of Hoover and Ree ͑p * = 11.70͒ ͓J. Chem. Phys. 49, 3609 ͑1968͔͒. Taking into account the strong size dependence of the free energy of the solid phase, we propose to introduce finite size corrections, which allow us to estimate approximately the free energy of the solid phase in the thermodynamic limit from the known value of the free energy of the solid phase with N molecules. We have also determined the free energy of a Lennard-Jones solid by using both the methodology of this work and the finite size correction. It is shown how a relatively good estimate of the free energy of the system in the thermodynamic limit is obtained even from the free energy of a relatively small system.

The Journal of Chemical Physics, 2007

We computed the phase diagram for a system of model anisotropic particles with six attractive pat... more We computed the phase diagram for a system of model anisotropic particles with six attractive patches in an octahedral arrangement. We chose to study this model for a relatively narrow value of the patch width where the lowest-energy configuration of the system is a simple cubic crystal. At this value of the patch width, there is no stable vapour-liquid phase separation, and there are three other crystalline phases in addition to the simple cubic crystal that is most stable at low pressure. Firstly, at moderate pressures, it is more favourable to form a body-centred cubic crystal, which can be viewed as two interpenetrating, and almost non-interacting, simple cubic lattices. Secondly, at high pressures and low temperatures, an orientationally ordered face-centred cubic structure becomes favourable. Finally, at high temperatures a face-centred cubic plastic crystal is the most stable solid phase.

The Journal of Chemical Physics, 2010

The phase diagram of model anisotropic particles with four attractive patches in a tetrahedral ar... more The phase diagram of model anisotropic particles with four attractive patches in a tetrahedral arrangement has been computed at two different values for the range of the potential, with the aim of investigating the conditions under which a diamond crystal can be formed. We find that the diamond phase is never stable for our longer-ranged potential. At low temperatures and pressures, the fluid freezes into a body-centred-cubic solid that can be viewed as two interpenetrating diamond lattices with a weak interaction between the two sublattices. Upon compression, an orientationally ordered face-centred-cubic crystal becomes more stable than the body-centred-cubic crystal, and at higher temperatures a plastic face-centered-cubic phase is stabilized by the increased entropy due to orientational disorder. A similar phase diagram is found for the shorter-ranged potential, but at low temperatures and pressures, we also find a region over which the diamond phase is thermodynamically favored over the body-centred-cubic phase. The higher vibrational entropy of the diamond structure with respect to the body-centred-cubic solid explains why it is stable even though the enthalpy of the latter phase is lower. Some preliminary studies on the growth of the diamond structure starting from a crystal seed were performed. Even though the diamond phase is never thermodynamically stable for the longer-ranged model, direct coexistence simulations of the interface between the fluid and the body-centred-cubic crystal and between the fluid and the diamond crystal show that, at sufficiently low pressures, it is quite probable that in both cases the solid grows into a diamond crystal, albeit involving some defects. These results highlight the importance of kinetic effects in the formation of diamond crystals in systems of patchy particles.

The Journal of Chemical Physics, 2010

Quantum path-integral simulations of the hydrate solid structures have been performed using the r... more Quantum path-integral simulations of the hydrate solid structures have been performed using the recently proposed TIP4PQ/2005 model. By also performing classical simulations using this model, the impact of the nuclear quantum effects on the hydrates is highlighted; nuclear quantum effects significantly modify the structure, densities, and energies of the hydrates, leading to the conclusion that nuclear quantum effects are important not only when studying the solid phases of water but also when studying the hydrates. To analyze the validity of a classical description of hydrates, a comparison of the results of the TIP4P/2005 model ͑optimized for classical simulations͒ with those of TIP4PQ/2005 ͑optimized for path-integral simulations͒ was undertaken. A classical description of hydrates is able to correctly predict the densities at temperatures above 150 K and the relative stabilities between the hydrates and ice I h . The inclusion of nuclear quantum effects does not significantly modify the sequence of phases found in the phase diagram of water at negative pressures, namely, I h → sII→ sH. In fact the transition pressures are little affected by the inclusion of nuclear quantum effects; the phase diagram predictions for hydrates can be performed with reasonable accuracy using classical simulations. However, for a reliable calculation of the densities below 150 K, the sublimation energies, the constant pressure heat capacity, and the radial distribution functions, the incorporation of nuclear quantum effects is indeed required.

The Journal of Chemical Physics, 2009

Path integral simulations have been performed to determine the temperature of the maximum in dens... more Path integral simulations have been performed to determine the temperature of the maximum in density of water of the rigid, nonpolarizable TIP4PQ/2005 model treating long range Coulombic forces with the reaction field method. A maximum in density is found at 280 K, just 3 K above the experimental value. In tritiated water the maximum occurs at a temperature about 12 K higher than in water, in reasonable agreement with the experimental result. Contrary to the usual assumption that the maximum in classical water is about 14 K above that in water, we found that for TIP4PQ/2005 this maximum is about 30 K above. For rigid water models the internal energy and the temperature of maximum density do not follow a linear behavior when plotted as a function of the inverse of the hydrogen mass. In addition, it is shown that, when used with Ewald sums, the TIP4PQ/2005 reproduces quite nicely not only the maximum in density of water, but also the liquid densities, the structure of liquid water and the vaporization enthalpy. It was shown in a previous work that it also reproduces reasonably well the density and relative stabilities of ices. Therefore the TIP4PQ/2005 model, while still simple, allows one to analyze the interplay between quantum effects related to atomic masses and intermolecular forces in water.

The Journal of Chemical Physics, 2009

With a view to a better understanding of the influence of atomic quantum delocalization effects o... more With a view to a better understanding of the influence of atomic quantum delocalization effects on the phase behavior of water, path integral simulations have been undertaken for almost all of the known ice phases using the TIP4P/2005 model in conjunction with the rigid rotor propagator proposed by Müser and Berne ͓Phys. Rev. Lett. 77, 2638 ͑1996͔͒. The quantum contributions then being known, a new empirical model of water is developed ͑TIP4PQ/2005͒ which reproduces, to a good degree, a number of the physical properties of the ice phases, for example, densities, structure, and relative stabilities.

The Journal of Chemical Physics, 2012

We present a local order parameter based on the standard Steinhardt-Ten Wolde approach that is ca... more We present a local order parameter based on the standard Steinhardt-Ten Wolde approach that is capable both of tracking and of driving homogeneous ice nucleation in simulations of all-atom models of water. We demonstrate that it is capable of forcing the growth of ice nuclei in supercooled liquid water simulated using the TIP4P/2005 model using over-biassed umbrella sampling Monte Carlo simulations. However, even with such an order parameter, the dynamics of ice growth in deeply supercooled liquid water in all-atom models of water are shown to be very slow, and so the computation of free energy landscapes and nucleation rates remains extremely challenging.

The Journal of Chemical Physics, 2013

In this work, we describe a procedure to evaluate the free energy of molecular solids with the GR... more In this work, we describe a procedure to evaluate the free energy of molecular solids with the GROMACS molecular dynamics package. The free energy is calculated using the Einstein molecule method that can be regarded as a small modification of the Einstein crystal method. Here, the position and orientation of the molecules is fixed by using an Einstein field that binds with harmonic springs at least three non-collinear atoms (or points of the molecule) to their reference positions. The validity of the Einstein field is tested by performing free-energy calculations of methanol, water (ice), and patchy colloids molecular solids. The free energies calculated with GROMACS show a very good agreement with those obtained using Monte Carlo and with previously published results. © 2013 AIP Publishing LLC. [http://dx.

The Journal of Chemical Physics, 2008

We consider the computation of the coexistence pressure of the liquid-solid transition of a syste... more We consider the computation of the coexistence pressure of the liquid-solid transition of a system of hard spheres from direct simulation of the inhomogeneous system formed from liquid and solid phases separated by an interface. Monte Carlo simulations of the interfacial system are performed in three different ensembles. In a first approach, a series of simulations is carried out in the isothermal-isobaric ensemble, where the solid is allowed to relax to its equilibrium crystalline structure, thus avoiding the appearance of artificial stress in the system. Here, the total volume of the system fluctuates due to changes in the three dimensions of the simulation box. In a second approach, we consider simulations of the inhomogeneous system in an isothermal-isobaric ensemble where the normal pressure, as well as the area of the ͑planar͒ fluid-solid interface, are kept constant. Now, the total volume of the system fluctuates due to changes in the longitudinal dimension of the simulation box. In both approaches, the coexistence pressure is estimated by monitoring the evolution of the density along several simulations carried out at different pressures. Both routes are seen to provide consistent values of the fluid-solid coexistence pressure, p = 11.54͑4͒k B T / 3 , which indicates that the error introduced by the use of the standard constant-pressure ensemble for this particular problem is small, provided the systems are sufficiently large. An additional simulation of the interfacial system is conducted in a canonical ensemble where the dimensions of the simulation box are allowed to change subject to the constraint that the total volume is kept fixed. In this approach, the coexistence pressure corresponds to the normal component of the pressure tensor, which can be computed as an appropriate ensemble average in a single simulation. This route yields a value of p = 11.54͑4͒k B T / 3 . We conclude that the results obtained for the coexistence pressure from direct simulations of the liquid and solid phases in coexistence using different ensembles are mutually consistent and are in excellent agreement with the values obtained from free energy calculations.

The Journal of Chemical Physics, 2007

We computed the complete phase diagram of the symmetrical colloidal electrolyte by means of Monte... more We computed the complete phase diagram of the symmetrical colloidal electrolyte by means of Monte Carlo simulations. Thermodynamic integration, together with the Einstein-crystal method, and Gibbs-Duhem integration were used to calculate the equilibrium phase behavior. The system was modeled via the linear screening theory, where the electrostatic interactions are screened by the presence of salt in the medium, characterized by the inverse Debye length, ͑in this work =6͒. Our results show that at high temperature, the hard-sphere picture is recovered, i.e., the liquid crystallizes into a fcc crystal that does not exhibit charge ordering. In the low temperature region, the liquid freezes into a CsCl structure because charge correlations enhance the pairing between oppositely charged colloids, making the liquid-gas transition metastable with respect to crystallization. Upon increasing density, the CsCl solid transforms into a CuAu-like crystal and this one, in turn, transforms into a tetragonal ordered crystal near close packing. Finally, we have studied the ordered-disordered transitions finding three triple points where the phases in coexistence are liquid-CsCl-disordered fcc, CsCl-CuAu-disordered fcc, and CuAu-tetragonal-disordered fcc.

Journal of Physical Chemistry C, 2007

A relatively simple equation of state is proposed for several forms of ice, whose parameters have... more A relatively simple equation of state is proposed for several forms of ice, whose parameters have been fitted to the results of extensive computer simulations using the TIP4P/Ice and TIP4P/2005 models of water. Comparison with available experimental data for ice Ih shows that both models reproduce the experimental density and isothermal compressibility to good accuracy over the entire range of thermodynamic stability except at low temperatures. The predictions for the thermal expansion coefficient are slightly worse but still reasonable. Results obtained with the TIP4P/2005 model are slightly better than those obtained with the TIP4P/Ice model. At temperatures below 150 K, the predictions of both models deviate significantly from experiment. As expected, at low temperatures, quantum effects become increasingly important, and classical simulations are unable to accurately describe the properties of ices. In fact, neither the heat capacity nor the thermal expansion coefficient go to zero at zero temperature (as they should be according to the third law of thermodynamics). Predicted compressibilites are however reliable even up to 0 K. Finally, the relative energies of the ices at 0 K have also been estimated and compared with the experiments. † Part of the "Keith E. Gubbins Festschrift".

The Journal of Physical Chemistry B, 2010

With the objective of establishing the importance of water flexibility in empirical models which ... more With the objective of establishing the importance of water flexibility in empirical models which explicitly include nuclear quantum effects, we have carried out path integral Monte Carlo simulations in water clusters with up to seven molecules. Two recently developed models have been used for comparison: the rigid TIP4PQ/ 2005 and the flexible q-TIP4P/F models, both inspired by the rigid TIP4P/2005 model. To obtain a starting configuration for our simulations, we have located the global minima for the rigid TIP4P/2005 and TIP4PQ/ 2005 models and for the flexible q-TIP4P/F model. All the structures are similar to those predicted by the rigid TIP4P potential showing that the charge distribution mainly determines the global minimum structure. For the flexible q-TIP4P/F model, we have studied the geometrical distortion upon isotopic substitution by studying tritiated water clusters. Our results show that tritiated water clusters exhibit an r OT distance shorter than the r OH distance in water clusters, not significant changes in the Φ HOH angle, and a lower average dipole moment than water clusters. We have also carried out classical simulations with the rigid TIP4PQ/2005 model showing that the rotational kinetic energy is greatly affected by quantum effects, but the translational kinetic energy is only slightly modified. The potential energy is also noticeably higher than in classical simulations. Finally, as a concluding remark, we have calculated the formation energies of water clusters using both models, finding that the formation energies predicted by the rigid TIP4PQ/2005 model are lower by roughly 0.6 kcal/mol than those of the flexible q-TIP4P/F model for clusters of moderate size, the origin of this difference coming mainly from the geometrical distortion of the water molecule in the clusters that causes an increase in the intramolecular potential energy.

The phase diagram of water has been calculated for the TIP4PQ/2005 model, an empirical rigid non-... more The phase diagram of water has been calculated for the TIP4PQ/2005 model, an empirical rigid non-polarisable model. The path integral Monte Carlo technique was used, permitting the incorporation of nuclear quantum effects. The coexistence lines were traced out using the Gibbs-Duhem integration method, once having calculated the free energies of the liquid and solid phases in the quantum limit, which were obtained via thermodynamic integration from the classical value by scaling the mass of the water molecule. The resulting phase diagram is qualitatively correct, being displaced to lower temperatures by 15-20K. It is found that the influence of nuclear quantum effects are correlated to the tetrahedral order parameter.

Physical Chemistry Chemical Physics, 2012

The melting point of ice I h , as well as the temperature of maximum density (TMD) in the liquid ... more The melting point of ice I h , as well as the temperature of maximum density (TMD) in the liquid phase, has been computed using the path integral Monte Carlo method. Two new models are introduced; TIP4PQ D2O and TIP4PQ T2O which are specifically designed to study D 2 O and T 2 O respectively. We have also used these models to study the "competing quantum effects" proposal of Habershon, Markland and Manolopoulos; the TIP4PQ/2005, TIP4PQ/2005 (D 2 O) and TIP4PQ/2005 (T 2 O) models are able to study the isotopic substitution of hydrogen for deuterium or tritium whilst constraining the geometry, while the TIP4PQ D2O and TIP4PQ T2O models, where the O-H bond lengths are progressively shortened, permit the study of the influence of geometry (and thus dipole moment) on the isotopic effects. For TIP4PQ D2O -TIP4PQ/2005 we found a melting point shift of 4.9 K (experimentally the value is 3.68K) and a TMD shift of 6K (experimentally 7.2K). For TIP4PQ T2O -TIP4PQ/2005 we found a melting point shift of 5.2 K (experimentally the value is 4.49K) and a TMD shift of 7K (experimentally 9.4K).

Physical Chemistry Chemical Physics, 2009

In this work the high pressure region of the phase diagram of water has been studied by computer ... more In this work the high pressure region of the phase diagram of water has been studied by computer simulation by using the TIP4P/2005 model of water. Free energy calculations were performed for ices VII and VIII and for the fluid phase to determine the melting curve of these ices. In addition molecular dynamics simulations were performed at high temperatures (440K) observing the spontaneous freezing of the liquid into a solid phase at pressures of about 80000 bar. The analysis of the structure obtained lead to the conclusion that a plastic crystal phase was formed.

Molecular Physics, 2011

In this work the path integral formulation for rigid rotors, proposed by Berne [Phys. Rev. Lett. ... more In this work the path integral formulation for rigid rotors, proposed by Berne [Phys. Rev. Lett. 77, 2638 (1996)], is described in detail. It is shown how this formulation can be used to perform Monte Carlo simulations of water. Our numerical results show that whereas some properties of water can be accurately reproduced using classical simulations with an empirical potential which, implicitly, includes quantum effects, other properties can only be described quantitatively when quantum effects are explicitly incorporated. In particular, quantum effects are extremely relevant when it comes to describing the equation of state of the ice phases at low temperatures, the structure of the ices at low temperatures, and the heat capacity of both liquid water and the ice phases. They also play a minor role in the relative stability of the ice phases. a eva.noya@iqfr.csic.es b cvega@quim.ucm.es ∞ J=0 J M =−J J K =−J

Molecular Physics, 2009

The so-called thermodynamic anomalies of water form an integral part of the peculiar behaviour of... more The so-called thermodynamic anomalies of water form an integral part of the peculiar behaviour of this both important and ubiquitous molecule. In this paper our aim is to establish whether the recently proposed TIP4P/2005 model is capable of reproducing a number of these anomalies. Using molecular dynamics simulations we investigate both the maximum in density and the minimum in the isothermal compressibility along a number of isobars. It is shown that the model correctly describes the decrease in the temperature of the density maximum with increasing pressure. At atmospheric pressure the model exhibits an additional minimum in density at a temperature of about 200K, in good agreement with recent experimental work on super-cooled confined water.

Journal of Physics: Condensed Matter, 2008

In this review we focus on the determination of phase diagrams by computer simulation with partic... more In this review we focus on the determination of phase diagrams by computer simulation with particular attention to the fluid-solid and solid-solid equilibria. The methodology to compute the free energy of solid phases will be discussed. In particular, the Einstein crystal and Einstein molecule methodologies are described in a comprehensive way. It is shown that both methodologies yield the same free energies and that free energies of solid phases present noticeable finite size effects. In fact this is the case for hard spheres in the solid phase. Finite size corrections can be introduced, although in an approximate way, to correct for the dependence of the free energy on the size of the system. The computation of free energies of solid phases can be extended to molecular fluids. The procedure to compute free energies of solid phases of water (ices) will be described in detail. The free energies of ices Ih, II, III, IV, V, VI, VII, VIII, IX, XI and XII will be presented for the SPC/E and TIP4P models of water. Initial coexistence points leading to the determination of the phase diagram of water for these two models will be provided. Other methods to estimate the melting point of a solid, as the direct fluid-solid coexistence or simulations of the free surface of the solid will be discussed. It will be shown that the melting points of ice Ih for several water models, obtained from free energy calculations, direct coexistence simulations and free surface simulations, agree within their statistical uncertainty. Phase diagram calculations can indeed help to improve potential models of molecular fluids. For instance, for water, the potential model TIP4P/2005 can be regarded as an improved version of TIP4P. Here we will review some recent work on the phase diagram of the simplest ionic model, the restricted primitive model. Although originally devised to describe ionic liquids, the model is becoming quite popular to describe the behaviour of charged colloids. Besides the possibility of obtaining fluid-solid equilibria for simple protein models will be discussed. In these primitive models, the protein is described by a spherical potential with certain anisotropic bonding sites (patchy sites).

The Journal of Chemical Physics, 2011

The expression for the quantum propagator for rigid tops, proposed by Müser and Berne [Phys. Rev.... more The expression for the quantum propagator for rigid tops, proposed by Müser and Berne [Phys. Rev. Lett. 77, 2638], has been extended to asymmetric tops. Path-integral Monte Carlo simulations are provided that show that the quantum propagator proposed in this work exactly reproduces the rotational energy of free asymmetric tops as evaluated from the partition function. This propagator can subsequently be used in path-integral simulations of condensed phases if a rigid molecular model is used. Noya, Vega, and McBride J. Chem. Phys. 134, 054117 (2011) ρ t,t+1 rot (β/P)

The Journal of Chemical Physics, 2007

The properties of ices I h , II, III, V, and VI at zero temperature and pressure are determined b... more The properties of ices I h , 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 I h 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.

The Journal of Chemical Physics, 2007

In this paper a new method to evaluate the free energy of solids is proposed. The method can be r... more In this paper a new method to evaluate the free energy of solids is proposed. The method can be regarded as a variant of the method proposed by Frenkel and Ladd ͓J. Chem. Phys. 81, 3188 ͑1984͔͒. The main equations of the method can be derived in a simple way. The method can be easily implemented within a Monte Carlo program. We have applied the method to determine the free energy of hard spheres in the solid phase for several system sizes. The obtained free energies agree within the numerical uncertainty with those obtained by Polson et al. ͓J. Chem. Phys. 112, 5339 ͑2000͔͒. The fluid-solid equilibria has been determined for several system sizes and compared to the values published previously by Wilding and Bruce ͓Phys. Rev. Lett. 85, 5138 ͑2000͔͒ using the phase switch methodology. It is shown that both the free energies and the coexistence pressures present a strong size dependence and that the results obtained from free energy calculations agree with those obtained using the phase switch method, which constitutes a cross-check of both methodologies. From the results of this work we estimate the coexistence pressure of the fluid-solid transition of hard spheres in the thermodynamic limit to be p * = 11.54͑4͒, which is slightly lower than the classical value of Hoover and Ree ͑p * = 11.70͒ ͓J. Chem. Phys. 49, 3609 ͑1968͔͒. Taking into account the strong size dependence of the free energy of the solid phase, we propose to introduce finite size corrections, which allow us to estimate approximately the free energy of the solid phase in the thermodynamic limit from the known value of the free energy of the solid phase with N molecules. We have also determined the free energy of a Lennard-Jones solid by using both the methodology of this work and the finite size correction. It is shown how a relatively good estimate of the free energy of the system in the thermodynamic limit is obtained even from the free energy of a relatively small system.

The Journal of Chemical Physics, 2007

We computed the phase diagram for a system of model anisotropic particles with six attractive pat... more We computed the phase diagram for a system of model anisotropic particles with six attractive patches in an octahedral arrangement. We chose to study this model for a relatively narrow value of the patch width where the lowest-energy configuration of the system is a simple cubic crystal. At this value of the patch width, there is no stable vapour-liquid phase separation, and there are three other crystalline phases in addition to the simple cubic crystal that is most stable at low pressure. Firstly, at moderate pressures, it is more favourable to form a body-centred cubic crystal, which can be viewed as two interpenetrating, and almost non-interacting, simple cubic lattices. Secondly, at high pressures and low temperatures, an orientationally ordered face-centred cubic structure becomes favourable. Finally, at high temperatures a face-centred cubic plastic crystal is the most stable solid phase.

The Journal of Chemical Physics, 2010

The phase diagram of model anisotropic particles with four attractive patches in a tetrahedral ar... more The phase diagram of model anisotropic particles with four attractive patches in a tetrahedral arrangement has been computed at two different values for the range of the potential, with the aim of investigating the conditions under which a diamond crystal can be formed. We find that the diamond phase is never stable for our longer-ranged potential. At low temperatures and pressures, the fluid freezes into a body-centred-cubic solid that can be viewed as two interpenetrating diamond lattices with a weak interaction between the two sublattices. Upon compression, an orientationally ordered face-centred-cubic crystal becomes more stable than the body-centred-cubic crystal, and at higher temperatures a plastic face-centered-cubic phase is stabilized by the increased entropy due to orientational disorder. A similar phase diagram is found for the shorter-ranged potential, but at low temperatures and pressures, we also find a region over which the diamond phase is thermodynamically favored over the body-centred-cubic phase. The higher vibrational entropy of the diamond structure with respect to the body-centred-cubic solid explains why it is stable even though the enthalpy of the latter phase is lower. Some preliminary studies on the growth of the diamond structure starting from a crystal seed were performed. Even though the diamond phase is never thermodynamically stable for the longer-ranged model, direct coexistence simulations of the interface between the fluid and the body-centred-cubic crystal and between the fluid and the diamond crystal show that, at sufficiently low pressures, it is quite probable that in both cases the solid grows into a diamond crystal, albeit involving some defects. These results highlight the importance of kinetic effects in the formation of diamond crystals in systems of patchy particles.

The Journal of Chemical Physics, 2010

Quantum path-integral simulations of the hydrate solid structures have been performed using the r... more Quantum path-integral simulations of the hydrate solid structures have been performed using the recently proposed TIP4PQ/2005 model. By also performing classical simulations using this model, the impact of the nuclear quantum effects on the hydrates is highlighted; nuclear quantum effects significantly modify the structure, densities, and energies of the hydrates, leading to the conclusion that nuclear quantum effects are important not only when studying the solid phases of water but also when studying the hydrates. To analyze the validity of a classical description of hydrates, a comparison of the results of the TIP4P/2005 model ͑optimized for classical simulations͒ with those of TIP4PQ/2005 ͑optimized for path-integral simulations͒ was undertaken. A classical description of hydrates is able to correctly predict the densities at temperatures above 150 K and the relative stabilities between the hydrates and ice I h . The inclusion of nuclear quantum effects does not significantly modify the sequence of phases found in the phase diagram of water at negative pressures, namely, I h → sII→ sH. In fact the transition pressures are little affected by the inclusion of nuclear quantum effects; the phase diagram predictions for hydrates can be performed with reasonable accuracy using classical simulations. However, for a reliable calculation of the densities below 150 K, the sublimation energies, the constant pressure heat capacity, and the radial distribution functions, the incorporation of nuclear quantum effects is indeed required.

The Journal of Chemical Physics, 2009

Path integral simulations have been performed to determine the temperature of the maximum in dens... more Path integral simulations have been performed to determine the temperature of the maximum in density of water of the rigid, nonpolarizable TIP4PQ/2005 model treating long range Coulombic forces with the reaction field method. A maximum in density is found at 280 K, just 3 K above the experimental value. In tritiated water the maximum occurs at a temperature about 12 K higher than in water, in reasonable agreement with the experimental result. Contrary to the usual assumption that the maximum in classical water is about 14 K above that in water, we found that for TIP4PQ/2005 this maximum is about 30 K above. For rigid water models the internal energy and the temperature of maximum density do not follow a linear behavior when plotted as a function of the inverse of the hydrogen mass. In addition, it is shown that, when used with Ewald sums, the TIP4PQ/2005 reproduces quite nicely not only the maximum in density of water, but also the liquid densities, the structure of liquid water and the vaporization enthalpy. It was shown in a previous work that it also reproduces reasonably well the density and relative stabilities of ices. Therefore the TIP4PQ/2005 model, while still simple, allows one to analyze the interplay between quantum effects related to atomic masses and intermolecular forces in water.

The Journal of Chemical Physics, 2009

With a view to a better understanding of the influence of atomic quantum delocalization effects o... more With a view to a better understanding of the influence of atomic quantum delocalization effects on the phase behavior of water, path integral simulations have been undertaken for almost all of the known ice phases using the TIP4P/2005 model in conjunction with the rigid rotor propagator proposed by Müser and Berne ͓Phys. Rev. Lett. 77, 2638 ͑1996͔͒. The quantum contributions then being known, a new empirical model of water is developed ͑TIP4PQ/2005͒ which reproduces, to a good degree, a number of the physical properties of the ice phases, for example, densities, structure, and relative stabilities.

The Journal of Chemical Physics, 2012

We present a local order parameter based on the standard Steinhardt-Ten Wolde approach that is ca... more We present a local order parameter based on the standard Steinhardt-Ten Wolde approach that is capable both of tracking and of driving homogeneous ice nucleation in simulations of all-atom models of water. We demonstrate that it is capable of forcing the growth of ice nuclei in supercooled liquid water simulated using the TIP4P/2005 model using over-biassed umbrella sampling Monte Carlo simulations. However, even with such an order parameter, the dynamics of ice growth in deeply supercooled liquid water in all-atom models of water are shown to be very slow, and so the computation of free energy landscapes and nucleation rates remains extremely challenging.

The Journal of Chemical Physics, 2013

In this work, we describe a procedure to evaluate the free energy of molecular solids with the GR... more In this work, we describe a procedure to evaluate the free energy of molecular solids with the GROMACS molecular dynamics package. The free energy is calculated using the Einstein molecule method that can be regarded as a small modification of the Einstein crystal method. Here, the position and orientation of the molecules is fixed by using an Einstein field that binds with harmonic springs at least three non-collinear atoms (or points of the molecule) to their reference positions. The validity of the Einstein field is tested by performing free-energy calculations of methanol, water (ice), and patchy colloids molecular solids. The free energies calculated with GROMACS show a very good agreement with those obtained using Monte Carlo and with previously published results. © 2013 AIP Publishing LLC. [http://dx.

The Journal of Chemical Physics, 2008

We consider the computation of the coexistence pressure of the liquid-solid transition of a syste... more We consider the computation of the coexistence pressure of the liquid-solid transition of a system of hard spheres from direct simulation of the inhomogeneous system formed from liquid and solid phases separated by an interface. Monte Carlo simulations of the interfacial system are performed in three different ensembles. In a first approach, a series of simulations is carried out in the isothermal-isobaric ensemble, where the solid is allowed to relax to its equilibrium crystalline structure, thus avoiding the appearance of artificial stress in the system. Here, the total volume of the system fluctuates due to changes in the three dimensions of the simulation box. In a second approach, we consider simulations of the inhomogeneous system in an isothermal-isobaric ensemble where the normal pressure, as well as the area of the ͑planar͒ fluid-solid interface, are kept constant. Now, the total volume of the system fluctuates due to changes in the longitudinal dimension of the simulation box. In both approaches, the coexistence pressure is estimated by monitoring the evolution of the density along several simulations carried out at different pressures. Both routes are seen to provide consistent values of the fluid-solid coexistence pressure, p = 11.54͑4͒k B T / 3 , which indicates that the error introduced by the use of the standard constant-pressure ensemble for this particular problem is small, provided the systems are sufficiently large. An additional simulation of the interfacial system is conducted in a canonical ensemble where the dimensions of the simulation box are allowed to change subject to the constraint that the total volume is kept fixed. In this approach, the coexistence pressure corresponds to the normal component of the pressure tensor, which can be computed as an appropriate ensemble average in a single simulation. This route yields a value of p = 11.54͑4͒k B T / 3 . We conclude that the results obtained for the coexistence pressure from direct simulations of the liquid and solid phases in coexistence using different ensembles are mutually consistent and are in excellent agreement with the values obtained from free energy calculations.

The Journal of Chemical Physics, 2007

We computed the complete phase diagram of the symmetrical colloidal electrolyte by means of Monte... more We computed the complete phase diagram of the symmetrical colloidal electrolyte by means of Monte Carlo simulations. Thermodynamic integration, together with the Einstein-crystal method, and Gibbs-Duhem integration were used to calculate the equilibrium phase behavior. The system was modeled via the linear screening theory, where the electrostatic interactions are screened by the presence of salt in the medium, characterized by the inverse Debye length, ͑in this work =6͒. Our results show that at high temperature, the hard-sphere picture is recovered, i.e., the liquid crystallizes into a fcc crystal that does not exhibit charge ordering. In the low temperature region, the liquid freezes into a CsCl structure because charge correlations enhance the pairing between oppositely charged colloids, making the liquid-gas transition metastable with respect to crystallization. Upon increasing density, the CsCl solid transforms into a CuAu-like crystal and this one, in turn, transforms into a tetragonal ordered crystal near close packing. Finally, we have studied the ordered-disordered transitions finding three triple points where the phases in coexistence are liquid-CsCl-disordered fcc, CsCl-CuAu-disordered fcc, and CuAu-tetragonal-disordered fcc.