Tetrahedrality, hydrogen bonding and the density anomaly of the central force water model. A Monte Carlo study.

Monte Carlo computer simulations in the canonical and grand canonical statistical ensemble were used to explore the properties of the central force (CF1) water model. The intramolecular structure of the H2O molecule is well reproduced by the model. Emphasis was made on hydrogen bonding, and on the tehrahedral, q, and translational, τ, order parameters. An energetic definition of the hydrogen bond gives more consistent results for the average number of hydrogen bonds compared to the one-parameter distance criterion. At 300 K, an average value of 3.8 was obtained. The q and τ metrics were used to elucidate the water-like anomalous behaviour of the CF1 model. The structural anomalies lead to the density anomaly, with a good agreement of the model's density with the experimental ρ(T) trends. The chemical potential-density projection of the model's equation of state was explored. Vapour-liquid coexistence was observed at sufficiently low temperatures.

Introduction

Water is one of the most important substances on Earth. It is a key player in a vast variety of biological, geological, environmental, engineering, and technological processes. The majority of interpretations of the experimentally observable properties of liquid water have been possible through modelling. Various models of water exist and can be broadly classified into three groups: quantum-chemical, atomistic, and coarse-grained [1]. In 1975, Stillinger, Lemberg and Rahman proposed an isotropic coarse-grained central force (CF) model of water [2-4] and later introduced some improvements [5]. In this model, water is regarded as a weak electrolyte, where oxygen and two hydrogens spontaneously form ion-triplets, i.e., a H2O molecule. A collection of three pair potentials describing the hydrogen-hydrogen, hydrogen-oxygen, and oxygen-oxygen interactions are responsible for a non-linear geometry of H2O molecule and its dipole moment. In a revised version of the model (CF1) by Haymet et al. [6], the set of potential parameters was adjusted in such a way as to bring the pressure of the model system at 25°C closer to the atmospheric. The model has also been used to study solvation of ions [7-10] and of hydrophobic solutes [11, 12] as well as at the interface with a planar wall (electrode) [13-15].

The most appealing feature of the central force models of water lies in the fact that a system of water molecules is treated as a mixture of partially charged particles representing oxygen and hydrogen atoms in the 1 : 2 number ratio. Compared to commonly used atomistic water models (for example SPC and TIP models) [1], no angular dependent terms are included in the CF model. This is an advantage for simulations as well as for theoretical manipulations (for example, the Ornstein-Zernike equation can readily be applied to the model, cf. references [6, 16–18]). Additionally, such a model can be used to potentially address the autoprotolysis of water molecules (2 H2O ⇌ H+ + H3O+) and consequently the pH of the solution.

Various CF models have been studied mainly by means of molecular dynamics (MD) computer simulations and integral equations. To our best knowledge, Monte Carlo computer simulations were used by Bresme [19], though not for the CF1 model proposed by Haymet et al. [6]. In this work, the CF1 water model was examined by means of both canonical and grand canonical Monte Carlo simulations. The emphasis was made on the geometry of the water molecule, hydrogen bonding, the tetrahedral and translational order parameters, the density anomaly and the equation of state.

The model

In the CF1 model, the electroneutrality of H2O molecule is ensured by assigning the partial charges of qH = 0.32983e0 to hydrogens and qO = −2qH to oxygen, where e0 designates the elementary charge. The oxygen-oxygen (UOO), hydrogen-hydrogen (UHH), and oxygen-hydrogen (UOH) effective pair potentials are defined as [6]: where r is the particle-particle separation distance. Inserting the numerical value of r in Ångströms, the unit of the potential functions (2.1–2.3) is kcal/mol (throughout the paper, energetic unit will be given in kJ/mol). All three potentials are displayed in figure 1.

Figure 1.

Distance dependence of the oxygen-oxygen (OO), hydrogen-hydrogen (HH), and oxygen-hydrogen (OH) pair potentials for the CF1 water model (cf. equations (2.1–2.3)).

