Molecular Simulations of the Vapor–Liquid Phase Interfaces of Pure Water Modeled with the SPC/E and the TIP4P/2005 Molecular Models

In our previous study [Planková et al., EPJ Web. Conf. 92, 02071 (2015)], several molecular simulations of vapor-liquid phase interfaces for pure water were performed using the DL POLY Classic software. The TIP4P/2005 molecular model was successfully used for the modeling of the density profile and the thickness of phase interfaces together with the temperature dependence of the surface tension. In the current study, the extended simple point charge (SPC/E) model for water was employed for the investigation of vapor-liquid phase interfaces over a wide temperature range from 250 K to 600 K. The TIP4P/2005 model was also used with the temperature step of 25 K to obtain more consistent data compared to our previous study. Results of the new simulations are in a good agreement with most of the literature data. TIP4P/2005 provides better results for the saturated liquid density with its maximum close to 275 K, while SPC/E predicts slightly better saturated vapor density. Both models give qualitatively correct representation for the surface tension of water. The maximum absolute deviation from the IAPWS standard for the surface tension of ordinary water is 10.4 mN ·m−1 and 4.1 mN ·m−1 over the temperature range from 275 K to 600 K in case of SPC/E and TIP4P/2005, respectively.


Introduction
Molecular simulations (MS) represent a computational technique that allows to study a large variety of physical and chemical processes on the molecular scale.MS works particularly well for fluid systems as the liquid phase and the gas phase can be treated simultaneously.The only difference between the phases is given by the local density.
There are two different methods of MS, the molecular dynamics (MD) and Monte Carlo.MD solves a simple Newton equation of motion for every atom in discretized time-steps while Monte Carlo uses the probability nature of the statistical mechanics.Although all the methods in MS stand on relatively simple basis, their more wide-spread use has emerged only in recent years.The main limitation of all MS methods are their enormous computational demands.However, as computers are getting faster and more powerful, MS can gradually solve larger systems over longer time periods and more complex physical models.
One of the most important substances investigated using the approaches of MS is water.Water is essential in daily life, industry, and biological and natural processes.However, it has a typical non-standard behavior and many anomalies, which are still hard to model theoretically.The non-trivial phenomena are caused especially by its polar character and by the formation of hydrogen bonds.
The pioneers in successful simulations of water were Barker and Watts [1] and Rahman and Stillinger [2] around 1970.For many problems which include modeling complex aqueous systems, relatively simple water models such as a e-mail: : barbora.plankova@gmail.comSPC or TIPXP-type (TIP3P, TIP4P, or TIP5P) can be used.Currently, the modifications of the original SPC and TIP4P models are considered as the most convenient water models; namely the SPC/E model [3,4] and the TIP4P/2005 model [5,6].
In this study, we continue MD simulations of water [7] focused on the investigation of properties at the vapor-liquid phase interface.A series of MS was performed with both the SPC/E and TIP4P/2005 models in order to evaluate the density profiles, the thicknesses of the phase interface, and the surface tensions.A "slab" geometry with water molecules was simulated at various temperatures from the triple point of water to the temperatures close to the critical point.The conditions and settings of MS were inspired by the works of Vega and Miguel [8] and Sakamaki et al. [9].All simulations were performed using DL POLY Classic 1.9 [10,11] running on a cluster with 4 × 2 Intel(R) Xeon(R) CPU E5645 (12 cores effectively for each CPU) @ 2.53GHz CPUs.New results were compared to both MS from the literature and the experimental data.

