Free-energy differences of OPC-water and SPC/HW-heavy-water models using the Bennett acceptance ratio.

Surface tension, vapor density of OPC-water and SPC/HW-heavy-water models have been estimated at low temperatures using the scaled model. The free-energy difference, - Δ F , of n-molecules and (n-1)-molecules plus a free probe has been calculated using the Bennett acceptance ratio with the aid of Monte-Carlo simulations. Our results show that the relation between the free-energy difference divided by k B T and the number of molecules to the power minus one-third is linear for n > 6 . Consequently, the surface tension can be extracted from the straight line slope, whereas the vapor density can be extracted from the intercept, which is proportional to the logarithmic ratio of liquid density to that of vapor density. By scaling the free-energy differences, for at least three different temperatures, to T C T - 1 , we estimated the critical temperature and hence the surface tension and the vapor density at a wide range of temperatures. The free-energy differences have been calculated at 240K, 260K, and 280K for OPC-water, and at 260K, 280K, and 300K for the SPC/HW-heavy water model.

Introduction

Water is the most important substance on earth due to its use and applications in our daily life [1, 2, 3]. Currently, more than 45 potential models have been proposed to study and produce all different properties of real water. The first proposed model started from the pioneering work of Bernal and Flower [4]. In the last few years, two promising proposed models were presented to estimate most of the thermodynamic properties of water and heavy water, namely, the OPC-water [5] and the SPC/HW-heavy-water [6] model. This paper aims to estimate the vapor density, the surface tension at relatively low temperatures, the excess surface entropy/k per molecule, and to propose a method to estimate the critical temperature. The standard method to estimate the vapor density and the surface tension using molecular dynamic simulation (MD) and Monte-Carlo simulations (MC) is to place the molecules in a slab placed between two vacuum boxes. The density profile is then either fitted to a hyperbolic tangent function or to an error function by guessing the vapor and liquid densities [7, 8]. The problem with this method lies in estimating the vapor density at low temperatures since the number of molecules in the vacuum is almost null, for this reason, one assigns a value of zero to the vapor density. To estimate the surface tension, many scientists have applied different methods to estimate it and compare their estimations with the experimental published data [9]. The standard method of estimating the surface tension is through evaluating the components of the pressure tensor. To get a better estimation for the surface tension, one might consider the tail correction [10, 11]. On the other hand, estimating the critical temperature might be a challenge. The standard method is to calculate the vapor-liquid phase diagram as a function of temperature, then make a fit to the points using Wegner expansion with constants depending on the substance [12, 13]. Other methods are briefly mentioned in ref [14].

In this work, we estimated the values of surface tension, , the vapor density, , of OPC-water and SPC/HW-heavy-water models at three different temperatures of each model using the scaled model [15]. We attempt to calculate, using this model, the energy difference between two systems each of which contains n interacting molecules except for one which has one of its molecules replaced by a free probe. These systems are composed of clusters with small number of molecules (up to 100 molecules at most) placed inside a sphere with volume up to five times the volume of molecules. The standard method of calculating the vapor density and the surface tension requires a high number of molecules to form a bulk, while applying the scaled model, the cluster can be as small as of two molecules. The free-energy differences are scaled to where is the critical temperature, which is treated as a variable that makes all the calculated free-energy differences at different temperatures collapse to a single line when it is plotted versus the number of molecules to the power minus one-third. By taking as a variable in the scaled free-energy differences, one might estimate the critical temperature as will be explained in the section of results and discussion. The surface tension then can be estimated at low temperatures by calculating the slope of the line, while the vapor density is extracted from the intercept of the line with the free energy. The intercept is related to the logarithmic ratio of , where is the liquid density, which is estimated by fitting the density profile to hyper tangent function using molecular dynamics.

Theory