Monte Carlo computer simulations

The behaviour of the CF1 water model was explored using Monte Carlo computer simulations in the canonical (N, V, T) and grand canonical (μ, V, T) statistical ensemble. N represents the number of particles, V the system’s volume, T the temperature, and μ the chemical potential. All simulations were performed in a cubic simulation box with the box-edge length L, determined from the desired starting density, ρ, and number of water molecules N. The initial configuration of the system was obtained by a random insertion of particles (oxygens and hydrogens in number ratio 1:2). Periodic boundary conditions were implemented along with the Ewald summation method (Ewald parameter α was set to 5/L) [20, 21].

200 water molecules (i.e., 200 oxygens and 400 hydrogens) were used in canonical Monte Carlo (CMC) simulations while the equilibrium number of water molecules in the grand canonical Monte Carlo (GCMC) simulations was determined by the chemical potential, μ. In the CMC simulations, the equilibration part involved 5 · 108 attempted particle displacement steps, followed by 1.5 · 109 steps for the production part. In a displacement step, a particle of the system (oxygen or hydrogen) was selected at random and moved to a new randomly selected position (max. displacement of 0.5 Å). Standard Metropolis sampling algorithm was employed [20, 21], with the acceptance probability for particle displacement being min [1, exp (−β [Unew − Uold])]. β = 1/kBT (kB is the Boltzmann constant), while Unew and Uold denote the energy of the system after and before the attempted particle move, respectively. In the GCMC simulations, a 1 · 108 steps long CMC simulation involving 200 water molecules was performed before the initiation of the grand canonical algorithm. The equilibration part in total involved 5 · 108 steps (1 · 108 CMC and 4 · 108 GCMC), while production runs were 1.5 · 109 steps long. In the GCMC part, one attempted particle displacement was followed by an attempt to insert or remove a water molecule (one oxygen and two hydrogens). For the deletion of a water molecule, a random oxygen atom was chosen alongside the two closest hydrogen atoms. These three particles were removed from the system. As for insertion, an oxygen atom was added to a random position inside the simulation box, then two hydrogens were added in a favourable geometry, so that the distance between the oxygen and hydrogen was in the interval 0.875 ⩽ rOH/ Å ⩽ 1.075 and the distance between the two inserted hydrogens was 1.3 ⩽ rHH/ Å ⩽ 1.7. The acceptance probability for water insertion was min [1, Y] and for the deletion it was min , where [22, 23]: Here, and denote the numbers of oxygens and hydrogens, respectively, in the state with less waters (i = N), i.e., before insertion or after removal, and the state with more particles (i = N + 1), i.e., before removal or after insertion. Parameter B is related to the chemical potential of H2O, μ, in the following way: where is the thermal de Broglie wavelength for hydrogens (i = H) or oxygens (i = O). h is the Planck constant and mi is the mass of the particle i = H or O.

The majority of simulations were performed at 300 K and water density of 1 g/mL. The temperature dependence of the studied quantities was explored up to 1000 K (at density of 1 g/mL), while density dependence at 300 K was studied in the range from 0.8 to 1.4 g/mL.

Results and discussion

In this section we present the results of Monte Carlo simulations of the CF1 water model. We studied the molecular geometry of the CF1 water, as well as the effect of temperature and density on it. We were also interested in hydrogen bonding. Two criteria were used to determine the average number of hydrogen bonds: the simplest one-parameter distance criterion based on the integration of the OH radial distribution function and an energy criterion based on the pair energy distribution function. Alongside this, using the tetrahedral and translational order parameters we studied the tetrahedrality of the model water and the concomitant density anomaly. Using the GCMC simulations, the temperature dependence of the model water density was calculated and compared with experimental data. Furthermore, the μ − ρ projection of the equation of state for the CF1 water model was examined in a broad temperature range.

Geometry of the CF1 water molecule