Water models
Despite the fact that water molecule is flexible and polarizable, simple rigid models are used regularly in many MS.Therefore, systems with large numbers of molecules can be modeled for a longer times.The most widespread simple models belong to the SPC and TIPXP families.A common feature of these models is that they locate the Lennard-Jones interaction site on the position of the oxygen and the positive charge on the hydrogen atoms.Nevertheless, they differ in the target properties, i.e. the charge distribution and the angle of H-O-H bonds.The basic three-site models TIP3P [12] and SPC [13] were developed in 1980's.An extended modification of the simple bond charge, i.e. the SPC/E model by Berendsen et al. [14], including empirical corrections to the polarization energy is very successful and reproduces vapor-liquid data very well.A modification of TIP3P is the TIP4P model, which tries to improve the electrostatic distribution within the water molecule [12] by adding an auxiliary atom.This atom is located close to the oxygen atom (only 0.15 Å distant), it has no mass and carries a negative charge.A further modification, called TIP4P/2005, was proposed by Vega and Abascal [15] in 2005.They tried to combine the good phase diagram of TIP4P with target properties of SPC/E improving the melting point.There are many other modifications of the TIP4P model, e.g., the extension TIP5P by Mahoney et al. [16] carrying the negative charge on two auxiliary atoms.However, the general performance of most of these models is not better than TIP4P/2005.In the present study, the two models TIP4P/2005 and SPC/E are employed to investigate vapor-liquid interfacial properties of water.Both the TIP4P/2005 and SPC/E models consider the Lennard-Jones interactions only between the oxygen atoms, because hydrogen atoms have a negligible mass compared to oxygen.The electrostatic interaction occurs between the hydrogen atoms and oxygen atom in the SPC/E model and between hydrogen atoms and the massless atom Q with the negative charge in the TIP4P/2005 model, respectively.A schematic representation of the SPC/E and TIP4P/2005 molecular models is shown in figure 1.The parameters of the models are summarized in Table 1.Both models provide a very good description of the vapor-liquid equilibria compared to the other molecular models popular for water, e.g., SPC, TIP3P, TIP4P/Ice [8].Moreover, TIP4P/2005 gives a reasonable description also for the compressibility, the maximum in liquid density (near room temperature), and thermal expansion coefficient.

Simulation methods
A series of MD simulations was performed with the NVT ensemble, i.e. for a constant number of molecules, a constant total volume of the system, and a constant temperature.The ensemble consisted of 1372 molecules, which is considerably more than in other studies, e.g., both Vega and Miguel [8] and Sakamaki et al. [9] considered  [17] using the reference equation of state for water by Wagner and Pruss [18].The liquid densities of the SPC/E and TIP4P/2005 models significantly differ from the experimental values at high temperatures close to the critical point.Consequently, at temperatures above 450 K, the box size had to be evaluated from an estimate of the liquid density from previous MS and not from the reference equation of state.After a homogeneous liquid phase was formed in the cubic box, the z-size of the simulation box was symmetrically expanded by a factor of 3. As a result, the liquid phase forms a slab in the centre of the simulation box and the expanded empty volumes on each side of the slab allow subsequent formation of the vapor phase.The slab configuration was simulated for 10 ns in all cases.An additional equilibration was not required as the liquid phase formed from the cubic box had convenient initial physical conditions.A constant time step of 2 fs was considered in all simulations.This value was also used in similar simulations by Chen and Smith [4] and Sakamaki et al. [9].An initial configuration of molecules within the cubic box had to be prepared as an input for the DL POLY program before the simulation.The configuration was prepared using our code developed in MATLAB [19].The initial cubic box was generated as a grid of oxygen atoms with face-centered cubic unit cell.The remaining atoms of the water molecules, i.e. the hydrogen atoms and in case of TIP4P/2005 also the massless charge atom Q, were afterwards placed with consideration of the molecular structure (figure 1) and then randomly rotated.A thin void of 1 Å was added around the cubic box on each side to prevent overlapping of atoms from the simulation box randomly.The last configuration of the simulation within the cubic box was subsequently used as the input configuration for the slab simulation.In both simulations (in the cubic box and in the slab geometry), the periodic boundary conditions were imposed on all walls of the simulation box.However, the periodic boundary conditions can cause problems when extending the simulation box in the z-direction, because the molecules close to the walls whose atoms are periodically reflected on the other side of the cubic box can split up.Consequently, each molecule influenced by the periodic conditions on the boundaries had to be joint together after the expansion of the simulation box such that the molecule would not be pulled apart.We developed another MATLAB code for this purpose.
The MD simulations do not fix the temperature implicitly, therefore a thermostat is needed to correct kinetic energies of the molecules and thus regulate the temperature of the system.For this purpose, we used the Nosé-Hoover thermostat with the relaxation constant of 100 fs.The Verlet velocity integrator was used.The cutoffs of the Lennard-Jones interactions and the van der Waals forces were set to 14.5 Å.The direct Ewald method with an automatic parameter optimization constant set to 10 −5 was used for electrostatic interactions.This value gave a convergence parameter of 0.192 Å −1 and k-vector indexes of 6 × 6 × 6 for the cubic box and 6 × 6 × 18 for the slab geometry.For higher temperatures closer to the critical point, k-vector indexes increased up to 7 × 7 × 7 and 7 × 7 × 29, respectively.

