The validity of the potential model in predicting the structural, dynamical, thermodynamic properties of the unary and binary mixture of water-alcohol: Methanol-water case Open Access

The study of alcohol-water mixtures has received high theoretical and experimental attention in the physical and chemical sciences over the past several years.^1–26 It is also commonplace in numerous engineering, industrial, medical and biological applications.^27,28 Many computational investigations have been undertaken to examine the dynamic and structural properties of pure methanol^23,29–33 and its aqueous solutions, and those properties seem to depend strongly on the molecular potential model used in the simulation.

The purpose of this work is to judge the best potential models for determining the general properties of the water-methanol mixtures. To validate the potential model, many thermodynamic, structural, mixing and transport properties should be reproduced quantitatively or qualitatively when compared to experimental published data. Each potential model was constructed and parameterized for different purposes.

Complementary experimental and theoretical studies on liquid water and its mixture with alcohol have been recently undertaken. Among them, Soetens and Bopp checked the ability of two flexible models (PHH (methanol) and BJH (water)) in producing mixing properties of water-methanol mixture which are already satisfactory for the pure liquid methanol and water.¹⁸ Many potential models have been developed to reproduce specific properties of methanol. van Leeuwen and Smit showed qualitative differences between five models of methanol (i.e. J1, J2, H1, H2 and van Leeuwen and Smit.³³ L Bianchi et al studied the structural properties of methanol, using six-site APR6 potential model and compared it with the extended three site H1 model and with the corresponding experimental values.²⁹ Edgar et al applied molecular dynamic simulations to study the thermodynamic, dynamical and dielectric properties of methanol-water mixture, using the OPLS-AA potential model combined with two different water models (SPCE and TIP4P).⁶ OPLS was developed by William L. Jorgensen et al for liquid simulations. In the OPLS-AA model, each atom has its own parameters including all the hydrogen atoms in the methyl group.³⁴ In our (forthcoming) paper we compared the performance of nine models of methanol (J1, J2, H1, H2, B3, OPLS-AA, TraPPE-UA, van Leeuwen and GROMOS) in reproducing the different properties of pure methanol such as surface tension, diffusion, density, and the radial distribution functions. TraPPE-UA treats the hydrogen atoms in the methyl group as a dead load.³⁵ All the mentioned properties were carried out at 11 temperature points between 200 and 300K with step of 10K. We found that both the TraPPE-UA and OPLS–AA potential models show good performance in predicting the mentioned properties. Consequently, our interest in this work is to compare the performance of these two models in predicting the above listed properties in the binary mixtures. All simulations were performed with fixed system size of 600 molecules in the unary system and a total of 1,000 molecules in the binary system. In the current work, our first goal is to check the performance of the TraPPE-UA and OPLS-AA force fields through testing two thermodynamic properties, i.e., density and surface tension, in the unary and in the binary systems. The dynamical property is also calculated to check the accuracy and performance of the two latter force fields in the same systems, using either the Einstein equation (EE) or the Green-Kubo (GK) method.^36,37 For the pure systems the diffusion coefficient was computed using the Einstein equation (EE). Our results are slightly different from those of Refs. 6 and 10. The conflicting and high fluctuating results for the diffusion coefficient obtained from Einstein equation (EE) in the binary system led us to use the Green-Kubo (GK) approach. In addition to calculating the above properties, we examined the validity of TraPPE-UA and OPLS-AA models by calculating the structural properties, specifically the radial distribution function of methanol and its aqueous solutions as a function of temperature for the unary system and as a function of methanol mole fraction in the binary system.

The paper is organized as follows: in section II the simulation methodology is introduced and details on the intermolecular potential models are given. In section III and IV we present and discuss the simulation results for the thermodynamic properties and dynamical properties of unary and binary systems, respectively. The study is completed by calculating the structural properties of the studied systems. The last section is devoted to conclusions and directions to future work.

All simulations were performed using GROMACS package.³⁸ The molecules are confined first inside a box of about 3x3x10 nm³ (see Figure 1) then the box is inserted in the middle of a bigger box of 3x3x30 nm³ to assure vacuum in both sides for the binary mixture with 1000 molecules. For the unary systems, the size of the slap was 3x3x3 nm³ while the whole simulation box is of 3x3x12 nm³ with a total number of molecules of 600. The density map of the binary water-methanol however was confined to a sphere of radius of 2nm. Our systems were equilibrated under pressure for 1ns using Berendsen algorithm³⁹ at 1 bar, followed by 4ns equilibration under a constant volume using Nose-Hoover thermostat.⁴⁰ The data were collected after 6 ns. The appropriate cutoff radius of interaction was taken from our previous work on methanol,⁴¹ which is 1.3 nm for OPLS-AA and 1.4 nm for TraPPE-UA for the Leonard Jones potential and for the particle mesh Ewald (PME) long range interaction. Periodic boundary conditions have been applied in all directions.

While analysing the data to compute the diffusion coefficient of binary systems, we found that the EE method did not converge to a straight line even when we increased the time of equilibration and time of collecting data to 30 ns. Hence, we were obliged to use the Green-Kubo method to calculate the diffusion coefficient, which is valid for relatively short time ∼100ps, but we had to save the data at every step of 2 fs, which requires a very large hard disk space. Periodic boundary conditions were applied in all directions.

In this work, we chose the best two potential models of methanol based on our unpublished work; the Transferable Potentials for Phase Equilibria unite-atom (TraPPE-UA) and the Optimized Potentials for Liquid Simulations all-atoms (OPLSS-AA) force fields. Regarding the water potential models, we chose the TIP4P-water and the SPC/E-water models due to their validity of reproducing most of the real water properties. In TraPPE-UA model, the methanol molecules are represented by three interacting sites located on the methyl group (CH₃) and on the hydrogen and oxygen of the hydroxyl group (OH). In OPLS-AA model, each atom in methanol is represented by one interacting site. In both OPLS-AA and TraPPE-UA force fields methanol is considered to be rigid. The total pair intermolecular potential between atoms i and j is the sum of Lennard-Jones (LJ) and Coulomb potentials

Where ϵ_ij and σ_ij 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: σ_ij = (σ_i + σ_j)/2 and $ϵij=ϵiϵj$⁠.

We have computed a number of thermodynamic, transport and structural properties of liquid methanol and its binary mixture with water for two different potential models of methanol and of water.

The most important quantity when one advising a new model for soft-matter is the vapor-liquid phase diagram. Usually most published model data start from room temperature up to an estimated critical temperature. At low temperatures, the only rigorous method of estimating the vapor density may be through scaled model.⁴² We followed the usual method of estimating liquid density by inserting all the molecules in a slab within a larger empty box and fitting the profile density to a hyperbolic tangent function after assigning zero to the vapor density

Where ρ_l, ρ_v are the bulk densities of the liquid and vapour respectively, z₀ is the position of the Gibbs dividing surface and d is the width of the interface. The density of pure methanol for both the TraPPE-UA and OPLS-AA is shown in Figure 2. The results indicate that the TraPPE-UA potential model shows a better agreement with the experimental values⁴³ compared to the OPLS-AA. This result is not surprising since TraPPE-UA was modeled originally to give the best vapor-liquid phase diagram.

The very good performance of the TraPPE-UA potential model in predicting the density of pure system of methanol led us to examine its accuracy and performance in its aqueous solution. Experimental results⁴⁴ and MD simulations of the total and partial densities from methanol-water binary mixture are displayed in Figure 3 including the densities of the pure systems which are shown at the end points. It is well-known that the density of methanol is lower than that of water; consequently, the density of the mixtures must decrease as the mole fraction of methanol increases regardless of the model. Both the TraPPE-UA and OPLS-AA potential models of methanol with the TIP4P-water give the trend in total density similar to the experimental value but the TraPPE-UA potential model shows a better agreement with the experimental results as shown in Figure 4. The total density of the mixture was also compared with the relation of mixing rules, i.e., ρ(x) = (1 − x)ρ_lw+x ρ_lM where ρ_lw, ρ_lM are the densities of pure liquid water and pure liquid methanol respectively, and x is the methanol mole fraction. From the results, one can deduce that it is not necessary to estimate the total density as a function of mole fraction, but rather to estimate them from the mixing rules and the calculated end points. One more point should be addressed here about the little discrepancy in calculating the total density, this discrepancy might due to the non-additive feature of the intermolecular potential, i.e., the Lorentz rule of mixing might be no longer satisfied,^45–48 and this non-additivity might affect the mixing properties of the mixtures.⁴⁹ 

At mole fraction of 0.3 for methanol, we were unable to estimate the liquid density of both water and methanol due to their peculiar profile shapes (the total density was easy to estimate, the density of water was estimated by taking the average of the density points, while the density of methanol was evaluated by subtracting the water density from the total density), refer to the right panel of figure 5. From Figure 5 we can predict that the structure at x = 0.1 is a core-shell (CS) (Figure 5i) which is a well-known structure in soft matter characterized by one phase surrounding the other phase. This is clear since methanol profile has peaks (like horns) near the borders of the slap indicating higher density. Increasing the mole fraction of 50% and above (Figure 5k,l), the profiles are smooth (no horns), indicating that methanol is totally homogeneous inside water and vice versa, and therefore called well-mixed (WM) structure. Near x = 0.3 (Figure 5j), the density distribution is not homogeneous having two horns with different values. To eliminate ambiguity about this structure, we performed another experiment by placing the molecules inside a sphere of radius of 4nm with simulation time of about 10 ns; then the density map was plotted to figure out the structure of the mixture as shown in Figure 5a-h. It is clear from the figure that the structure of methanol at x < 0.3 is core-shell and for x > 0.5 the structure is well-mixed, but at x = 0.3 specifically, the structure is neither. It can be explained that both solutions start to segregate; the density on the right side in Figure 5b is clearly higher than that of the left side. The new structure is phase-separated (PS) which looks like Russian-doll (RD). It is worth mentioning that the number density contour maps and the density profiles of methanol around the symmetry axis (z = 0) for the above results were computed at different mole fractions x_MEOH=10, 30, 50 and 90% of methanol at T=300K.

From Figure 5(a-d) we found that the density profiles at all mole fractions of methanol seems to be approximately symmetric about the midpoint L_z/2. This gives an indication that the mixture (i.e. inhomogeneous system) is properly equilibrated. We also observed that the preferred structure at low mole fractions of methanol is the core-shell; this structure is converted to a well-mixed structure by increasing the mole fraction of methanol. At x_MEOH=0.3, at the first sight the distribution is somewhat similar to the phase separated structure (PS). This led us to examine whether this structure is regular or not at different temperatures. In the current work, we chose different temperatures 210, 240, 270 and 300 and checked the distribution of methanol into water throughout density contour map plot at x_MEOH=0.3, 0.5, 0.7, and 0.9. We found that the PS structure appears again at T=270K for x_MEOH=0.7. This confirms that this structure does not regularly happen at all temperatures and mole fractions of methanol. Figure 6 illustrates this unusual PS structure at the latter concentrations of methanol at T=270K. In our forthcoming papers in binary mixtures of ethanol and propanol with water using the same force fields studied in this paper, we will show the effect of carbene (CH2) on the structure of mixtures as a function of mole fraction and temperature. The Russian-Doll structure was clearly demonstrated in recent work of Hrahsheh et al⁵⁰ on binary pentanol-water mixture. More details about the distribution of different type of molecules in the mixture are discussed in the supplementary material.

Surface tension is the thermodynamic property, which plays an essential role in nucleation and very important in industry. Due to high fluctuations in pressure, we ought to equilibrate our systems for a long time, sometimes for 10 ns, and then collect the data for another 10 ns. The surface tension has been calculated using the following relation:

Where P_αα is the αα component of the pressure tensor, L_z is the box length on z direction, ϵ and σ are the Lennard-Jones parameters, and r_c is the cutoff radius.

The effect of temperature on the surface tension was examined for the unary system of methanol using TraPPE-UA and OPLS-AA potential models. It is clear from Figure 7 that the surface tension of methanol decreases with temperature in both models, which agrees with the results of other authors.^41,42 From Figure 7, we notice that the TraPPE-UA model shows a better agreement with the experiment⁵¹ at relatively high temperatures; whereas OPLS-AA model shows a better agreement with the experimental results at very low temperatures.

Surface tension was also calculated at different mole fractions of methanol from x_MEOH=0.02 to 0.1 with steps of 0.02 followed by steps of 0.1 up to x = 1. The comparison between TraPPE-UA and OPLS-AA force fields in reproducing the surface tension in the binary mixture is shown in Figure 8. From the figure, it is evident that both the TraPPE-UA and OPLS-AA potential models fit the experimental data⁵² very well when they are mixed with TIP4P-water, but the TraPPE-UA model shows a better overall agreement especially at a very low concentration of methanol. It is worth mentioning that the fast drop in the surface tension when increasing the methanol composition from 0 to 20% is mostly due to the methyl group (CH3). We propose an explanation to this behavior with the aid of the density contour maps of methanol as follows: the structure of the drop is a core-shell, where most of the molecules on the outer-shell are mainly methanol making methanol a dominating factor in surface tension. Moreover, it is shown in the supplementary material that the methyl groups are on the surface while the hydroxyl groups point inward the surface, which indicates again that methanol is controlling surface tension. Even though by increasing the methanol concentration the structure changes to a well-mixed, but the outer most shell of the sphere mainly consists of methanol which explains why the surface tension of the mixture changes slowly to its value of methanol, please refer to the supplementary material.

It is well-known that the diffusion coefficient is a sensitive quantity to the size of the box and time of the simulated system. In our systems, the diffusion coefficient has been strongly affected by the period of simulation time. In principle, the Einstein relation and the Green-Kubo methods should produce the same results. As we will see later, even though the Einstein relation was perfect in reproducing the diffusion coefficient in the unary systems, it fails badly to estimate the coefficient for the binary systems even with simulation time up to 30 ns. This will be explained below.

The self-diffusion coefficients were calculated from the time dependence of the mean square displacement (MSD) using Einstein equation (EE)

Where D is the diffusion coefficient, r(t) and r(0) are the positions of the molecule at time t and 0, respectively.

The Einstein equation (EE) has been used to compute the diffusion coefficient in both the unary and binary mixture systems of methanol with water. For pure methanol the self-diffusion coefficient was calculated in the temperature range from 200 to 300K. All mean square displacement versus time plots were linear; this means that our results for diffusion coefficient converge and the time of simulation is sufficient. The MDS for the two potential models were plotted versus $|<r(t)2−r(0)2>|$ which yields a straight line with the slope equals to 6D. Figure 9 shows the diffusion coefficient as a function of temperature compared to experimental data.⁵³ It is clear from the figure that the TraPPE-UA and OPLS-AA potential models follow the experimental trends. Both models predict higher values than the experimental results, but the TraPPE-UA model shows a better agreement over a wide range of temperatures up to 280K. The above conclusion is confirmed when estimating the activation energy of diffusion through Arrhenius plot. Figure 10 shows Arrhenius plot of diffusion by assuming that the diffusion is proportional to $∼e−EkBT$ where E is the activation energy of diffusion. It is clear that the TraPPE-UA model not only predicts the experimental data, but also gives the closest activation energy of diffusion.

Regarding the diffusion coefficients of species in binary mixtures, the Einstein equation fails to reproduce the experimental data. Figure 11 shows the variation of MSD as a function of time up to 6ns. From the figure, we see that the MSD is linear, but in three different regions, the slope changes from relatively high value to a lower value then back again to higher value. The data should be collected from the last region since EE is valid for the long time. For this reason, we increased the time to 20ns, but the linearity of MSD as a function of time did not converge to a one linear function. After collecting the data, the diffusion is plotted as a function of mole fraction of methanol (Figure 12), and the results were completely random.

Instead of repeating all our MD simulations over a long time (maybe > 50 ns) to get a sensible value of diffusion coefficients, we used Green-Kubo (GK) method which requires a short simulation time. Even though the simulation time required to calculate the diffusion coefficients using GK is short, which is an advantage of applying GK, we ought to save every frame at 2fs up to 100 ps which requires a very large disk space.

The Green-Kubo formula (GK) is the method of choice for calculating the diffusion coefficient in short time simulations, which is given by:

Where v_k(t) is the center of mass velocity vector of molecules at some time t, v_k(0) is the same but at time t = 0, ⟨…..⟩ denotes the ensemble average, and N_i is the total number of molecules of type i. To calculate diffusion coefficient we applied the criterion of Zlenko.⁵⁴ Figure 13 shows the self-diffusion coefficient of methanol and water for all different models combinations. It is evident that the MD simulated diffusion coefficient produced from the OPLS-AA/SPCE-water adequately fits the experimental data over the entire range of composition. This indicates that combining methanol with SPCE-water is more valid in assessing methanol diffusion coefficients.

We now proceed to investigate how the molecular structural features of the pure methanol changes as the temperature increases for both models throughout the position of the first peak, the coordination number, the structure function and the distinct radial distribution function. All the radial distribution functions RDFs (i.e. g_oo(r), g_oH(r), g_CO(r), g_CC(r), g_HC(r), g_HH(r)) were calculated for the OPLS-AA and TraPPE-UA potential models. Figures 14 and 15 show the radial distribution function of O-O for TraPPE-UA and OPLS-AA respectively. From the figures we notice that the position of the first maximum does not change with temperature and it is about 2.7A° as the experiment proposes⁵⁵ and the height of the first maximum peak decreases as the temperature increases. This behaviour is expected since by decreasing temperature, the molecules will come closer to each other, and the network becomes more structured, and as a consequence the coordination number will increase as well. The first peak of g_oo(r) shows a reduction in of magnitude as the temperature increases as a result of breaking the hydrogen bonds in the liquid methanol for both potential models. The coordination number is of about 2 at 300K and it increases slowly as the temperature decreases. But generally speaking the value is slightly lower than 2, which suggests that the network of methanol even at low temperature is a chain network. The distinct radial distribution function will be explained later when we discuss the mixture of methanol-water, but most of the partial distribution function (PDF) contributions come from g_OO and g_CC as will be explained in the following paragraphs. For the other radial distribution functions, please refer to the supplementary material of this manuscript. We noticed from the figures that the TraPPE-UA potential is more structured, and this is not unexpected since the force field contains three atoms, while the OPLS-AA has 6 atoms.

Figure 16 shows the evolution of the first peak with the temperature increase. The position of the first peak in the OPLS_AA and TraPPE-UA models is nearly around 2.7A°. From the figure we notice that the height of the first peak changes linearly with temperature, but the variation becomes high for the TraPPE-UA force fields, starting from about a value of 5.3 at T = 200K to a value of 4 at T = 300, while the value of the peak in OPLS-AA changes from 3 to 2.5 as the temperature changes from 200 to 300K.

In the light of the above results, we can say that the g_oo(r) is more structured in the TraPPE-UA models as the temperature decreases (i.e. the interaction between methanol molecules is stronger in this potential model).

For the binary mixtures, it is worth showing the RDF of O-O of all force fields used in this work at 300K to understand what is going on when increasing the concentration of methanol. Figure 17 shows such behavior at x = 0 (pure SPCE and TIP4P-water) and at x = 1 (pure methanol of TraPPE-UA and OPLS-AA). The corresponding coordination numbers will be shown later, but it is about 4 for water which is tetrahedral as suggested in the literature.

The corresponding distinct RDF is shown in Figure 18; on the other hand, the distinct RDF is calculated as follows:

Where ρ_m is the total number density of all atoms and S_d(k) is the intermolecular contribution as a function of wavenumber k which is given as:

Where f_α is the X-ray form factor approximated to be equal to the total number of electrons in each group, which is a valid approximation for low values of k, $Cα=Nα∑iNi$ is the concentration of atoms of type α. and $ŝαβ$ is defined as:

We notice that both potentials reproduce the experimental curve for water models. For methanol, both models reproduce the experimental data but with advantage for the OPLS-AA regarding the first minima, while the TraPPE-UA works better for predicting first maxima. We notice from the figure that the first peak occurs about 2.6A°, which means that the O-O PDF is responsible for the appearance of the first peak, while the C-C RDF is responsible for the second and third peaks; see the positions of maximum points of the g_CC in the supplementary material. Even though all models reproduce the experimental data quite well, the small deviation from experimental data might be due to approximating the form factor to be the same number as that of electrons in the atom. Although such approximation is valid only for small values of k, and in calculating the distinct RDF, the integration is taken over the whole k-space. The effect of this approximation will be shown later when we discuss in detail the total structure factor.

Figure 19 shows the oxygen-oxygen PDF of binary mixtures at a methanol mole fraction of 0.6. It is clear that when the water oxygen is involved, the results do not depend on the model chosen for methanol since the RDF of pure water is almost the same for both SPCE and TIP4P-water. For O-O PDF of methanol, we notice that the first maximum peak depends greatly on the model chosen, and the first maximum peak is higher for TraPPE-UA than for the OPLS-AA as for pure systems, but the height decreases sharply by increasing water concentration.

Figure 20 shows the variation of oxygen-oxygen RDF functions as the mole fraction of methanol increases. The figure shows only the OPLS-AA methanol mixed with TIP4P-water (left panel (a-d)), and mixed with SPCE-water (right panel (f-h)). For TraPPE-UA, please refer to the supplementary material. These figures are important to calculate the coordination number which might be thought as the number of hydrogen bonds which will be discussed in detail later. As an example regarding the coordination number, if the methanol behaves as solute the O_M-O_W at x = 0.2 indicates that the hydroxyl oxygen has 2.3 water molecules within a 3.34A° radius, and the carbon is surrounded by 12 water molecule within a 5.2A° as can be figured out from the C_M-O_W. On the other hand if the methanol behaves as a solvent (x = 0.8 as an example), the O_M-O_W indicates that the hydroxyl oxygen has 0.65 water molecule within a 3.42A° radius, and the carbon (from C_M-O_W) is surrounded by 1.6 water molecules within radius of 4.9A°. It is also clear that by increasing the methanol concentration, the first maximum grows substantially while the second maximum remains almost unchanged. One might notice that by adding more methanol, the first minimum becomes lower, which is an indication that the first shell is separated from the second shell.

Figure 21 shows how the maximum value of the first peak of oxygen-oxygen RDF changes with mole fraction of methanol. It can be seen that the relation between the values of the first maxima of the O-O RDF is linearly proportional to the increase of methanol concentration up to x = 0.7 regardless of the type of oxygen. This does not necessarily mean that the coordination number should increase as we will see later.

Figure 22 shows the corresponding distinct RDF of the above figures but with TraPPE-UA methanol compared to experimental data at the end points (x = 0 or 1, see Figure 18). The figure shows how the structure of water (x = 0) changes slowly by increasing the concentration of methanol till it reaches the pure structure of methanol (x = 1). By increasing x, the second and third peaks start to appear, and this is a clear indication that the second and third peaks come from carbon.

The total structure factor S(K) obtained from the MD simulation for the mixture of methanol-water at T = 300K is shown in Figure 23 and compared with the structure factors obtained from X-ray and neutron diffraction.⁷ The structural change is obvious in the mixture as the mole fraction of methanol increases. From the figure we notice that the S(K) for pure water is almost the same as the experimental data beyond the third maximum, and the positions of the first three maxima are predicted well by the models. As the methanol concentration increases, the structure function coincides better regarding the first maxima, but the second maxima deviate to a lower value till x = 0.61. We also see that the phase of the third maxima starts to deviate outwards. This means that the position of the C-C RDF also deviates to the right compared with the experimental data. The most awkward result is at x = 0.8. We notice that all mixed models almost coincide with each other, but the height of the first maxima is higher than the experimental data, and the second maximum is lower. At x = 1 for pure methanol, the results are almost perfect when compared to experimental data, except for the phase shift of the first minima. The total scattering structure function of experimental pure water shows a double-peak first maximum, which is something SPCE and TIP4P-water could not reproduce in contrast to results of Galicia-Andres’ et al.⁷ Upon adding methanol, the double-peak first maximum occurred at x = 0.2, and by increasing the concentration of methanol, the first maxima increase slowly and the second maxima decrease till it becomes a shoulder. Same behavior also occurred regarding the fourth maxima. For the pure methanol, the TraPPE-UA shows a better agreement over the OPLS-AA regarding the appearance of the shoulders, but both models show a shift in the first minima inward. The changes in the total scattering structure function is due to the change of the network of molecules from tetrahedral for pure water to short chain-like for methanol passing through cross-linked network at x∼0.7 at which the simulated results show the biggest discrepancy with experimental data. The main reason for not reproducing the results of Galicia-Andres’ et al is due to the X-ray form factor. In our calculations we approximated the form factor to be equal to the number of electrons in each group, which is valid, as explained earlier, for low values of wave number (k); while the form factor in Galicia-Andres’ et al work was calculated using a sum of five Gaussians with coefficients taken from Ref. 56 as:

The results of Galicia-Andres’ et al were reproduced when the form factor is substituted as a Gaussian function of wave number. In this scheme, see the supplementary material for more details, we notice that the appearance of the doubly first peak in pure water coincide with the experimental data and TIP4P-water shows a better agreement over the SPCE-water. Increasing the methanol mole fraction, TraPPE-UA is better regarding the first maxima while OPLS-AA is better regarding the second maxima when compared to experimental data up to x < 0.5. For x > 0.5, OPLS-AA is better in reproducing the experimental data regarding the first two maxima but TraPPE-UA works better for predicting the second minima. For pure methanol, both models work very well compared to experimental data but with advantageous to TraPPE-UA.

The complexity in the structure of methanol-water mixture is due to the hydrophilic head in methanol structure, this makes it capable of forming hydrogen bonds with itself and with other polar molecules like water. The number of hydrogen bonds can be approximated through the coordination number, which is defined as the area under the RDF curve up to the first minimum of the g_OO, as:

Where ρ is the number density and the integration starts from r = 0 and ends at first minimum of the site RDF g_αβ(⁠$r′$⁠). The studied coordination numbers are (O_W-O_W), (O_M-O_M), (O_W-O_M) and (O_M-O_W), where the O_W and O_M correspond to the water oxygen and methanol oxygen respectively.

Figure 24 shows the number of hydrogen bonds between the O_W-O_W and O_M-O_M as a function of mole fraction of methanol. Figure 24(a) shows the effect of water type using OPLS-AA for methanol, while Figure 24(b) shows the effect of methanol type when mixed with SPCE-water. From the figures we noticed that the type of water or methanol does not change the hydrogen bond, and that the variation of hydrogen bonds with increasing the mole fraction of methanol is linear. One might deduce that the network structure is tetrahedral when x = 0 (pure water) and the network structure starts to change slowly by increasing x till we get a chain network when x = 1 (pure methanol). Even though the total network structure depends highly on the number of hydrogen bonds, one might deduce from the figures by ignoring the mixing hydrogen bonds that the network structure at x∼0.7 is cross-linked (n_HB = 3), where each type of molecules share roughly the same contribution of 1.5 coordinated molecules.

Figure 25 shows all mixing hydrogen bonds as a function of methanol mole fraction. We notice that both coordination numbers calculated from O_W-O_W and O_M-O_W RDF’s decrease as x increases. While for O_W-O_M and O_M-O_M, the hydrogen bonds increase as x increases. Since the variation is linear, one might deduce that the relation between the number of hydrogen bonds of O_W-O_W and O_M-O_M follows the mixing rules as $nHBtS=xnM+1−xnW$⁠, where $nHBtS$ is the total number of hydrogen bonds of the same type (MM or WW), n_M is the number of methanol hydrogen bonds and n_W is the number of water hydrogen bonds, see Figure 24(b). The same is true for the total mixing hydrogen bonds $nHBtM$ which can be written in terms of water-methanol hydrogen bonds $nW′$ and methanol-water hydrogen bonds $nM′$ as $nHBtM=xnW′+1−xnM′$⁠. One might notice that for pure water, the molecules form four hydrogen bonds. Adding methanol i.e., starting from left to right, the four coordinated molecules decrease linearly with increasing the concentration of methanol, while the two and three coordinated molecules increase with increasing x_MEOH. If one starts from right to left, i.e., adding water to methanol, number of hydrogen bonds decreases linearly, and at the same time increases in the three coordinated molecules species (sees O_M-O_W).

We end this section by giving a general review about the above results. For the unary systems of methanol, we saw that by decreasing the temperature, the first maximum increases, and the distinct RDF shows that the TraPPE-UA coincides better with the experimental results regarding the first maxima, while OPLS-AA agrees well with the experimental results regarding the first minima. For the binary systems, our two different methods of calculating the structure factor through the form factor shows that the use of Gaussian function for the wave number (k) as a substitution for the form factor gives better agreement with experimental data rather approximating the form factor to be the number of electrons in the molecule. Again, the TraPPE-UA shows better agreement with experimental values when calculating the structure factor at 300K. In calculating the number of hydrogen bonds, both models almost give the same values.

A series of molecular dynamic simulations for pure liquid methanol and its mixture with water have been performed using two different methanol models (TraPPE-UA and OPLS-AA) with two different water (SPCE and TIP4P) potential models. The results are compared with the experimental thermodynamic, dynamical and structural data in order to examine the validity of presented potential models in producing the structural, dynamical, and thermo-dynamical properties in both the unary and binary mixture of methanol-water. Density and surface tension have been computed in order to reveal the performance of the above mentioned models in the pure and mixed methanol with water.

The total density of the mixtures as well as the number of hydrogen bonds of the same type and of different types follow the relation of mixing rules.

Regarding the structure factor of mixtures which depends on calculating the form factor, we showed that better agreement with experimental data can be constructed by calculating the form factor using the sum of five Gaussians rather assuming it to equal number of electrons in each group.

The overall conclusion is that the TraPPE-UA potential model works well in reproducing surface tension, density and diffusion in both the unary and binary mixture of methanol with water. Unfortunately, the results of diffusion based on Einstein equation (EE) randomly fluctuated over an entire range of methanol mole fractions regardless of potential model we used. The results obtained by Green-Kubo (GK) method give the trend in diffusion coefficient similar to the experimental results in all systems (i.e. the unary and mixed system). Furthermore, the OPLS-AA model of methanol combined with SPCE water is more valid in producing methanol diffusion coefficient and structure function (for x > 0.5).

Finally, it is worth mentioning that the validity of potential model in predicting the desired properties in unary system does not mean that the model is valid in predicting the same properties in the binary mixture.

See supplementary material for the structure of methanol inside the mixture, and for the radial distribution functions of carbon with other atoms and the effect of Gaussian form factor to the total structure factor.

The authors would like to thank Professor Akram Rousan for his valuable comments and proof reading the article. This work has been financially supported by Jordan University of Science and Technology under the grant #27/2016.