In figure 2a we show the radial distribution functions (RDFs) for a CF1 water model at density 1 g/mL and temperature 300 K, while figure 2b displays the running coordination numbers estimated from the corresponding RDFs, i.e. , where ρ is the number density of species j. The first peak in the oxygen-hydrogen RDF corresponds to the intramolecular O-H bond, while the first peak in the hydrogen-hydrogen RDF corresponds to the interaction between two hydrogens of a given water molecule. The intramolecular OH and HH coordination numbers are 2.0 and 1.0, respectively (indicated by a grey line in figure 2b) in excellent agreement with the composition of the water molecule. The average O-H bond length is lOH = 0.962 Å, and the average separation distance between the two hydrogens of the H2O molecule is 1.496 Å. From these two values, the HOH bond angle is θ = 102.1°. The lOH is within the uncertainty interval of the experimentally determined O-H bond length for water in the liquid state, measured by neutron diffraction, (0.970 ± 0.005) Å [24], and approximately 0.5% larger than the experimental value of a free water molecule in the gas phase, (0.9572 ± 0.0003) Å [25]. The H-H distance is approximately 3.5% shorter than the experimentally determined, (1.55 ± 0.01) Å [24], which leads to smaller bond angle than observed in experiments, i.e., 106.1° ± 1.8° for liquid water and 104.52° ± 0.05° for water in the gas phase. From the un-normalized bond angle distribution function, , shown in figure 2c one can see that the reported bond angle of the CF1 water corresponds to the most probable value, although the distribution is quite broad and somewhat asymmetric towards larger θ. In table 1, the values of lOH and θ for some atomistic water models frequently employed in computer simulations are given. The parameters of TIPxP (x = 3, 4, 5) water models are equal to experimental values for an isolated water molecule, while ST2 and SPC models assume a perfect tetrahedral angle and O-H bond length of 1 Å.

Figure 2.

The oxygen-hydrogen (blue triangles and lines), hydrogen-hydrogen (green stars and lines) and oxygen-oxygen (red circles and lines) radial distribution functions (panel a) and the corresponding running coordination numbers (panel b) for the CF1 water model with density 1 g/mL and temperature 300 K. The height of the first peak of OH RDF is ~ 39 (not shown). Panel c shows the unnormalized distribution of HOH bond angles, . Symbols show the data obtained by canonical Monte Carlo simulations while lines are only guides to the eye.

Table 1.

Geometric parameters (O-H bond length and HOH bond angle) of some commonly used water models [1], and the experimental values for a water molecule in the gas [25] and liquid phase [24].

Water (model)	lOH, Å	θ, °
ST2	1	109.47
SPC, SPC/E	1	109.47
TIP3P, TIP4P, TIP5P	0.9572	104.52
CF1	0.962	102.1
Experiment (gas)	0.9572	104.52
Experiment (liquid)	0.970	106.1

Experiments show that lOH and θ remain unchanged within the experimental error in a broad temperature range (from 298 K to 473 K) [24]. We simulated the CF1 model at 1 g/mL in the interval from 240 K to 500 K and observed no changes in lOH and θ as well. The differences were less than 0.3%. In addition, the intramolecular geometry also remained the same at ρ = 0.9 and 1.1 g/mL in the reported temperature range.

Hydrogen bonding in CF1 water

The number of hydrogen bonds per water molecule can be approximately estimated from the OH coordination numbers determined at distances which correspond to the first and second minimum in the oxygen-hydrogen RDF (cf. figure 2a). For 1 g/mL at 300 K, integrating the gOH(r) up to the second minimum gives the coordination number of 4.14. Subtracting the value of the intramolecular OH coordination number (2.00) one obtains 2.14. The number of hydrogen bonds (donors and acceptors) between a test water molecule with its surrounding water molecules is therefore equal to 2 × 2.14 = 4.28. Experimental estimate for the number of hydrogen bonds per water molecules is approximately 3.5 [26]. Using the one-parameter distance definition to determine the average number of hydrogen bonds of the CF1 water model leads to a larger number compared to experiment. In this work, we investigated the one-parameter distance criterion because such hydrogen bond definition is often used in analysing the results of integral equations. However, in simulations involving atomistic water models, commonly extended “distance-only” criteria are used where the O-O and O-H distance together with OH-O angle are taken into account as well [26]. Various cut-offs for the hydrogen bond donor-acceptor distance and the hydrogen-donor-acceptor angle have been defined [26]. Employing such a distance-angle definition for the hydrogen bond usually leads to a better agreement with experimental data.