Statistical mechanics provides a mechanism that connects the macroscopic properties of a substance to the microscopic states in the substance. The connection is through the bridge equation, which relates the Helmholtz free energy to the partition function, where is the Helmholtz free energy, is the Boltzmann constant, is the temperature, and is the partition function. Since calculating the free energy is quite impossible, one might calculate the difference in free energy. Many attempts have been made to calculate the free-energy difference of states, such as the exponential averaging (EXP) [16], the thermodynamic integration (TI) [17], the umbrella integration (UI) [18], the umbrella sampling (US) [19], and the weighted histogram analysis method (WHAM) [20, 21]. Shirts et al pointed out that most of these methods lack a standard test, and what is called the Bennett acceptance ratio (BAR) [22] is more efficient and rigorous than other methods. Many researchers used this methods. In this method, the difference in free energy between two states iswhere , is an arbitrary number, and is the Fermi function. refers to the system which consists of -molecules plus the free probe, and refers to the system which consists of -molecules.The optimal value of has proven that when the averaging Fermi functions' ratio equals one, leaves us with . In our simulations, we used the same number of independent configurations, and in this case the variance in free energy is

Following the scheme of Hale [23], the total interaction of ensemble is , and that of ensemble is , where is the interaction energy of all the molecules, is the interaction of the probe with molecule cluster, the value of is taken to be which is very small following the analysis of Hale [23]. Also, Hale showed that the difficulty of connecting the optimal value to the free energy comes from the simulation volume of the free probe. Hale assumed that the volume that contains n-molecules is given by , where is the number of molecules in the cluster, and is a constant that should be kept constant in all simulations as proposed by Lee [24]. In our simulations, we varied to be 5, 6, and 7 for small number of molecules, and all the calculated free energies were similar, so we fixed to be 5 in all our simulations. Following Hale, the difference in free energy between ensembles and is given bywhere is the free-energy difference divided by [25]. As explained by Hale and DiMattio [26], the free-energy difference can be written aswhere , and is related to the surface tension as . If we plot versus , then represents the extrapolation at which the value of could be estimated by knowing the value of , and the slope would be at which the surface tension can be estimated.

Molecular model and simulation details

Over the decades, scientists have tried to find a single model that describes and produces all the thermodynamic properties of water. Up to this moment, there are more than 45 different models of water, such as TIP4P [27], SPC/E [28, 29], TIP3P [30], and TIP4P2005 [31]. None of these known models produce all thermodynamics properties; some are good in calculating the critical temperature; others are good in calculating the phase diagram; some are good in calculating the surface tension. The difference between these models lies in describing the water molecule as 3-point charge and 4-sites, 3-point charge and 3-sites, 3-point charge and 5-sites, and so on, with different segment diameter and energy depth, charges, angles, and bond length.

In this section, we will be using a new promising model called the OPC-water model [5, 32] (Optimal-Point-Charge model) and SPC/HW-heavy-water model [6], and all their parameters will be given in the following section.

One of the essential methods to study many-body problems in the liquid state is the Monte-Carlo simulation or the molecular dynamic simulations presented in simulation details.

Molecular interaction potential function

The OPC-water model is a four-site model, the oxygen atom carries an amount of negative charge equals to -1.3582e located at a distance from the atom equals to 0.1594Å, and each of the hydrogen atoms carries an amount of positive charge equals to +0.6791 located at the atoms. Regarding the Lennard-Jones interaction, the well's energy depth is , and the segment diameter is for the oxygen atom only. The hydrogen-oxygen-hydrogen angle equals to , and the bond length between hydrogen-oxygen is |OH| = 0.8724Å, as shown in Figure 1.

Figure 1

The OPC-water model structure.

The SPC/HW-heavy-water model is a three-site model with one Lennard-Jones interaction site on the oxygen atom with energy depth , and segment diameter . In this model, the values of the charges of atoms are as follows: for oxygen atom, it has a negative charge equals to -0.87e that is located at the atom, and for the hydrogen atoms each one has a positive charge equals to +0.435e that is also located at the hydrogen atoms. The angle between hydrogen-oxygen-hydrogen is , and the bond length between hydrogen and oxygen equals |OH| = 0.1nm, as in Figure 2.

Figure 2

The SPC/HW-heavy water model structure.

