Molecular dynamics simulations of the properties of water-methanol mixtures. Effects of force fields
M. Cruz Sanchez H. Dominguez O. Pizio

A puzzle on the existence of Lyapunov functions for Limit cycle system





    [
    
===========================================================================




Isothermal-isobaric  molecular dynamics simulations 
are used to examine the microscopic structure and some properties
of water-methanol liquid mixture. 
The TIP4P/2005 and SPC/E water models are combined 
with the united atom TraPPE and the all-atom force field model for methanol. 
Our principal focus is to evaluate the quality of predictions of
different combinations of model force fields concerning the composition dependence of 
basic properties of this system.
Specifically, we explored the composition effects on density, excess
molar volume and excess entropy, as well as on the surface tension and static
dielectric constant. In addition, the
structural properties are described in terms of the coordination numbers 
and the average number of hydrogen bonds between molecules of constituent species.
Finally, the composition dependence of self-diffusion coefficients of the species
is evaluated. All theoretical predictions are tested with respect to experimental data.


ater-methanol  mixtures, 
mixing properties, surface tension, molecular dynamics simulations

61.20.-p, 61.20-Gy, 61.20.Ja  





§ INTRODUCTION


This manuscript is the first part of our two-stage project that involves water-methanol
mixtures. Namely, in this first part, we present a set of results coming from the 
isothermal-isobaric (NPT)
computer simulation concerning composition changes of density, excess mixing volume and
entropy, first coordination numbers of species and average number of hydrogen bonds.
Moreover, we explore the behavior of the surface tension and the dielectric constant  
on composition, as well as the self-diffusion coefficients of species.

The forthcoming, second part of the project, is devoted to the exploration of changes
of all the properties mentioned above, brought by addition of NaCl salt to water-methanol
solvent. Again, the NPT computer simulation technique will be applied. One of
the principal issues we would like to address is to evaluate the validity and quality of 
theoretical results coming from different combinations of force fields with respect 
to experimental data.  

For specific purposes of our project, it is worth mentioning that molecular dynamics computer
simulations have been widely applied to mixtures of water with various solvents.
Most frequently studied seem to be the mixtures of water with alcohols, specifically
the water-methanol mixtures,
see, e.g., <cit.> 
and references therein for a rather comprehensive account of the previously applied modelling. 
Our present report has been inspired by most recent contributions
concerning the system in question <cit.>.
On the other hand, important experimental observations concerning this system and used 
in the present study have been discussed in <cit.>.

The majority of computer simulation studies cited above have been focused on the application
of a single combination of force fields describing each constituent species, water and
methanol with the exception of <cit.>. This latter work involved two methanol
force fields, namely the TraPPE model <cit.>  and OPLS-all atom model <cit.>.  
As concerns water, the SPC/E and TIP4P models were chosen, see, e.g., <cit.> and <cit.>,
respectively.  Formally, we apply similar strategy. However, in contrast to <cit.>, 
we apply the TIP4P/2005 water model <cit.> rather than the TIP4P version <cit.>.
It is known that the the TIP4P/2005 version provides a better performance of the
microscopic structure of water <cit.>  and yields even better description of various 
properties compared to the TIP4P, see, e.g., table 2 of reference <cit.>.
Moreover, we are specifically interested in using the TIP4P/2005, because it has been quite recently applied to parametrize the description of properties of NaCl aqueous solution, 
see <cit.>. It opens up the possibility to study NaCl solutions with
water-methanol solvent in subsequent work. On the other hand, the gained
knowledge would permit to explore complex solutes in mixed solvents tuning with 
confidence their solubility on solvent composition. 
Some specific molecules of interest in medicinal chemistry, for example curcumin,  
are marginally soluble in water, but dissolve well in, e.g., alcohols or 
dimethylsulfoxide <cit.>.
Additional comments concerning the
methodological difference of our procedure and technical details in comparison with reference <cit.>
are given in the body of the manuscript below.

To summarize, the principal objective of the present report
is to investigate a set of properties of water-methanol mixtures in the entire  interval 
of composition. All the properties are validated by comparison with experimental results.



§ MODELS AND SIMULATION DETAILS


Preliminarily discussing the simulation methodology, for the sake of convenience for the reader, in table <ref>
we list some of the models used in previous studies of water-methanol mixtures and 
related to the present work. Within the united atom model for
methanol, like OPLS/UA and TraPPE (see table <ref>), the CH_3 group is considered as a single
site. On the other hand, the all atom models explicitly involve all the hydrogens of the methanol molecule.




 §.§ Technical details


In this work, we explore the SPC/E model <cit.> and the 
TIP4P/2005 model <cit.> for water.
For methanol, we used two models, namely
the united atom model-TraPPE <cit.>,
and all atom model <cit.>, denominated as  MET/TraPPE, and MET/AA, respectively.
The Lorentz-Berthelot combination rules are used to determine the cross parameters.


The long-range electrostatic interactions were handled by the
particle mesh Ewald method implemented in the GROMACS software package (fourth
order, Fourier spacing equal to 0.12) with a precision of 10^-5.
The nonbonded interactions were cut off at 1.4 nm.
The van der Waals tail correction terms to the energy and pressure were taken into account.
In order to maintain the geometry of the water and methanol  molecules, the LINCS
algorithm was used.

Our calculations were performed in the isothermal-isobaric (NPT) ensemble
at 1 bar, and at a temperature of 298.15 K. We used the GROMACS
software package <cit.>, version 5.1.2.
Concerning the procedure, a periodic cubic simulation box was set up for each system. 
The GROMACS genbox tool was employed to randomly place all particles in the 
simulation box.
To remove the possible overlaps of particles introduced by the procedure of
preparation
of the initial configuration, each system underwent energy
minimization using the steepest descent algorithm implemented in the GROMACS
package. Minimization was followed by a 50 ps NPT equilibration run at 298.15 
K and 1 bar using a timepstep of 0.25 fs.
We used the Berendsen thermostat and
barostat with τ_T = 1 ps and τ_P = 1 ps during equilibration.
Constant value of 4.5·10^-5 bar^-1 for the compressibility of the mixtures was
employed. In the case of pure methanol solvent, the compressibility
was taken to be 1.2·10^-4 bar^-1.
The V-rescale thermostat and Parrinello-Rahman
barostat with τ_T = 0.5 ps and τ_P = 2.0 ps
and the time step 2 fs were used during production runs.
Statistics for each mole solvent composition and various ions
concentration for any of the properties were collected over
several 10 ns NPT runs, each starting from the last configuration of the
preceding run. The time extension for each series of calculations will be mentioned below
in the appropriate place but not less than 70 ns.

While exploring the composition changes, 
the total number of molecules 
was kept fixed at 3000. 
The composition is described by the mole fraction of methanol molecules X_m, X_m=N_m/(N_m+N_w).



§ RESULTS AND DISCUSSION


In order to make things clear from the very beginning, we study four models of water-methanol
liquid mixtures described by the following force fields: TIP4P/2005-MET/TraPPE,
TIP4P/2005-MET/AA, SPC/E-MET/TraPPE and SPC/E-MET/AA. In all cases, the Lorentz-Berthelot
combination rules are used to obtain cross interaction parameters.
The behavior of mixture density on composition, X_m,  is given in figure <ref>.  
Our NPT computer simulation results are supplemented by the experimental data <cit.>. 
It can be seen that the density is perfectly well described by the TIP4P/2005-MET/TraPPE
model. The SPC/E-MET/TraPPE combination of force fields yields only a bit worse
predictions in the interval of intermediate compositions. Two models, the TIP4P/2005-MET/AA
and SPC/E-MET/AA, are less accurate regarding the density dependence on composition, apparently
because the density of pure liquid methanol is not described accurately within the MET/AA model. 
These trends of behavior have been mentioned recently in <cit.>
exploring TIP4P/Ew water model <cit.> combined with methanol models with a different 
degree of sophistication.
On the other hand, the density dependence on composition from the NPT simulations 
by Soetens and Bopp
also agrees with the experimental data perfectly well <cit.>. It is worth mentioning 
that the authors used well established BJH <cit.> and PHH <cit.> 
flexible models for water and methanol, respectively.

The results presented in figures 2 and 3 of <cit.> are slightly off all the observations mentioned above concerning the density dependence on methanol 
molar fractions. Namely, in figure 2 
of that work, the liquid
density of pure methanol from TraPPE and OPLS-AA models at 300 K looks too low compared to the 
values reported in TraPPE database and by other authors. Moreover, the discrepancy between
the simulation predictions and experimental data given in figure 3 of <cit.> (the y-scale
of this figure seems to be not very appropriate to appreciate the total density) is too
big compared to other works. Apparently, trends of behavior of density from this work <cit.>
are affected by the number of particles, system size, method of calculations using 
a rather small liquid slab surrounded by vacuum and by simulation time.
To summarize the discussion of our figure <ref> and the observations of other authors, it seems 
that the dependence of density on composition is one of the important properties, 
but not too demanding computationally, 
if the system size and time of simulations are properly chosen.

Next, in the spirit of previous works from our laboratory <cit.>  and 
of the development by Soetens and Bopp <cit.>, we turn our attention to the mixing 
properties. The excess molar volume, Δ V_mix, is defined as,

    Δ V_mix=V_mix-(1-X_m)V_w-X_mV_m ,

where V_mix is the volume of the mixture at a certain composition, V_w and V_m refer
to the molar volumes of pure water and pure methanol, respectively. Our simulation results 
for Δ V_mix are shown in figure <ref> (a).
As in the case of density, cf. figure <ref>, here we observe again that the TIP4P/2005-MET/TraPPE model
yields the best agreement with experimental results <cit.>. The entire set of experimental points is
well reproduced, slight overestimation of the magnitude of Δ V_mix is observed in the
interval of intermediate compositions and when methanol content is higher than of water.
Nevertheless, the experimentally observed  minimum at X_m=0.5 
is reproduced by simulations. If methanol is described
by all-atom model (MET/AA), the agreement between simulation results and experimental points
deteriorates, the minimum shifts to X_m≈ 0.6. Quite similar dependence 
of Δ V_mix on composition (to the TIP4P/2005-MET/AA force field) has been obtained 
from simulations of flexible BJH-PHH model, cf. figure 2 of reference <cit.>. If one combines the SPC/E
water with each of two methanol models in question, the values for Δ V_mix in almost entire 
range of composition are underestimated, figure <ref> (a). This kind of inaccuracy has been observed 
previously and was discussed in detail in references <cit.>.

Energetic aspects of mixing are given by the excess enthalpy of mixing, Δ H_mix. It is
defined similarly to equation <ref>,

    Δ H_mix=H_mix-(1-X_m)H_w-X_mX_m ,

where H_mix refers to the mixture entalphy whereas H_w and H_m describe the entalphy of each
species in pure state at the same temperature and pressure. From the comparison of the
simulation results in figure <ref> (b) with experimental data <cit.>, we conclude that none of the employed
force field combinations reproduce the experimental trends perfectly well. Concerning the magnitude of
Δ H_mix, the TIP4P/2005-MET/TraPPE and SPC/E-MET/AA models provide a better
description than the TIP4P/2005-MET/AA and SPC/E-MET/TraPPE. However, all the models
of this study predict the minimum of Δ H_mix at X_m in the interval between 0.5 and 0.6
rather than the experiment that shows the corresponding minimum at X_m=0.3. Similar kind
of inaccuracy has been observed and discussed for a set of previously studied 
models <cit.>. Unfortunately, the excess entalphy of mixing from NPT simulations
has not been presented in reference <cit.> for the flexible BJH-PHH model to provide a critical
evaluation of the present modelling. Three data points previously reported for the
excess potential energy from the NVE simulations of the combination of flexible models, 
see table 1 of reference <cit.>, do not provide a definite answer in this respect, unfortunately.
 
The microscopic structure of water-methanol mixtures has been discussed in terms of
evolution of the pair distribution functions on composition in many occasions. However,
quantitative insights about the changes of structure on X_m are usually described in terms
of the first coordination numbers of the species and of the number of hydrogen bonds. Here,
in order to avoid unnecessary repetition, we adopt a similar point of view.
The first coordination numbers of the species are defined from the running coordination numbers
as common,

    n_ij(r)=4ρ_j ∫ _0 ^r  g_ij(R)R^2  R ,

where the upper limit of integration in this equation is taken to correspond to the first mi­ni­mum 
r_min of the pair distribution function g_ij(r) of the species i and j.
A detailed description of the behavior of this property, 
for a set of models related to the present study, was given in <cit.>.
Very recent reports also concern the evolution of the first coordination numbers, 
see, e.g., figure 7 of reference <cit.> and figures 24 and 25 of reference <cit.>.

In order to evaluate the average number of H-bonds between molecules, the corresponding utility
of GROMACS software is applied with default distance — angle criterion. However, for distance
cutoff we used the first minimum of the corresponding pair distribution function,
see also discussion of this issue in references <cit.>. 
Averaging is performed over a piece of or the entire simulated trajectory.



Changes of n_ij=n_ij(r_min) on X_m for oxygens that belong to water or methanol,
OW-OW and OW-OM, and resulting from our simulations, are shown in figure <ref>(a). 
The coordination number n_OW-OW monotonously,
almost linearly,
decreases with an increasing methanol mole fraction starting from the value ≈ 4.5 for pure water.
On the other hand, n_OW-OM increases up to ≈ 3  while approaching X_m=1.
Two combinations of force fields, TIP4P-2005-MET/TraPPE and SPC/E-MET/TraPPE, yield
almost the same results. The crossing point between n_OW-OW(X_m) and n_OW-OM(X_m)
describing “inversion” in the type of predominant neighbors around water molecule 
occurs at X_m≈ 0.65. The crossing point at the same composition characterizes 
the BJH-PHH combination of flexible models for two species, figure 7 of <cit.>. 
On the other hand, evolution of the average number of hydrogen bonds ⟨ n_HB⟩ per water molecule
on the methanol mole fraction from our simulations is given in figure <ref> (a) as well.
Actually the dependences of ⟨ n_HB⟩ on X_m follow the trends of behavior of the first 
coordination numbers. The crossing point describing the predominant number of H-bonds 
between two kinds of molecular species is observed at X_m≈ 0.65. Again, 
two water models (TIP4P/2005 and SPC/E), if combined with MET/TraPPE, yield very similar
predictions concerning the trends of behavior of ⟨ n_HB⟩ on X_m.   

The “inversion” or crossover point in n_OW-OW(X_m) and n_OW-OM(X_m) or ⟨ n_HB⟩
is located close to the minimum of the excess mixing volume, cf. figure 2 (a), as it was
discussed in <cit.>. Finally, we would like to mention that the approximation of
average hydrogen bond numbers by the corresponding coordination numbers is not very
appropriate in the case OW-OW, see figure <ref> (a), in contrast to what was claimed in <cit.>.

Concerning the behavior of methanol coordination number and the average number of hydrogen bonds
in a mixture with varying composition, figure <ref> (b), we would like to mention the following trends.
We explored two methanol models, the MET/TraPPE and MET/AA, combined with TIP4P-2005 water
model. The models yield a similar behavior. In pure methanol, the coordination number is
around three and the average number of hydrogen bonds between methanol molecules is only
slightly less. The methanol coordination smoothly decreases with a decreasing X_m. On the other hand,
the cross coordination number OM-OW increases with a decreasing X_m, the corresponding
hydrogen bonds number follows this behavior. The “inversion” of the composition of the surrounding of
a methanol molecule on average occurs at X_m≈ 0.65. In the region of low methanol mole
fractions, methanol molecule incorporates the hydrogen bonded structure of water and on average forms a bit more than two hydrogen bonds with water molecules.

An excellent detailed analysis of the coordination numbers and hydrogen bonds network 
topology was 
performed in reference <cit.>, though the energetic definition of hydrogen bonds was
used in that work. Nevertheless, all trends of behavior deduced from geometric
criterion for hydrogen bonding qualitatively agree with the predictions described by using
energetic definition.




One of the properties representing a quite demanding test of the employed force field
is the surface tension that describes the changes of surface free energy upon 
changing the surface area. 
Therefore, we have undertaken additional calculations focused
on the evaluation of surface tension for a set of models under study. Relevant
experimental results are available <cit.>. Moreover, this property
has been explored in detail in several computer simulation 
studies of water-methanol system, see, e.g., <cit.>.

The simulations aiming at surface tension calculations at each point of composition
axis in terms of X_m, have been performed by using the final configuration of particles 
in the box from the NPT run.
Next, the box edge along z-axis was elongated by a factor of 3, generating a box with
liquid slab and two liquid-mixture – vacuum interfaces in the x-y plane,
in close similarity to the procedure used in reference <cit.>.
The total number of particles, 3× 10^3, seem to be reasonable as it provides the area of 
the x-y face of the liquid slab large enough to avoid size effects. 
The elongation of the liquid slab along z-axis is sufficient as well.
To be more specific, we would like just to mention that the box dimensions 
were 4.63× 4.63 × 13.89 nm^3 for X_m=0.1 up to 5.75 × 5.75 × 17.25 nm^3 for X_m=0.9,
for the TIP4P/2005-MET/TraPPE model calculations. These numbers favourably compare 
with the cut-off distance for non-bonded interactions 1.4 nm.  
The executable molecular dynamics file was modified by
deleting fixed pressure condition, just the  V-rescale thermostatting with the same
parameters as in the NPT runs has been preserved. Other corrections have not been used.

The values for the surface tension, γ, follow from the combination
of the time averages for the components of the pressure tensor,

    γ = L_z ⟨[P_zz - 1/2(P_xx+P_yy)]⟩/ 2 ,

where P_ij (i,j = x,y,z) are the components of the pressure tensor,
and ⟨ ... ⟩ denotes the time average.
We performed a set of NVT runs, not less than 5–6, each with the time duration of 10 ns,
and obtained the results for γ making the block average.



Our results for the surface tension and the excess mixing surface tension are given
in two panels of figure <ref>. It is known that the TIP4P/2005 model leads to a better 
prediction for the surface tension of pure water in comparison with SPC/E <cit.>.
This is also documented in figure <ref> (a).
Our result for pure methanol within MET/TraPPE model agrees perfectly well with reference <cit.>.
On the other hand, the application of the same procedure with MET/AA model yields a
slightly higher value for pure methanol, very close to the experimental
result, figure <ref> (a). Small discrepancy of our result with the value given in the
supplementary material to <cit.> can be attributed to slightly different 
technical details of simulations, e.g., cut-off, corrections to the
electrostatics for a slab system, etc.

Overall  trends of behavior of γ (X_m) in the entire composition interval
are qualitatively correctly reproduced by all four combinations of the force fields.
A better description at low values of X_m is reached if the TIP4P/2005 model is
invloved. For high values of X_m, the application of the MET/AA model leads to a
better agreement with the experimental data. The values calculated by us are lower 
than the experimental data <cit.> at all X_m, due to the defficiency of performance of 
models for each component, if the TIP4P/2005 and MET/TraPPE are used. With the 
MET/AA model of methanol, the surface tension values are higher than the experimental
points at intermediate and high methanol mole fractions. Interestingly, the surface
tension substantially decreases even if a small amount of methanol is added to water.
This behavior can be attributed to the tendency for methanol to be located close to
the interface with vacuum. Our findings are in agreement with the results reported
in <cit.> for the same system but with a smaller number of particles.


If one focuses on the deviation
of the surface tension from ideal mixing behavior, figure <ref> (b), the TIP4P/2005-MET/TraPPE
model is the closest to the experimental predictions. The present
satisfactory modelling can be seemingly used with confidence to other more complex systems.
It is worth mentioning that experimental results predict a maximum absolute value of
deviation from ideality at X_m≈ 0.3 and computer simulation results
reproduce this behavior. This particular composition actually coincides
with the maximum deviation of mixing entalphy, cf. figure <ref> (b), that has not been
predicted by the models in question. Seemingly, the 
local composition fluctuations with specific orientations of molecules (yielding
nonideality of surface tension) are missing in the bulk phase to provide a 
correct behavior of nonideality of entalphy.
At present, it is difficult to offer a recipe how to improve the performance of nonpolarizable 
rigid models for surface tension and for nonideality of entalphy, though. 



In order to explore how the particles of different species are distributed in the box when
X_m changes, we have plotted the number density profiles of water and methanol 
oxygens in the entire slab plus vacuum system in  figure <ref>. At a low value of X_m,
see X_m=0.1 as an example, water is uniformly distributed in the slab. Methanol species
are distributed uniformly in the inner part of the slab, though we observe a high
maximum for methanol oxygens distribution  at each liquid slab — vacuum interface. 
At intermediate composition, X_m=0.5, both species are almost uniformly distributed
inside the liquid slab, but the maxima of O_m observed for X_m=0.1 disappear, in expence 
of the growing amount of methanol oxygens in the part of each interface exposed to vacuum. 
The two profiles describing the intermediate composition, X_m=0.5, are not perfectly
symmetric with respect to the box center. This situation occurs very rarely
but can have consequences, if one pretends to obtain the precise density of species 
from this kind of procedure. Specifically, in reference <cit.>
the authors failed to estimate the liquid density both for methanol and water at
X_m=0.3.  
Finally, at X_m=0.9 in figure <ref>, it can be seen that the interface predominantly contains methanol molecules 
whereas water molecules “hide” in the inner part of the liquid slab. Only a small
fraction of water molecules penetrates the interfacial region.

The final part of the manuscript is concerned with the description of the 
self-diffusion coefficients of species and changes of the dielectric constant
with composition of water-methanol mixtures.
The self-diffusion coefficients of water and methanol were calculated 
 from the mean-square displacement (MSD) of a particle via Einstein relation,

    D_i =1/6lim_t →∞/ t⟨| 
    r_i(τ+t)- r_i(τ)| ^2⟩ ,

where i refers to water or methanol and τ denotes the time origin. 
Default settings of GROMACS were used for the separation of the time origins. 
Moreover, the  fitting interval (from 10% to 50% of the 
analyzed trajectory) has been used to calculate D_m and D_w.
Moreover, a special care has been taken to fitting for the cases with a small 
number of particles on the extremes along the X_m axis. 
A set of our results is given in  figure <ref>. 



The results for D_w, are qualitatively correct, 
the  TIP4P/2005 model slightly underestimates D_w, it yields 2.1 for pure
water, whereas the SPC/E overestimates D_w, it yields ≈ 2.6, see figure <ref> (a).
Concerning the self-diffusion coefficient of pure methanol, it 
can be seen that MET/TraPPE underestimates D_m, it gives ≈ 2.25, whereas
the all atom model, MET/AA, substantially overestimates D_m, figure <ref> (b). This tendency
has been already discussed in reference <cit.>. These inaccuracies prohibit
quantitatively correct predictions for the composition changes of the self-diffusion
coefficients of species. 

Namely, the dependence D_w(X_m) is qualitatively correct. The TIP4P/2005, if combined
with MET/TraPPE or with MET/AA methanol model, leads to a very good prediction of D_w
up to X_m≈ 0.3, for higher values of X_m simulation data deviate from 
experimental results <cit.>. 
The miminum value for D_w from simulations is in the interval 
X_m between 0.4 and 0.5 whereas the experiment predicts this minimum at 0.3.
The minimum value of D_w from simulations along the composition axis coincides with the minimum of 
excess mixing volume at X_m=0.5. At high values of X_m, the growth of D_w is similar 
in simulations and in experiment <cit.>. 
If water is descibed in the framework of the SPC/E
model, the resuts for D_w(X_m) substantially overestimate this self-diffusion 
coefficient in a wide range of composition. Solely at high values of X_m, the agreement
with experimental values becomes more acceptable. It is worth mentioning that we
were unable to get a better shape of D_w(X_m) within the SPC/E-MET/TraPPE model. 
Apparently, there are two minima, one at X_m=0.4 and anoter at X_m=0.7. It is
difficult to establish (without additional exploration of various properties) if
this behavior is related to the clustering of species at local scale as it has
been discussed in the experimental study <cit.>. Here, we would like
just to mention that similar evolution of D_w on X_m was reported in
our recent work <cit.> in the study of SPC/E model combined by OPLS/UA model 
of methanol <cit.>.  

Concerning the trends of behavior of the self-diffusion of methanol species in mixtures
of different composition D_m(X_m), we would like to emphasize the following. 
The best combination of force fields is provided by the TIP4P/2005-MET/TraPPE model.
It describes the changes of the function D_m(X_m) pretty well in the entire composition 
range, figure <ref> (b). The minimum value of D_m is described at a correct place, X_m≈ 0.3.
All other combinations of the force fields exhibit defficiencies either due to the water model
like SPC/E or due to all-atom modelling of methanol at MET/AA level.  
An overall most satisfactory picture of the dependence of the self-diffusion coefficients 
of two species, therefore, results if the TIP4P/2005-MET/TraPPE model is used.
We believe that alternative calculations of the self-diffusion coefficients by 
applying velocity autocorrelation functions should lead to similar conclusions.

Our final concern is the evolution of the dielectric constant with composition.
The long-range, asymptotic behavior of correlations between molecules possessing a
permanent dipole moment is determined by the dielectric 
constant, ε.
Usually,  long molecular dynamics runs are necessary to obtain
reasonable values for ε,
because it is calculated from the time-average of the fluctuations of the total
dipole moment of the system <cit.>,

    ε=1+4/3k_BTV(⟨ M^2⟩-⟨ M⟩^2),

where k_B is the Boltzmann constant and V is the volume of the simulation box.



The lines from our NPT simulations of ε are shown in figure <ref> (a).
An overall behaviour of ε(X_m) is that it decreases
starting from a high value for pure water
to a lower value corresponding to pure methanol with an increasing X_m.
As it follows from the comparison of the simulation results and experimental 
data <cit.>,
all four combinations of the force fields substantially underestimate the values 
for ε in the entire composition range. The reason is that the
dielectric constant for two constituents, water and methanol, is
essentially underestimated.
It is highly probable to improve the dependence of the static dielectric constant
on composition by applying the model for water specifically parametrized to reproduce
the dielectric constant for water <cit.>. Simultaneously, it would require
parametrization of the force field for methanol, see, e.g., <cit.>.
However, common experience is that  parametrization of a single property leads to worse 
predictions for other properties. Therefore, this issue requires additional
computation efforts. 

On the other hand, a sensitive test is provided by comparison of the excess
dielectric constant,
Δε^mix = ε_m- [X_mε_m+(1-X_m)ε_w],
with the experimental predictions <cit.>.
Experimental points indicate a negative deviation from
ideality in the entire composition range, figure <ref> (b).
Maximal (negative) deviation from the ideal type of behaviour
reported from the experimental measurements is at X_m≈  0.45.
The simulations results  reproduce the position of a minimum approximately; namely,
the minimum is in the interval 0.15–0.4, dependent on the combination of the force fields.
These trends are in accordance with observations concerning the deviations from ideality
of thermodynamic properties.
The magnitude of the excess static dielectric constant
is overestimated if the SPC/E model for water is used and is
underestimated if TIP4P/2005 is involved.
The TIP4P/2005-MET/TraPPE model is the most close to the experimental data combination of models. The excess dielectric
constant curve as function of chemical composition of the mixture
can be related to the excess refractive index measurements, see, e.g., <cit.>.
Consequently, a complementary comparisons with experiments would be desirable in future
work.






§ SUMMARY AND CONCLUSIONS


This work has been principally inspired by the necessity to evaluate
various properties of water-methanol mixtures using a specific set of 
models for each component. The choice of the force fields of this study 
permits extensions necessary to describe ionic solutions with such combined solvents
in the spirit of very recent contributions from the research
laboratory of C. Vega <cit.> concerning 
aqueous NaCl solutions with a novel force field for ions.

In the present work, extensive NPT  molecular dynamics simulations were conducted to study
thermodynamic, dynamic and structural properties of water-methanol
mixtures. Two different water and two methanol nonpolarizable models 
were combined and simulated with the purpose of testing their predictions
for an ample set of properties in the entire range of compositions. 
Comparisons with the available experimental data were performed.
Considering the scope of the models, the predictions obtained for
the mixtures appear to be qualitatively correct, particular properties
were a bit better described at low methanol compositions, i.e., for
water-rich compositions. However, as a general trend, it is observed
that the best predictions are given with the
water (TIP4P/2005)-methanol(TraPPE) mixtures since these models reasonably well predict several properties of the pure components.
The results indicate that a good agreement with laboratory experiments
could be obtained when both force fields, of the two components in the mixture,
are good. In fact, the dielectric constant is not well predicted by any of the
simulated mixtures since none of the selected models (water and methanol)
predicts correctly that property.





§ ACKNOWLEDGEMENTS

M. Cruz and O.P. are grateful to M. Aguilar for technical support of this work
at the Institute of Chemistry of the UNAM. M. Cruz acknowledges support of CONACyT 
of Mexico for Ph.D. scholarship.




   


99



wensink  Wensink E.J.W., Hoffmann A.C., van Maaren P.J., van der Spoel D., J. Chem. Phys., 2003, 119, 7308, 
 10.1063/1.1607918.

galicia1 Galicia-Andrés E., Pusztai L., Temleitner L., Pizio O.,
J. Mol. Liq., 2015, 209, 586, 
 10.1016/j.molliq.2015.06.045.

galicia2 Galicia-Andrés E., Dominguez H., Pusztai L., Pizio O., Condens. Matter Phys., 2015, 18, 43602, 
 10.5488/CMP.18.43602.


perera1  Perera A., Sokolić F., Almásy L., Koga Y., J. Chem. Phys., 2006, 124, 124515, 10.1063/1.2178787.

perera2 Perera A., Zoranić L., Sokolić F., Mazighi R., J. Mol. Liq., 2011, 150, 52,
10.1016/j.molliq.2010.05.006.

bako Bakó I., Pusztai L., Temleitner L., Sci. Rep., 2017, 7, 1073, 10.1038/s41598-017-01095-7.

bopp1 Palinkás G., Bakó I., Heinzinger K., Bopp P., Mol. Phys., 1991, 73, 897, 10.1080/00268979100101641.

bopp2 Soetens J.-C., Bopp P.A., J. Phys. Chem. B, 2015, 119, 8593, 
10.1021/acs.jpcb.5b03344.

jordan Obeidat A., Abu-Ghazleh H., AIP Adv., 2018, 8, 065203,
10.1063/1.5025575.

matsumoto Matsumoto M., Takaoka Y., Kataoka Y., J. Chem. Phys., 1993, 98, 1464,
10.1063/1.464310.

chang Chang T.-M., Dang L.X., J. Phys. Chem. B, 2005, 109, 5759,
10.1021/jp045649v.

takamuku Takamuku T., Yamaguchi T., Asato M., Matsumoto M., Nishi N.,
Z. Naturforsch., A: Phys. Sci., 2000, 55, 513, 

10.1515/zna-2000-0507.

wakisaka Wakisaka A., Ohki T., Faraday Discuss., 2005, 129, 231,
10.1039/B405391E.

mikhail Mikhail S.Z., Kimel W.R., J. Chem. Eng. Data, 1961, 6, 533, 10.1021/je60011a015.

mcglashan McGlashan M.L., Williamson A.G., J. Chem. Eng. Data, 1976, 21, 196, 
10.1021/je60069a019.

lama Lama R.F., Lu B.C.-Y., J. Chem. Eng. Data, 1965, 10, 216, 
10.1021/je60026a003.

vazquez Vazquez G., Alvarez E., Navaza J.M., J. Chem. Eng. Data, 1995, 40, 611,
10.1021/je00019a016.

trappe Chen B., Potoff J.J., Siepmann J.I., J. Phys. Chem. B, 2001, 105, 3093, 10.1021/jp003882x.

jorgensen2 Jorgensen W.L., Maxwell D.S., Tirado-Rives J., J. Am. Chem. Soc., 1996,
118, 11225, 10.1021/ja9621760.

spce Berendsen H.J.C., Grigera J.R., Straatsma T.P., J. Phys. Chem., 1987, 
91, 6269, 10.1021/j100308a038.

jorgensen3  Jorgensen W.L., Chandrasekhar J., Madura J.D., Impey R.W., Klein M.L., J. Chem. Phys., 1983, 79,  926, 

10.1063/1.445869.

vega-tip Abascal J.L.F., Vega C., J. Chem. Phys., 2005, 123, 234505, 10.1063/1.2121687.

pusztai Pusztai L., Pizio O., Sokolowski S., J. Chem. Phys., 2008, 129, 184103, 10.1063/1.2976578.

vega-pccp Vega C., Abascal J.L.F., Phys. Chem. Chem. Phys., 2011, 13, 19663, 10.1039/C1CP22168J.

benavides1 Benavides A.L., Aragones J.L., Vega C., J. Chem. Phys., 2016, 144, 124504, 10.1063/1.4943780.

madrid Benavides A.L., Portillo M.A., Chamorro V.C., Espinosa J.R.,
Abascal J.L.F., Vega C., J. Chem. Phys., 2017, 147, 104501, 10.1063/1.5001190.

bagchi  Hazra M.K., Roy S., Bagchi B., J. Chem. Phys., 2014, 141, 18C501, 10.1063/1.4895539.

roccatano Samanta S., Roccatano D., J. Phys. Chem. B, 2013, 117, 3250, 10.1021/jp309476u.

patsahan Patsahan T., Ilnytskyi J.M., Pizio O., Condens. Matter Phys., 2017, 20, 23003, 10.5488/CMP.20.23003.

horn Horn H.W., Swope W.C., Pitera J.W., Madura J.D., Dick T.J., Hura G.L., Head-Gordon T., J. Chem. Phys., 2004, 120, 9665,
10.1063/1.1683075.

jorgensen Jorgensen W.L., J. Phys. Chem., 1986, 90, 1276,
10.1021/j100398a015.

haughney Haughney M., Ferrario M., McDonald I.R., J. Phys. Chem., 1987, 91, 4934,
10.1021/j100303a011.

guevara Guevara-Carrion G., Nieto-Draghi C., Vrabec J., Hasse H., J. Phys. Chem. B, 2008, 112, 16664, 
 10.1021/jp805584d.

vrabec2 Schnabel T., Srivastava A., Vrabec J., Hasse H., J. Phys. Chem. B, 2007, 111, 9871, 10.1021/jp0720338.

hasse1 Kohns M., Horsch M., Hasse H., Fluid Phase Equilib., 2018, 458, 30, 10.1016/j.fluid.2017.10.034.

pozar Požar M., Kerasidou A., Lovrinčević B., Zoranić L., Mijaković M., Primorac  T., 
Sokolić F., Teboul V., Perera A., J. Chem. Phys., 2016, 145, 144502, 10.1063/1.4964487.





gromacs Van der Spoel D., Lindahl E., Hess B., Groenhof B., Mark A.E.,
Berendsen H.J.C., J. Comput. Chem., 2005, 26, 1701, 10.1002/jcc.20291.

washbrun Washbrun E.W. (Ed.), International Critical Tables of Numerical Data, Physics, Chemistry and Technology, Knovel, New York, 2003.

bopp Bopp P., Jancsó G., Heinzinger K., Chem. Phys. Lett., 1983, 98, 129, 10.1016/0009-2614(83)87112-7.

palinkas Palinkas G., Hawlicka E., Heinzinger K., J. Phys. Chem., 1987,
91, 4334, 10.1021/j100300a026.

nezbeda González-Salgado D., Nezbeda I., Fluid Phase Equilib., 2006, 240, 161, 10.1016/j.fluid.2005.12.007.

old-palinkas Pálinkás G., Bakó I., Z. Naturforsch., A: Phys. Sci., 1991, 46, 95, 10.1515/zna-1991-1-215.

galicia3 Galicia-Andrés E., Dominguez H., Pusztai L., Pizio O., J. Mol. Liq., 2015, 212, 70, 

10.1016/j.molliq.2015.08.061.

chihaia Bakó I., Megyes T., Bálint S., Grósz T., Chihaia V., Phys. Chem. Chem. Phys., 2008, 10, 5004, 
 10.1039/b808326f.

fischer Fischer N.M., van Maaren P.J., Ditz J.C., Yildrim A.,
van der Spoel D., J. Chem. Theory Comput., 2015, 11, 2938, 
10.1021/acs.jctc.5b00190.

biscay Biscay F., Ghoufi A., Malfreyt P., J. Chem. Phys., 2011, 134, 044709, 
10.1063/1.3544926. 

derlacki Derlacki Z.J., Easteal A.J., Edge A.V.J., Woolf L.A., Roksandic Z., J. Phys. Chem., 1985, 89, 5318, 
 10.1021/j100270a039.

martin Neumann M., Mol. Phys., 1983, 50, 841, 10.1080/00268978300102721.

albright Albright P.S., Gosting L.J., J. Am. Chem. Soc., 1946, 68, 1061, 10.1021/ja01210a043.

alejandre Alejandre J., Chapela G.A., Saint-Martin H., Mendoza N., 
Phys. Chem. Chem. Phys., 2011, 13, 19728, 
 10.1039/C1CP20858F.

dominguez Salas F.J., Méndez-Maldonado G.A., Núñez-Rojas E., Aguilar-Pineda G.E., Domínguez H., Alejandre J., J. Chem. Theory Comput., 2015, 11, 683, 10.1021/ct500853q.

gofurov Gofurov Sh., Ismailova O., Makhmanov U., Kokhkharov A., Int. J. Chem. Mol. Nucl. Mater. Metall. Eng., 2017, 11, 330.





































































 




	label1  쳿,  i  i . i, i, 
	label2  ,  i  i . i, i, 





	     -   
	          -.
	   TIP4P/2005  SPC/E     '   TraPPE        . 
	              
	    . ,      ,
	  '   ,         .
	 ,            '
	   . ,      糿
	.         .
	
	,  ,  ,    




