In figure 3, the temperature dependence of the intermolecular part of the gOH(r) (panel a) and the running OH coordination number (panel b) are shown. We see that upon increasing the temperature of the system, the second minimum in the OH RDF slightly shifts towards larger r-values. The changes are, however, small: at 300 K, the minimum was found at 2.43 Å, while at 400 and 600 K it was at 2.48 Å. Figure 3b displays the sensitivity of the intermolecular OH coordination number on the gOH(r) integration end-point, r. The running coordination numbers, nOH(r), for different temperatures cross at r = (2.45 ± 0.02) Å. Very small differences in determining the integration end-point lead to opposing temperature trends in nOH (indicated by arrows in figure 3b). In table 2, the average number of hydrogen bonds, determined by a distance criterion (r = 2.5 Å), 〈H-bond〉, is given. Basically, no difference in the disruption of the hydrogen bonding network upon increasing temperature can be observed (all values of 〈H-bond〉 are approximately 4.3). Integration of the OH RDF is therefore not a very reliable way of determining the number of hydrogen bonds in the case of CF1 water model. It is also sensitive to the bin-width used in collecting the RDFs during the MC simulations (in our case, the bin-width was 0.05 Å).

Figure 3.

The temperature dependence (300, 400, 600 K) of the intermolecular part of the oxygen-hydrogen radial distribution function (panel a), and the corresponding running OH coordination numbers (panel b). The density of the CF1 water model was 1 g/mL.

Table 2.

The temperature dependence of the average number of hydrogen bonds estimated from the OH RDF (one-parameter distance criterion, 2.5 Å), 〈H-bond〉, and from the pair energy distribution function (with H-bond energy −9.0 kJ/mol), 〈H-bond〉, of the CF1 water model with density of 1 g/mL.

T / K	300	350	400	600
〈H-bond〉r	4.28	4.32	4.35	4.33
〈H-bond〉E	3.79	3.75	3.20	2.27

We have therefore estimated the number of hydrogen bonds also from the energetic definition. In figure 4, the pair energy distribution function, , of the CF1 water model at various temperatures (300, 400, and 600 K) is given. The shows a minimum which fades-off with an increasing temperature. The minimum is located at −8.7, −9.4, and −11.3 kJ/mol for 300, 400, and 600 K, respectively. The values are somewhat less negative than for the atomistic water models. For example, TIP3P water model at 300 K has the pair energy minimum at approximately −10.5 kJ/mol [27]. The minimum in the is a convenient energetic definition of a hydrogen bond [26]. We have rather arbitrarily chosen that any pair of CF1 waters having energy at least −9.0 kJ/mol is considered to be hydrogen bonded. Table 2 gives the temperature dependence of the average number of hydrogen bonds estimated from such an energetic criterion, 〈H-bond〉. The average number of hydrogen bonds at 300 K, determined in this way, is 3.79. This is much closer to the experimentally determined value for liquid water (~ 3.5 [26]) than the value estimated from the gOH(r). The increase of temperature also leads to the correct trend, namely the hydrogen-bonding network weakens with increasing temperature (average number of hydrogen bonds decreases with temperature). From the results presented in this subsection, we conclude that in case of the CF1 water model, the pair energy criterion of the average number of hydrogen bonds gives more consistent results compared to the distance criterion.

Figure 4.

The unnormalized pair energy distribution function, , of the CF1 water model at 1 g/mL and 300, 400 and 600 K.

Tetrahedrality of the CF1 water