The total pair intermolecular potential between atoms i and j is the sum of Lennard-Jones (LJ) and Coulomb potentialswhere and are the Lennard-Jones (LJ) parameters for interaction between atoms of different types. The Lorentz-Berthelot rules are used for the interaction between the unlike atoms: /2 and .

Simulation details

Using the BAR method, we applied the Monte-Carlo simulation to calculate the free-energy differences [33]. We wrote a code for the four-site OPC-water model and the three-site SPC/HW-heavy-water model. The free-energy differences of OPC-water were calculated for clusters consist of 2–15, 18, 20, 25, 30, 35, 50, 75 and 100 molecules at different temperatures: T = 240K, 260K, and 280K. The free-energy differences for the SPC/HW-heavy-water were calculated for clusters consisting of 2, 3, 4, 5, 6, 7, 8, 10, 12, 14, 18, 20, 50 and 75 molecules at different temperatures: T = 260K, 280K, and 300K by using Metropolis Monte-Carlo algorithm. Regarding the simulation volume, number of equilibration and number of Monte-Carlo steps, as well as the step move in three dimensions are the same as in the work of one of the co-authors, ref [14].

Also, we applied molecular dynamics simulation using GROMACS package [34] to calculate the surface tension from low temperature up to a point near the critical temperature, from which we estimated the critical temperature. Also, we have calculated the liquid density of the OPC-water and SPC/HW-heavy-water models at all temperatures between 200K and 300K with a step size of 10K. This value is very important to calculate vapor density. The intercept of the free-energy differences versus the number of molecules raised to the power minus one-third is proportional to the logarithmic function of the liquid density ratio to that of the vapor density. Each system is equilibrated for 1ns using NPT ensemble. For NPT we used Berendsen algorithm [35] for which the temperature and the pressure is kept at 1 atm, and the time step used is 1fs. After applying NPT for 1ns, the system is then equilibrated for another 1ns using NVT ensemble with Nose-Hoover thermostat [36]. The data then are collected for another run of 4ns. In our simulations, periodic boundary conditions are used in all dimensions with cutoff radius of , and the size of our system is 1000 molecules.

Results and discussion

The whole idea behind scaling the free-energy difference between ensemble A and ensemble B to is to see if all the free-energy differences lie on the same line. If that is the case, one can estimate many thermodynamic properties such as the surface tension, the excess surface entropy/k per molecule, and the vapor density at any temperature. To ensure that the scaling method works fine, we need to simulate the studied system at least for three different temperatures. One might estimate the critical temperature from the free-energy differences by scaling them to and changing the value until all the lines coincide with each other. This method usually overestimates the true experimental value [37, 38]. Another method is to apply molecular dynamic simulations to estimate the critical temperature of the system. This method is more accurate since we are not studying real systems but models such as the OPC-water model.

In our work, we use the molecular dynamic simulations to estimate the critical temperature and the liquid densities of the OPC-water model. To ensure we have liquid and vapor phases, we insert 1000 molecules of the OPC-water inside a box surrounded by a vacuum. To get the critical temperature, which is the most important ingredient in this work, we varied the temperature from 200K up to 660K, then fitted the surface tension points at high temperature to a line and extract the value of the critical temperature at which the surface tension is zero, and we get a value of 693.3K for the critical temperature. Regarding the liquid densities at low temperatures, we fitted the density profiles to a hyperbolic tangent function as [39].where and are the bulk densities of the liquid and vapor, respectively, is the place of the Gibbs dividing surface, and d is the width of the interface. In the fitting process, we set to zero and extract the value of .

The surface tension using molecular dynamic simulations has been calculated using the following relation [40]:where is the component of the pressure tensor, is the box length along direction, and are the Lennard-Jones parameters, and is the cutoff radius.

Using Eq. (4) after scaling it to , one can estimate the surface tension. Also, since the intercept is related to the system's density in either phase, liquid, or vapor densities, one can estimate the vapor density at any temperature by knowing the value of liquid density at that temperature. The liquid density has been estimated using Eq. (6). The vapor density is estimated from the intercept of the free-energy differences, which is related to .

