Molecular simulation of water vapor – liquid phase interfaces using TIP 4 P / 2005 model

Molecular dynamics simulations for water were run using the TIP4P/2005 model for temperatures ranging from 250 K to 600 K. The density profile, the surface tension and the thickness of the phase interface were calculated as preliminary results. The surface tension values matched nicely with the IAPWS correlation over wide range of temperatures. As a partial result, DL POLY Classis was successfully used for tests of the new computing cluster in our institute.


Introduction
Water is perhaps the most studied substance in the world due to its importance in daily life, industry or physical, chemical or biological processes.Due to its many anomalies and non-standard behavior it is however very hard to model.The non-trivial phenomena are caused by its polar character and consequently by its formation of hydrogen bonds.
As for numerical simulations, water has been studied extensively over past decades with the pioneering works by Barker and Watts [1] and Rahman and Stillinger [2].Despite its compexity, simple models of water are widely used, because they work surprisingly well.This moreover allows computation of more extensive systems consisting of thousands of molecules.Many simple models are developed, but as the best of them are considered SPC/E model and TIP4P/2005 model [3,4].
Motivation of our work is to shed light on the discrepancies between experiments and simulations, e.g. the second inflection point of water [5], to reproduce measurements of the surface tension of supercooled water as well as to enhance our theoretical work concerning nucleation rates predictions [6] and capillary waves modeling [7].
In this paper we publish preliminary results of part of our newly formed simulation group.Primary objective was to get the software executing the simulation procedure working in our new cluster in the direction of our interests.We used DL POLY Classic program on 4 × 24 Intel(R) Xeon(R) CPU E5645 @ 2.40GHz CPUs, debugging and early computations were executed on 4 AMD Phenom (tm) 9600 Quad-Core GPUs.
We simulated a water slab for various temperatures from 250 K towards the critical point (which was not exceeded) stopping at T = 600 K. To have a relevant reference point, we followed specifications of papers [8] and [9].a e-mail: barbora.plankova@gmail.com

TIP4P/2005 model
Although water is flexible and polarizable, simple rigid non-polarizable models are used in most of simulations.The reason is that some of them represent experimental data very well, which allows one to model larger or longerrunning systems.All the simple models from SPC or TIPXP families have in common that they locate the positive charge on the hydrogen atoms and a Lennard-Jones interaction site on the position of the oxygen.They however differ in bond geometry (angle of H-O-H bonds), charge distribution or their target properties.Three-site models TIP3P [10] and SPC [11] were developed in 1980s.An SPC modification, SPC/E model [12], is very much used nowadays.SPC and SPC/E models differ in a way that SPC/E includes some empirical corrections to a polarization energy.This model became very successful and reproduces data very well.
To improve the electrostatic distribution within the water molecule, an auxilliary atom M was introduced in TIP4P [10] model.This atom carries the negative charge, is massless and close to oxygen atom (0.15 Å).Probably the best model is TIP4P/2005, proposed by Vega and Abascal [13] in 2005.They tried to combine good phase diagram of TIP4P with target properies of SPC/E improving the melting point.It has been shown [3,4] that TIP4P/2005 behavior is even better than that of SPC/E.Therefore, in this work we use the TIP4P/2005 model.Another extension TIP5P was proposed by Mahoney et al. [14] to carry the negative charge on two auxiliary atoms.However, the performance is not better than of TIP4P/2005.TIP4P/2005 model includes Lennard-Jones interactions between oxygen atoms only; hydrogens have negligible mass compared to them, which makes the simulation easier.Other interaction is electrostatic which occurs between hydrogen H and charge M atoms.TIP4P/2005 parameters are listed in table 1.

Simulation methods
On our new cluster we installed the DL POLY Classic (former DL POLY 2) and DL POLY 4 programs.We performed computations using both, but then we decided to use DL POLY Classic on this problem.The reason was that DL POLY 4 is designed for larger systems and its usage was less efficient.
We performed a number of molecular dynamics simulations within NVT (canonical) ensemble, i.e. number of particles, volume and temperatures were fixed.For an initial configuration of molecules, we used phase centered cubic mesh where four molecules are put to a corners of one cube, these cubes are put next to each other with a blank space of a size of one cube.Therefore, number of molecules N is restricted to numbers 4 × N 3 c , where N c is a positive integer (see table 2).We chose the number of molecules of 1372 which is quite close to 1024 molecules considered by Sakamaki et al. [8].
The simulation was performed as follows: first, a liquid cubic box was run for 50 ps, then the z-size of the box was expanded to approximately 3× the original proportion and run for 10 ns to provide reliable data for surface tension determination.Sizes of the box were calculated depending on the NIST [15] values of the water saturated liquid density for particular temperature.For supercooled region, a constant box size corresponding to 300 K system was used.Periodic boundary conditions were used in all directions.This imposed a small problem in extension of the system in z-direction -molecule can be on boundaries of the simulation box, therefore, simple extension of the axis would "break" the molecule.Between the two simulations, a Fortran code was run to rotate the molecules so that none of them is on the boundary of the z-direction.Timestep of the simulation was chosen as 2 fs (same as in [8]) with velocity Verlet integrator.To maintain constant temperature, we used the Nosé-Hoover thermostat with relaxation constant 100 fs.Cutoffs of Lennard-Jones interactions and van der Waals forces were set to 14.5 Å, close to the value in [8].For electrostatic interactions, direct Ewald method was used, with automatic parameter optimization constant set to 10 −5 .This value gave a convergence parameter 0.192 Å −1 and k-vector indeces 6 × 6 × 6 for a cubic box, 6 × 6 × 18 for slab geometry.For higher temperatures, the values got up to 7 × 7 × 7 and 7 × 7 × 29 respectively.Density was computed as a histogram in z-direction in every step and averaging through the time.Examples of the density profiles converted to g/cm 3 for 300 and 500 K can be seen in figure 1 depicted by the solid lines.Snapshots of the simulations for two temperatures are shown in figure 2.
The density profiles were divided to two halves approximately in the centre of the simulation box.One half of the density profile was subsequently fitted to a hyperbolic tan-  gent density profile model, where ρ L and ρ V are fitted bulk densities, z 0 is the position of the Gibbs dividing surface of the interface, d is the parameter for the thickness of the interface.For z and ρ we used the simulation data to be fitted.The results of the fit are depicted as the dashed lines in figure 1.The fitted values were used by evaluating the surface tension γ in the following manner: In equation ( 2), L z denotes the box size in z direction, P ii is the ii-th component of the pressure tensor, and σ are the Lennard-Jones parameters for oxygen atom, and r c is the cutoff for the Lennard-Jones potential.Their values as well as other model parameters are summarized in table 1. Second term in equation ( 2) corresponds to the Lennard-Jones tail correction [16].