Water is a tetrahedral liquid [1, 28, 29].The first nearest neighbour water molecules correspond to the first peak in the gOO(r). The second peak in the oxygen-oxygen RDF, however, is argued to indicate the presence of tetrahedral order of waters in the liquid state [6, 30]. Figure 5 shows how the gOO(r) changes with temperature. By increasing the temperature the height of the first and second peaks diminishes. The location of the second maximum slightly shifts towards larger distances (indicated by the arrow) and becomes only marginally pronounced at high temperatures (600 K). This indicates that the tetrahedral order of CF1 water diminishes with an increasing temperature. As also seen from the average number of hydrogen bonds, 〈H-bond〉 (table 2), the hydrogen bonding network in the CF1 water model weakens with temperature.

Figure 5.

The temperature dependence of the oxygen-oxygen pair distribution function, gOO(r), of the CF1 water model at 1 g/mL.

A variety of different order parameters have been proposed to characterise the local structure of water in the liquid state. Here, we explore the temperature and density dependence of the so-called angular tetrahedral order parameter, defined as [31, 32]: where ψ denotes the angle between the vectors joining the oxygen atom of a selected test water and its nearest neighbour oxygen atoms i and j. It describes the angular ordering of the first hydration shell waters. For an ideal gas, the average value of q equals zero, while for a perfect tetrahedral arrangement it equals one. High values of q can be found for an angularly ordered but radially disordered hydration shell as well as in cases where there is not a clear separation between the first and second shell. By cooling the system, it becomes more tetrahedral, although it still differs from the ice-like configuration.

In figure 6a, we show the distribution of the tetrahedral order parameter, , of the CF1 water model at 1 g/mL for four different temperatures (240, 300, 400 and 600 K). At low temperatures (240 and 300 K), a shoulder-like distribution is observed. Such a shape is characteristic of cases where compatible fractions of two distinct local environments exist: a tetrahedral and nontetrahedral [32]. The most probable value of q [indicated by the maximum on the ] is 0.81 at 240 K and moves to 0.77 by increasing the temperature to 300 K. A bi-modal distribution observed at 240 and 300 K does not mean that in the liquid there exist two types of arrangements that would be energetically favoured. It only means that the arrangement of four nearest neighbours of the central water molecule can be predominantly ice-like structured (the peak at high q values) or unstructured (smaller peak at lower q values) [32, 33]. Increasing the temperature even further (400 and 600 K) leads to the loss of strong tetrahedrality and at 600 K a completely uni-modal distribution is observed with the most probable q value of 0.49. In figure 6b, the temperature dependence of the average tetrahedral order parameter, 〈q〉, is shown. The 〈q〉 monotonously decreases with an increasing temperature, in line with temperature trends in gOO(r). It is worth noting that at 240 K the system is still in the liquid state rather than a solid ice.

Figure 6.

Distribution of the tetrahedral order parameter, , for the CF1 water model at 1 g/mL for four different temperatures (240, 300, 400 and 600 K) (panel a) and the temperature dependence of the average tetrahedral order parameter, 〈q〉, at 1 g/mL (panel b).

Density anomaly of the CF1 water model

For all liquids, the local order is determined by the pair correlations. To evaluate the degree of ordering due to pair correlations between oxygen atoms, the so-called translational order parameter was calculated as [34, 35]: where gOO(r) denotes the oxygen-oxygen pair correlation function, and rc is an appropriate cut-off distance. In our case, rc was chosen as a half of the simulation box edge-length. For an ideal gas, τ is zero. τ increases as the system becomes more ordered. The degree of correlation between the average tetrahedral order parameter, 〈q〉, and the average translational order parameter, 〈τ〉, can be used to explore the structural anomalous region in the density-temperature plane. This region approximately encloses the region of the thermodynamic (density) anomaly [32].