OPC-water model

Figure 3 shows the surface tension as a function of temperature using molecular dynamics. Since the surface tension vanishes at the critical temperature, we fitted the surface tension data at high temperatures to a line. Using this method, the intercept of the line with the temperature axis is just the critical temperature. The estimated critical temperature using the MD simulations is 692.233K compared to the experimental value of 647.096K.

Figure 3

The surface tension results of the OPC-water model using the Molecular Dynamics (MD), versus temperatures T.

Table 1 below shows the free-energy differences with a different number of molecules at different temperatures of 240K, 260K, and 280K. We plotted the free-energy differences versus the number of molecules to the power minus one-third . Figure 4(a) shows the results of the free-energy differences at three different temperatures, and Figure 4(b) shows the same results after scaling the free-energy differences to . We noticed that all the free-energy differences curves scale nicely to one curve and all the free energy calculations lie within error of 1–3%. We also noticed that the new curve becomes a straight line when the number of molecules inside the cluster exceeds six molecules. We can estimate the surface tension, the excess surface entropy/k per molecule, and vapor density at any temperature by calculating the slope and the intercept of the line from the straight line. Table 2 shows these important values. The table also shows these values as well as the excess surface entropy/k per molecule at different critical temperatures. Mistakenly researchers used the experimental critical temperature when scaling the free-energy differences [23, 25, 26], while one has to use the critical temperature of the model itself. The table also shows the estimated critical temperature using our proposed method. Our results of calculating the excess surface entropy/k per molecule is close to the experimental value of 1.5 [41].

Table 1

The free-energy differences of the OPC-water model at 240K, 260K, 280K.

n	−δF
	T = 240K	T = 260K	T = 280K
2	7.31	6.41	5.61
3	12.31	10.41	9.21
4	13.61	13.21	11.81
5	12.51	11.21	10.01
6	12.71	11.41	10.01
7	13.91	12.31	10.51
8	14.31	12.61	11.41
9	14.51	12.81	11.71
10	14.8	13.01	11.51
11	14.91	13.21	11.91
12	15.11	13.21	12.11
13	15.61	13.31	12.21
14	15.41	13.71	12.01
15	15.91	13.31	12.11
18	15.81	14.11	12.51
20	15.61	13.81	12.51
25	16.01	14.61	12.21
30	16.61	14.61	12.81
35	16.21	14.61	13.21
50	17.11	15.01	13.21
75	16.91	14.61	13.61
100	16.61	16.41	14.41

Figure 4

(a) The free-energy differences versus the number of molecules , at 240K, 260K, and 280K. (b) The free-energy differences scaled to , versus at 240K, 260K, and 280K

Table 2

The values of the slope and the intercept of the scaled free-energy difference of OPC-water model to different estimated critical temperatures.

Method name	T_C	Slope	Intercept	Ω
TC:MD	692.23	-5.81321	10.58414	1.80311363
TC:MC	765	-4.97243	9.056296	1.54232355
TC: experiment	647.096	-6.49618	11.82449	2.01495268

Figure 5 shows the surface tension values of our Monte-Carlo results compared to that of molecular dynamic simulations and experimental results. We noticed from the figure that the MC results overestimated surface tension values compared to MD results regardless of the critical temperature used in the calculations. This result is not surprising since we do not estimate the surface tension of a bulk but rather of small clusters.

Figure 5

Surface tension of MC results compared to the surface tension of MD results at different low temperatures.

Figure 6 shows the density profile as a function of temperature. The liquid density has been estimated using the MD simulations which was used to estimate the vapor density from the scaled free energy.

Figure 6

The density profile as a function of temperature.

SPC/HW-heavy water model

Figure 7 shows the surface tension as a function of temperature using molecular dynamics. Since the surface tension vanishes at the critical temperature, we fitted the surface tension data at high temperatures to a line. Using this method, the intercept of the line with the temperature axis is just the critical temperature. The estimated critical temperature using the MD simulations is 671.966K compared to the experimental value of 643.847K.