Results
A series of 15 simulations at various temperatures spaced by 25 K was performed for each of the two water models.One simulation of the slab geometry running on 22 CPU cores took on average 48 hours in the case of the SPC/E model and 62 hours in the case of the TIP4P/2005 model.The difference in the computational time is caused by the additional massless atom Q that needed to be included in the TIP4P/2005 model.
The density profile ρ(z) at a given temperature was computed as an average from the histograms of positions of water molecules in the z-direction calculated at each time step.Properties of the phase interface can, according to Chapela et al. [20], be determined from the correlation of the averaged density profile by a hyperbolic tangent function given in the following manner In equation ( 1), z marks the coordinate perpendicular to the phase interface, ρ(z) denotes the density in Å −3 and z 0 , d, ρ L , and ρ V are adjustable parameters.z 0 is the position of the Gibbs dividing surface of the interface, d is a parameter of the phase interface thickness, ρ L denotes the saturated liquid density, and ρ V is the saturated vapor density.We note that Chapela et al. [20] used the term 2(z − z 0 )/d in the hyperbolic tangent, while the subsequent studies, e.g., [3,8,9], considered only (z − z 0 )/d.To assure compatibility with other MS studies of water, we decided to use the later definition in equation (1) as well.
Figure 2 shows the saturated liquid and vapor densities obtained with the SPC/E model compared to other MS and to experimental data.The experimental data are represented by predictions of the reference equation of state IAPWS-95 [18] obtained from the TREND package [17].The new simulations are in good agreement with other studies [3,8,9].The saturated vapor density ρ V obtained with the SPC/E model is in a quite good agreement with the IAPWS-95 this study with the literature data [8,9] and the reference equation of state IAPWS-95 [18].
equation of state, while the saturated liquid density ρ L shows a clear deviation.The SPC/E liquid density is lower by 0.6% and 5.3% than the experimental data at temperatures 250 K and 500 K, respectively.We note that the difference becomes less pronounced when comparing the densities related to the reduced temperature T r = T/T c .
As will be discussed later, near Table 2, both molecular models have lower critical temperature T c than the physically correct value employed in IAPWS-95.The melting point temperature of both water models considered in this work is below the real triple point of water 273.16K.The triple point values are 213 K for the SPC/E model and 249 K for the TIP4P/2005 model [9].Consequently, the vapor-liquid equilibrium could be calculated down to 250 K for both models.
The densities for the TIP4P/2005 model are shown in figure 3. The saturated liquid density obtained with TIP4P/ EFM 2015 02136-p.3 2005 lies closer to IAPWS-95 than the SPC/E model data (see figures 2 and 3).The TIP4P/2005 liquid density is lower only by 0.2% and 1.9% than the experimental data at temperatures 250 K and 500 K, respectively.Moreover, the TIP4P/2005 model provides relatively good representation of the maximum in the density of liquid water at 277 K.The maximum of approximately 0.999 g/cm 3 lies close to 275 K in our simulations.The saturated vapor density obtained with TIP4P/2005 is slightly worse than in the case of SPC/E.
The temperature dependence of the 10-90 thickness of the phase interface t related to thickness parameter d from equation (1) as t = 2.1972d is shown in figure 4. The 10-90 thicknesses calculated with SPC/E and TIP4P/2005 are in good mutual agreement over the whole temperature range from 250 K to 600 K.The results of the new simulations are compared to other MS and the experimental data by Matsumoto and Kataoka [21].As can be seen, the new results agree with the previous simulations at low temperatures below 400 K.However, with the increasing temperature, the thickness from the new simulations becomes considerably higher than in other MS, studies which are in relatively good mutual agreement over the entire temperature range.The reason for this difference might be the higher amount of molecules.Our ensemble consists of 1372 molecules compared to Alejandre et al. [3] with 1000 molecules or Vega and Miguel [8] with 1020 molecules in ensemble.As will be shown later in Table 3, number of molecules in the simulated ensemble has a systematic influence on the calculated results for the saturated liquid density and the interface thickness.Results of the new simulation compared to the literature data [3,8] and experimental data [21].
The experimental data for the interface thickness t that is shown in figure 4 were obtained by Matsumoto and Kataoka [21] using elipsometry, a measuring technique that allows to determine information about the interface thickness and the dielectric constant from the polarization of the light reflected at the phase interface.As can be seen in figure 4, the MS agree well with the experimental data at temperatures below 300 K.At higher temperatures of 350 K and 450 K, both the new simulations and the previous MS studies underpredict the experimental data for t by several Ångströms.
The surface tension was computed in a conventional way as a sum of the pressure tensor (or initial) term γ pt and the tail correction to the long-range Lennard-Jones interactions γ tail that corrects the cutoff truncation error, ( The pressure tensor term can be determined in the following manner where L z denotes the box size in z direction, P ii is the ii-th diagonal component of the pressure tensor.The tail correction γ tail is given by where and σ are the Lennard-Jones parameters for oxygen atom, r c is the cutoff for the Lennard-Jones potential and ρ L and ρ V are the bulk densities determined from equation (1).[3,4,8,9] and the IAPWS standard [22].Dash-anddotted black line and dashed red line are results of relation (5) with parameters taken from Table 2 for Vega and Miguel [8] and this study, respectively.
The results for the surface tension obtained with the SPC/E model are shown in figure 5.The new simulations are in a very good agreement with the previous studies by Chen and Smith [4], Vega and Miguel [8], and Sakamaki et al. [9].The surface tension predicted in 1995 by Alejandre et al. [3] is much higher compared to other MS studies and lies remarkably close to the experimental data represented by the IAPWS standard for the surface tension of ordinary water [22].Moreover, the data by Alejandre et al. [3] show relatively large scatter compared to other studies.Alejandre et al. [3] performed rather short simulations with 250 ps EPJ Web of Conferences 02136-p.4equilibration time and subsequent production runs of 250 ps and 125 ps for samples with 512 and 1000 molecules, respectively.The time step was set to 2.5 fs, i.e. a similar value as in other simulations.However, the overall simulation time between 375 and 500 ps seems to be insufficient in order to obtain correct and reproducible results.For example, Vega and Miguel [8] stated that a minimum time for the surface tension simulations is above 0.5 ns.They performed their simulations for 1.5 to 2.0 ns, which is still at least five times less than in the case of Sakamaki et al. [9] and this study.Chen and Smith [4] concluded that a minimum time for an accurate surface tension simulations is around 2 to 5 ns.The longer simulation time of 10 ns in the new simulations may also explain the differences in the increasing 10-90 interface thickness t shown in figure 4. Unfortunately, Sakamaki et al. [9] did not provide the results for t from their 10 ns simulations, which would be an interesting comparison with our results.
The TIP4P/2005 model gives a better prediction for the surface tension of water compared to the SPC/E model.The difference between the models can be seen from the deviations of the SPC/E and TIP4P/2005 results from the black line representing the IAPWS standard [22] in figures 5 and 6.The new results for TIP4P/2005 are in good agreement with our previous simulations [7] and the results by Vega and Miguel [8] and Sakamaki et al. [9].[7][8][9] and the IAPWS standard [22].Dashand-dotted black line and dashed red line are results of correlation (5) with parameters taken from Table 2 for Vega and Miguel [8] and this study, respectively.
The surface tension data shown in figures 5 and 6 were correlated by the following equation where B and b are adjustable parameters and T c is the critical temperature of the molecular model.We continue using this equation for comparison purposes.We note that Eq. ( 5) is of the same form as the IAPWS standard [22], except for the exponent, whose IAPWS value 1.256 is    [9] and this study.We note that values in Table 2 are rather sensitive to the amount and quality of the correlated data points.The values for T c , B, and b given in the IAPWS standard [22] are also given in Table 2 for comparison.
The influence of the ensemble size was thoroughly investigated with the SPC/E model at a constant temperature of 350 K. Four simulations were performed with the numbers of molecules varying from 500 to 2048. Figure 7 shows the density profiles close to the centre of the simulation box for all four simulations.The thickness of the EFM 2015 02136-p.5 liquid slab is continuously increasing as the amount of water molecules gets higher.Table 3 contains results for the saturated liquid and vapor densities, the thickness of the phase interface, and the surface tension.The liquid density slightly increases in a larger ensemble.The vapor density and the surface tension show small fluctuations with no obvious dependency on the number of simulated molecules.The only significant influence of the ensemble size can be seen for the 10-90 thickness of the phase interface t, which continuously increases approximately by 17% with increasing the number of molecules from 500 to 2048.The observed dependence of the interface thickness on the ensemble size can also explain differences in t shown in figure 4. The new simulations were made with 1372 molecules, while Alejandre et al. [3] and Vega and Miguel [8] used a maximum number of molecules of 1000 and 1024, respectively.

Conclusions
Properties of the vapor-liquid phase interfaces of pure water were studied using the approaches of MD.Two commonly used molecular models of water, SPC/E and TIP4P/ 2005, were used to model the density profiles and the surface tensions of water at temperatures ranging from 250 K to 600 K. Results of our new simulations are in good agreement with most of the previous MS studies.Both models were found to provide relatively good results for the vapor-liquid phase equilibria and the densities of water at temperatures up to 500 K.At higher temperatures, both models start to deviate from the experimental data as they underpredict the critical point of water.SPC/E gives slightly better reproduction of the vapor density while TIP4P/2005 is more accurate on the liquid side.The TIP4P/2005 model gives higher values for the surface tension lying closer to the experimental data correlated by IAPWS [22] than the SPC/E model.An influence of the number of molecules in the ensemble on the 10-90 thickness of the phase interface was observed.The interface thickness increases with the larger amount of molecules in the simulation.
On the basis of the previous studies by Vega and Miguel [8] and Sakamaki et al. [9] and our new results we conclude that the TIP4P/2005 molecular model is convenient for modelling the interfacial properties and surface tension of pure water in the temperature range from 250 K to 550 K, although a fully quantitative prediction cannot be expected of these relatively simple models.

DOI: 10
.1051/ C Owned by the authors, published by EDP Sciences,

Fig. 2 .Fig. 3 .
Fig.2.Saturated liquid and vapor densities for water obtained from MS with the SPC/E water model.Comparison of results from this study with the literature data[3,8,9] and the reference equation of state IAPWS-95[18].

Fig. 4 .
Fig. 4. The 10-90 thickness of the interface given as t = 2.1972d, where the thickness parameter d was determined from relation (1).Results of the new simulation compared to the literature data[3,8] and experimental data[21].

Fig. 5 .
Fig.5.Surface tension of water obtained from MS with the SPC/E water model.Comparison of results from this study with the literature data[3,4,8,9] and the IAPWS standard[22].Dash-anddotted black line and dashed red line are results of relation(5) with parameters taken from Table2for Vega and Miguel[8] and this study, respectively.

Fig. 6 .
Fig.6.Surface tension of water obtained from MS with the TIP4P/2005 water model.Comparison of results from this study with the literature data[7][8][9] and the IAPWS standard[22].Dashand-dotted black line and dashed red line are results of correlation(5) with parameters taken from Table2for Vega and Miguel[8] and this study, respectively.

Fig. 7 .
Fig. 7. Density profiles of the SPC/E model at a constant temperature 350 K for simulations with various amount of water molecules N. Black points are results of MS, lines represent correlation (1).
[9].Parameters B and b from the new simulations lie between the values of the other two MS studies.The values for T c , B, and b for the TIP4P/2005 model are very close for all three studies.The only difference can be seen in the critical temperature, where the value by Vega and Miguel[8] is 5 K above by Sakamaki et al.

Table 1 .
Parameters of the molecular models SPC/E and TIP4P/2005 for water: oxygen O, hydrogen H, massless atom with charge Q. M denotes molar mass, q is charge in units of elementary charge, and σ are Lennard-Jones parameters.
in other studies, the size of the cubic box was not constant as it varied with temperature.The size of the box was determined from the temperature dependence of the saturated liquid density calculated from the TREND package

Table 2 .
Critical temperature T c and coefficients B and b from equation (5) correlated to the surface tension data obtained from MS. * exponent has different value of 1.256 in the IAPWS standard closer to the current estimate 1.260 of the universal critical exponent [23].The values for T c , B, and b are given in * correlated in this study to the data by Sakamaki et al.[9]*

Table 2 .
[8] critical temperature for SPC/E is approximately 4 K lower than the values provided by Vega and Miguel[8]and Sakamaki et al.

Table 3 .
Influence of the number of molecules N in the simulation with the SPC/E model at a constant temperature 350 K on the following results: ρ L liquid density, ρ V vapor density, t 10-90 interface thickness, γ surface tension.