We have performed canonical Monte Carlo simulations of the CF1 water model at 300 K for different densities, ranging from 0.8 to 1.4 g/mL. The dependence of 〈q〉 and 〈τ〉 on the density is given in figures 7a and 7b, respectively. The shape of 〈q〉 on ρ is characteristic of tetrahedral liquids: it exhibits a maximum at densities approaching the density of the ice-like structures. In our case, the maximum is located at 0.95 g/mL. The translational order parameter has a shallow maximum at low densities followed by a deeper minimum at 1.15 g/mL. The maximum in 〈τ〉 at 0.95 g/mL coincides with the maximum in 〈q〉. The locus of maxima and minima in the density dependence of 〈τ〉 is characteristic of the structurally anomalous region [36]. In addition, such a region also implies a linear correlation between the 〈q〉 and 〈τ〉 [32, 36], as seen for the CF1 water model in figure 7c (the anomalous region is marked with a grey rectangle).

Figure 7.

Density dependence of (a) average tetrahedral order parameter, 〈q〉, and (b) average translational order parameter, 〈τ〉. Correlation of 〈q〉 with 〈τ〉 is given in panel (c). Arrows indicate increasing density (from 0.8 to 1.4 g/mL). Data apply for T = 300 K.

Next, we explore how the CF1 water model captures the correlation between the local structural order and density anomaly. It is well known that liquid water at ambient pressure exhibits a density anomaly, i.e., the density reaches a maximum value at around 277 K [1]. The origin of the anomalous behaviour lies in increased structural fluctuations upon cooling water down to the Widom line [37]. We used GCMC simulations to explore the temperature dependence of CF1 model density in the range from 273 to 373 K. A linear approximation for the dependence of the chemical potential, μ, on temperature [38] was assumed in the whole temperature range. In figure 8, we compare our calculations with the experimental data for liquid water. Calculated density at a given temperature, T, corresponds to the following chemical potential: μ = −589.0 kJ/mol − 0.08 kJ/mol K · T. We see that a linear approximation for μ(T) gives good agreement between the calculated and experimental density of water. The difference between the calculated and experimental density is only 0.25% at 277 K, 0.65% at 298 K, and around 1.3% at 373 K. These differences are smaller than the uncertainties in Monte Carlo calculations, which are of the order of 1.5–2%. It is also worth mentioning that in the entire temperature range studied, the relative differences in experimental densities are small. For example, experimental densities for water at 277 and 373 K differ by approximately only 4.2%. The CF1 water model also exhibits a small density maximum that can be viewed as the density anomaly (true density anomaly, however, should be explored at constant pressure conditions). Frequently used atomistic models of water (for example SPC or TIP3P) do not give such good results (data for TIP3P water model are given in figure 8).

Figure 8.

Temperature dependence of the water density. Displayed are experimental data for liquid water (blue x’s), CF1 water model (red circles) and TIP3P water model (green stars). For the CF1 model, the results were obtained by GCMC simulations using a linear dependence of the chemical potential on temperature. Data for TIP3P and experiment was taken from [39].

Preliminary insights into the chemical potential-density projection of the equation of state for the CF1 water model

Here, we also present some preliminary insights into the behaviour of the μ − ρ projection of the equation of state of the CF1 water model in a limited temperature interval. A more extensive investigation of the model’s equation of state — namely exploration of the phase transition, determination of the coexisting densities, construction of the pressure-temperature projection of the phase diagram — calls for a separate study.

Employing the GCMC simulations, we determined the dependence of the system’s density, ρ, on the chemical potential, μ, at different values of temperature (ranging from 300 to 650 K). Results are shown in figure 9a: at sufficiently high temperatures (above 400 K), the dependence of the ρ on μ is monotonous, i.e., the density increases by increasing the chemical potential. The ρ(μ) is a single-valued sigmoid-like function (a given value of the chemical potential determines a unique density of the system). Such a behaviour is characteristic of a system above the vapour-liquid critical point. By decreasing the temperature, the change in density (from gas-like to fluid-like) becomes more steep, until two distinct branches occur: one corresponds to the gas phase and the other to the liquid phase. For a given value of the chemical potential, two densities (low and high) can be found, and this value of μ corresponds to the liquid-vapour coexistence. Hysteresis loop upon approaching the phase transition from the gas phase densities is observed.