Figure 7

The surface tension results of the SPC/HW-water model using the Molecular Dynamics (MD), versus temperatures T.

Table 3 shows the values of the free-energy differences of the SPC/HW-heavy-water model at three different temperatures as the number of molecules inside the cluster varies from two molecules up to 75 molecules.

Table 3

The free-energy differences of the SPC/HW-heavy-water model at 260K, 280K, 300K.

N	−δF
	T = 260K	T = 280K	T = 300K
2	6.81	6.01	5.21
3	11.31	9.81	8.71
4	15.31	12.71	11.51
5	11.01	10.41	9.61
6	11.51	10.01	9.21
7	12.31	10.91	9.81
8	12.11	10.81	10.11
10	13.01	11.51	10.31
12	13.71	11.41	10.71
14	13.81	11.81	10.81
18	14.21	12.21	11.11
25	14.71	12.61	11.21
50	15.01	13.11	11.71
75	15.31	13	11.81

Table 4 shows the intercept, the slope, and the excess surface entropy/k per molecule at different critical temperatures of the SPC/HW-heavy-water model. Our results of calculating the excess surface entropy/k per molecule is close to the experimental value of 1.5 [41].

Table 4

The values of the slope and the intercept of the scaled free energy of SPC/HW-water model to different estimated critical temperatures.

Method name	T_C	Slope	Intercept	Ω
TC:MD	671.97	-6.07216	11.07387	1.88343345
TC:MC	775	-4.73935	8.735688	1.47002986
TC: experiment	643. 89	-6.57645	11.94701	2.03985086

Figures 8(a) and (b) show the free-energy differences as a function of the number of molecules raised to the power minus one-third at three different temperatures before and after scaling. We noticed that this model's results scale nicely to the real critical temperature of heavy water of 643.89K [17]. One can calculate the surface tension and the vapor density of heavy-water at any low temperature from this figure.

Figure 8

(a) The free-energy differences versus the number of molecules at 260K, 280K, and 300K. (b) The scaled free-energy differences to .

It is worth mentioning that the scaled free-energy difference of heavy water is much better than the scaled free-energy difference of water. Most probably, this is due to the natural structure of each model. The OPC-water model is a four-site model, while the SPC/HW-heavy-water model is a three-site model. This means that the number of Monte-Carlo steps needed to simulate the water model should be greater than the number taken in this work.

Figure 9 shows the density profile as a function of temperature. The liquid density has been estimated using the MD simulations which was used to estimate the vapor density from the scaled free-energy difference.

Figure 9

The density profile as a function of temperature.

Conclusion

We calculated the free-energy differences of the OPC-water model and the SPC/HW-heavy-water model in this work. The free-energy differences are calculated between ensembles A and B. Ensemble A consists of ()-molecules plus a free probe, while ensemble B consists of n-molecules of the same type. The free-energy differences are calculated at three different temperatures, then scaled to . The scaled free-energy differences show a straight line when plotted versus the number of molecules raised to the power minus one-third. We extracted some important thermodynamic properties from the straight line such as the surface tension, the excess surface entropy/k per molecule, and the vapor density at low temperatures.

The scaled free-energy difference of heavy water shows better scaling than that of the water model. Since the water model is a four-site model, this suggests that the water model simulations should be rerun with higher Monte-Carlo steps.

Regarding the surface tension results, our simulations using molecular dynamics show that the OPC-water model is more accurate around the room temperature, while the SPC/HW-heavy-water model is more accurate around the critical temperature.

Declarations

Author contribution statement

Khetam Khasawinah; Zain Alzoubi & Abdalla Obeidat: Conceived and designed the experiments; Performed the experiments; Analyzed and interpreted the data; Contributed reagents, materials, analysis tools or data; Wrote the paper.

Funding statement

Abdalla Obeidat was supported by () [674/2020].

Data availability statement

Data will be made available on request.

Declaration of interest's statement

The authors declare no conflict of interest.

Additional information

No additional information is available for this paper.