Results
Insipired by the work of Sakamaki et al. [8], we performed molecular dynamics simulations of a water slab enclosed by the vapor for temperatures T = 250 K, 270 K, 275 K, 300 K, 350 K, 400 K, 450 K, 500 K, 550 K, and 600 K. Dimensions of the box varied from 34.53 Å corresponding to a density of liquid water at 300 K which was used also for systems at temperatures below 300 K up to 39.83 Å for 600 K.We decided to vary the simulation box sizes, since it corresponded to reality better than box with constant sizes.Moreover for 600 K, small box forced the almost incompressible liquid water to too low volume which caused an "explosion" of the liquid cube, shifted the slab and made the system unnecessarily non-equilibrium at start of the slab simulation.
Given by the parameters stated in previous section, computations ran approximately 3 days if running on all 24 CPUs.Melting temperature for TIP4P/2005 is T m = 249 K, therefore even for lowest temperature of 250 K we did not encounter any of the liquid water during the simulation.
Figure 1 shows density ρ varying with the z coordinate averaged through the simulation.The structure of the simulation box can be seen: there is vapor on the left and right side of the z axis, while the liquid phase is allocated in the centre of the box.For the purposes of parameters fitting, only the right half of the profile was necessary.As the temperature increased, the bulk liquid density decreased and bulk vapor density increased.Figure 2 shows an example of two configurations for two temperatures (300 K, 550 K).The liquid phase persisted in the centre of the simulation box, while the vapor phase gradually expanded into the vacuum space after the box got stretched in the z-direction.As can be seen at low temperatures, the molecule escaping from the liquid phase into the vapor phase was rather rare event.On the other hand at the elevated temperatures, the vapor phase got significantly denser.
Surface tension computed using equation ( 2) is shown in figure 3 as circles, the IAPWS values [17,15] are shown as a solid line.Simulated values nicely describe the reference data.
Figure 4 shows the thickness of the profile t, known also as the "90-10" thickness defined as t = 2.1972d, where d is the thickness parameter obtained from equation ( 1).An exponential trend can be observed which agrees with the temperature scaling along the critical point (critical scaling theory) [18,19].

Conclusions
In this work we performed molecular dynamics simulations for water for various temperatures ranging from 250 K to 600 K.As a water model we used TIP4P/2005 which is probably the best rigid non-polarizable water model available at the moment.We computed surface tension as well as the interface thickness using the hyperbolic tangent fit as preliminary results.
Main goal of this paper was to test our new computing cluster on molecular simulations and to investigate the performance of the DL POLY Classic software.Surface tension data for water fitted nicely to the IAPWS tabulated values [17,15].In future, we would like to perform more simulations in the supercooled region of liquid water using the TIP4P/2005 and the SPC/E models to compare simulated results with our recent experiments.[5].
Also, we would like to model the so-called capillary waves contribution to the surface tension, i.e. to simulate, how the thermal motion of molecules affects (lowers) the surface tension for planar phase interface.New molecular simulations would support our theoretical work [6,7,[20][21][22].

Fig. 1 .
Fig. 1.Density profiles ρ as functions of z-coordinate for systems having temperatures T = 300 K and 500 K. Solid lines are timeaveraged profiles obtained from the simulations, dashed lines are fits of the right-hand profiles (liquid -vapor) to the hyperbolic tangent profile, equation (1).

Fig. 2 .
Fig. 2. Snapshots of two configurations during the simulations for two temperatures, a) T = 300 K and b) 550 K. Liquid cubes are in the middle, vapor is on the left and right side (not so apparent for lower temperature).

Fig. 4 .
Fig. 4. Thickness of the profile t as a function of temperature T predicted in this work using the tangent hyperbolic fit of parameter d and reevaluation t = 2.1972d.

Table 1 .
Simulation parameters of TIP4P/2005 water molecule atoms: oxygen O, hydrogen H, mass-less charge M. m is molar mass, q charge in units of elementary charge, and σ are Lennard-Jones parameters.

Table 2 .
Number of molecules N permitted by the phase centered cubic mesh algorithm of the water molecules initial configuration.N c is arbitrary positive integer.