Figure 9.

(a) The μ − ρ projection of the equation of state of the CF1 water model, and (b) the density fluctuations, 〈ρ2〉 − 〈ρ〉2, at various temperatures.

The fluctuations in the system’s density, 〈ρ2〉 − 〈ρ〉2, as a function of the chemical potential, μ, are shown for different temperatures in figure 9b. We can see that by decreasing the temperature (from 650 down to 300 K), the function starts to exhibit sharper peaks. The increasing fluctuations in density are an indication that by further decreasing the temperature, the system will undergo the phase transition (in this case the vapour-liquid phase transition). The chemical potential at which the function diverges corresponds to the phase transition. At this point, we need to stress that our system was rather small, composed of approximately only 200 water molecules. This fact prohibits one to adequately address particle fluctuations.

The standard GCMC simulations employed in this work along with the small size of the system do not allow for a precise determination of the critical point. Further temperature examination of the ρ(μ) dependence for a system with a larger number of water molecules is needed as well as the use of more sophisticated sampling techniques (for example [40]) to construct the density-temperature projection of the phase diagram. The only report on the phase diagram of a modified version of the central force model used in this work has been reported in [41]. The author used nonequilibrium molecular dynamics simulation to calculate the vapour-liquid coexistence by explicitly simulating the vapour-liquid interface. The initial configuration of 540 water molecules arranged in a face-centered cubic lattice was subjected to a heat flux along one of the simulation box directions, leading to the generation of the liquid and vapour phases. To obtain the equilibrium phases, the temperature gradient was removed and constant temperature along the box was maintained. From local densities in different layers, the coexisting densities were determined. Good agreement of the vapour-liquid coexistence curve with experimental data for water was observed. Scaling law with the Ising exponent was used to obtain the critical temperature (Tc = 630.4 K) and critical density (ρc = 0.337 g/mL) of the model, and critical pressure (pc = 337 bar) was determined by using the Clausius-Clapeyron equation. From our calculations, we can only conclude that the critical temperature of the employed CF1 water model is quite below the experimental value for water (647.1 K). The reasons need to be addressed in the future: it can be due to a small size of the system used in our simulations, due to the used GCMC methodology, or perhaps due to differences in the model (cf. figure 1 of [41] for comparison of potentials used in this work and by Bresme).

Conclusions

We have studied the CF1 model of water using the canonical and grand canonical Monte Carlo computer simulations. We have shown that the CF1 model reproduces well the molecular geometry of the water molecule. Our focus was on the hydrogen bonding of the model water molecules and on the temperature trends in the average number of hydrogen bonds. By using the energetic criterion (i.e., the minimum in the pair energy distribution function) instead of the commonly used one-parameter distance criterion (integration of the OH radial distribution function), we obtained a better agreement of the average number of hydrogen bonds with experiment as well as correct temperature trends. We also studied the evolution of the tetrahedral and translational order parameters as a function of water density. Structural anomalies were observed at 300 K leading to correct temperature trends in the density of water (so-called density anomaly). The μ − ρ projection of the equation of state allowed us to address the vapour-liquid phase transition.

In the future, a more detailed analysis of hydrogen bonding in the CF1 water model is needed. This includes an extension of the one-parameter distance criterion to distance-angle criteria, which are often used in explicit water computer simulations. In addition, knowledge of the μ − ρ projection of the equation of state in a wide range of phase space is extremely important, for example, to study adsorption phenomena. Our results of the equation of state of the model are only preliminary. Simulations of systems with a much larger number of CF1 water molecules will be required, together with the use of more sophisticated simulation techniques. This will allow the determination of the coexisting densities and the construction of the p − T phase diagram. Comparison with the experimental phase diagram of water will provide opportunities to modify the model’s potential functions to achieve a better agreement of the model with real water.