Revisiting the composition dependence of the properties of water-dimethyl sulfoxide liquid mixtures. Molecular dynamics computer simulations

We have revisited the composition dependence of principal properties of liquid water-DMSO mixtures by using the isobaric-isothermal molecular dynamics computer simulations. A set of non-polarizable semi-flexible models for the DMSO molecule combined with the TIP4P-2005 and TIP4P/$\varepsilon$ water models is considered. We restrict calculations to atmospheric pressure, 0.1013 MPa, and room temperature, 298.15 K. Composition trends of density, excess mixing volume and excess mixing enthalpy, partial molar volumes and partial molar enthalpies of species, apparent molar volumes are reported. Besides, we explore composition trends of the self-diffusion of species, the static dielectric constant and the surface tension. Evolution of the microscopic structure of the mixture with composition is analyzed in terms of radial distributions functions, coordination numbers and the fractions of hydrogen-bonded molecules. We intend to capture the peculiarities of mixing the species in the mixture upon the DMSO molar fraction and the anomalous behaviors, if manifested in each of the properties under study. The quality of several combinations of the models for species is evaluated in detail to establish the possibility of necessary improvements.


Introduction
Liquid mixtures of dimethyl sulfoxide (DMSO) and water are of much practical importance as solvents and reaction media in organic chemistry and chemical engineering [1][2][3]. Besides, they are widely used in several biological applications, cryoprotection being just one example, see e.g., [4][5][6].
Thus, it is not only desirable but also necessary to have an accurate description of various properties of these mixtures dependent on composition, , temperature, , and pressure, . Common experimental practice is to study the composition effects at a room temperature and not far from it. Moreover, in many studies the atmospheric pressure is assumed. Much less studies are focused on temperature and pressure trends of the behavior of water-DMSO liquid mixtures. Unfortunately, experimental results do not cover all the cases of interest. In other words they do not provide comprehensive knowledge about the dependence of even some basic properties of the system in question, upon changes of these thermodynamic variables. As documented in reference [7], the experimental knowledge and understanding of the microscopic structure and dynamic properties of the systems in question mainly follows from the application of neutron scattering, nuclear magnetic resonance, dielectric relaxation and vibrational spectroscopy methods. In order to interpret the experimental observations in detail and to get ampler insights, one is forced to resort to computer simulations methodology. all atoms and/or groups is assumed as a sum of Lennard-Jones (LJ) and Coulomb terms. A comprehensive description of the parameters of models of this type is given in table 1 of [40] and in table 1 of [13]. In all cases, the models are considered as rigid. Some basic properties under specific conditions concerned with computer simulations setup are also reported in table 2 of [40]. A set of parameters for each model follows from fitting to a few target properties, e.g., the heat of vaporization and charge distribution in order to yield a dipole moment of the molecule in the gas phase, which is the case to construct for example the P2 model for DMSO [10]. Experimental predictions for the intra-molecular structure and auxiliary ab initio calculations were used in the design of DMSO models as well. Several attempts to better describe a wider set of target properties of DMSO at ambient conditions were undertaken via parametrization of charges and LJ parameters still keeping the model rigid [11,15,32]. In order to improve the rigid P2 model, Benjamin proposed to include the effects of intra-molecular flexibility [28]. The parameters for intra-molecular potential describing harmonic bonds and angles were chosen to reproduce the gas phase normal mode frequencies. It was shown however [29], that the model essentially overestimates the liquid DMSO density at room temperature. Thus, it is improbable to expect a good level description of the density of water-DMSO mixtures over the entire composition range.
Both rigid or flexible DMSO model can also be constructed by using a well established OPLS database [41]. The OPLS and P1 and P2 models [10] are characterized by identical molecular geometry, though the LJ parameters are different. The site charges are the same for the OPLS and P2 models. In the construction of a flexible model, one can take parameters from the OPLS database. However, to keep the molecular geometry, it is then recommended to complement the DMSO model by the improper dihedral angle S-C-O-C see reference [42].
In the present study we apply the DMSO united atom models with four sites, O, S, CH 3 . The CH 3 is abbreviated as C in what follows. The models are non-polarizable and denominated as P1, P2 [10] and OPLS (see e.g., table II of reference [43]). The interaction potential between all atoms and/or groups is a sum of Lennard-Jones (LJ) and Coulomb terms. All the parameters are given in table 1 for the sake of convenience of the reader. In addition, after testing several possibilities, we assumed C-S and S-O bonds rigid whereas the O-S-C and C-S-C are chosen flexible, the harmonic constant is taken from the OPLS database. Moreover, the improper dihedral S-C-O-C is included into the DMSO models. The parameters are taken from reference [42]. Our tests indicate that this element is crucial to yield a reasonable value for pure DMSO liquid density in each of the models involved in this study. In summary, we applied semi-flexible models for DMSO molecule. For water, the TIP4P-2005 model was applied [44]. In many cases we used the TIP4P/ model [45] as well. It was designed to improve the inadequacy of the prediction of the static dielectric constant of water by TIP4P-2005. Lorentz-Berthelot combination rules were used to determine the cross parameters for the relevant potential well depths and diameters for the P1-water and P2-water models, whereas the geometric combination rule was assumed for the OPLS-TIP4P-2005 model.
Molecular dynamics computer simulations of DMSO-water mixtures were performed in the isothermalisobaric (NPT) ensemble at atmospheric pressure 1 bar and at temperature 298.15 K. We used GROMACS package [46] version 5.1.2. The simulation box in each run was cubic, the total number of molecules of both species is fixed at 3000. Composition of the mixture is described by the molar fraction of DMSO molecules, = /( + ). As common, periodic boundary conditions were used. Temperature and pressure control was provided by the V-rescale thermostat and Parrinello-Rahman barostat with 33202-3 = 0.5 ps and = 2.0 ps, the timestep was 0.001 ps. The value of 4.5 × 10 −5 bar −1 was used for the compressibility of mixtures.
The non-bonded interactions were cut-off at 1.1 nm, whereas 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 the precision 10 −5 . The van der Waals correction terms to the energy and pressure were applied. In order to maintain the geometry of water molecules and DMSO intra-molecular bonds rigid, the LINCS algorithm was used.
After preprocessing and equilibration, consecutive simulation runs, each for not less than 10 ns, with the starting configuration being the last configuration from the previous run, were performed to obtain trajectories for the data analysis. The results for each property were obtained by averaging over 7-10 production runs.

Density of water-DMSO mixtures on composition
As we mentioned in the introductory section, there were several experimental reports concerning the evolution of density of DMSO-water mixtures upon changing composition. We used experimental data by Egorov and Makarov at room temperature = 298.15 K, and atmospheric pressure [47]. The experimental data at = 298.15 K witness that the density of water-DMSO mixtures increases quite rapidly with starting from pure water density (997 kg/m 3 ) and reaches a weakly pronounced maximum around ≈ 0.6. Next, with a further increase of , i.e., in the interval where the DMSO amount within the mixture predominates, the density slightly decreases with and reaches the pure liquid DMSO density, = 1095 kg/m 3 , figure 1 (panel a). It is reasonable to interpret the presence of a maximum density along the composition axis as an evidence of the stronger water-DMSO interaction, compared to the water-water inter-molecular interaction [47].
Molecular dynamics computer simulation results for P1, P2 and OPLS DMSO models combined with TIP4P-2005 water model in the entire composition range are given in figure 1 a by red, blue and magenta lines, respectively. We observe that qualitative trends of behavior of density on composition are correctly  reproduced by the P1-TIP4P-2005 and P2-TIP4P-2005 models. The maximum of density from simulations of these models is captured. Moreover, the growth and decay of density in the interval of compositions with small DMSO amount and in the interval of small water concentrations, respectively, are reproduced by both models as well. By contrast, the OPLS-TIP4P-2005 model predicts a monotonous growth of the mixture density with an increasing in the entire composition interval. Thus, the simulations of all model combinations show that the density increases upon adding foreign species to water. If water is added to pure DMSO solvent, P1 and P2 predict a growth of density. The OPLS model leads to an opposite trend, it yields a decreasing density, or in other words at constant pressure the volume of principal simulation box increases. If water enters DMSO structure within P2 or P1 models, the attraction dominates, then the density increases. On the other hand, if water enters the "sea" of OPLS molecules, repulsion prevails and consequently the box increases or equivalently the density decreases. The

Excess mixing volume, partial molar and apparent molar volumes of species on composition
However, it is of importance not only to describe the trends of behavior of the given property on composition, but to correctly capture the deviation from ideality as well. These insights follow from the behavior of the excess density or the excess mixing volume, for example. The excess mixing volume is defined as follows, Δ mix = mix − − (1 − ) , where mix , and refer to the molar volume of the mixture and of the individual components, DMSO and water, respectively.
In order to elucidate the composition trends of the behavior of the excess mixing volume, we plot the simulation results and experimental data at 298.15 K in figure 2. Experimental data show that Δ mix is negative and exhibit a minimum at The maximum deviation from the ideal mixing at ≈ 0.4 is interpreted as trends for the formation of "complexes" or "associates" of the 2DMSO-3H 2 O type [18]. On the other hand, it is documented that the maximum deviation from ideality of density and of other properties on composition can be characterized in terms of percolation phenomena [48].
In order to discern the contribution of each species of the mixture into the mixing volume, it is common to resort to the excess partial molar volumes. They are defined as, = − , where denotes the partial molar volume of species ( or ). The partial molar volumes result from the excess mixing volume, Δ mix as follows [49], where = 1 − . The excess partial molar volumes of DMSO and water species are plotted in the panels a and b of figure 3. The experimental data are taken from figures 5 and 6 of reference [49]. The best agreement with the experimental results for is obtained for the P2-TIP4P-2005 model. Namely, the shape of the curve and the minimum for are reproduced quite well. Predictions from the P1-TIP4P-2005 model are reasonable as well. However, the location and the value for the minimum are less accurate in comparison to P2-TIP4P-2005 model. Concerning the predictions of different models w.r.t. the , the situation is a bit worse. Both models, the P2-TIP4P-2005 and P1-TIP4P-2005, describe the behavior of the excess partial molar volume well at a high water content. On the other hand, for DMSO-rich mixtures, the simulation data deviate from the experimental results. The experimental data show that the minimum of is weakly pronounced and is located at ≈ 0.1. By contrast, the simulation data predict much stronger pronounced minimum at a higher water content. In summary, the dissolution of water species in the DMSO solvent is not accurately described. The worse predictions for and come out from the application of the OPLS-TIP4P-2005 model. Similar insights into the geometric aspects of mixing on composition, both from experiments and simulations, can be obtained by resorting to the notion of the apparent molar volume of species rather than the partial molar volumes. The apparent molar volume for each species [49] according to the definition is: ( ) = + Δ mix /(1 − ) and ( ) = + Δ mix / . We elaborated the experimental density data from references [47,50] and the results from our simulations to construct the plots shown in panels c and d of figure 3. These plots confirm that the P2-TIP4P-2005 model provides a reasonable description of the composition behavior for ( ) in water-rich mixtures while the P1-TIP4P-2005 model is satisfactory for DMSO-rich mixtures.
This kind of composition behavior of apparent molar volumes of DMSO species in water-rich mixtures can be related to the experimental results for abnormal intensity of scattered light [51]. The minimum of  the apparent molar volume of DMSO indicates the hydrophobic effect in the systems in question due to enhanced concentration fluctuations within certain composition interval.

Energetic aspects of mixing of DMSO and water molecules
Next, we proceed to the energetic manifestation of mixing trends. The experimental data are scarce in this aspect, we used the results from references [52][53][54][55] to explore the mixing enthalpy upon composition of the water-DMSO mixtures.
To begin with, the simulation results for the enthalpy on composition for different combinations of models are given in panel a of , that the P1 and P2 DMSO models combined with the TIP4P-2005 water substantially over-predict the mixing enthalpy. Changing the water model from TIP4P-2005 to TIP4P/ does not yield an improvement in the considered aspect. This seems to be a quite general failure of a set of non-polarizable water and DMSO combinations of models as documented in the comprehensive study [13], see figures 3 and 5 of that reference. On the other hand, from our figure 4 b, we observe that the OPLS-TIP4P-2005 model underestimates the mixing enthalpy but the results are not far from the experimental points. The excess partial molar enthalpies can be obtained by using equations (3.1) and (3.2), just the volume should be substituted by enthalpy. The agreement between the experimental data and predictions of a set of models in question can be termed as qualitative and in some cases as quantitative for water-rich or DMSO-rich mixtures. The OPLS-TIP4P-2005 model yields the most reasonable coincidence with the experimental data. However, the energetic aspects of mixing of species, in terms of and simultaneously, are better described for water-rich mixtures.

Comments on the evolution of the microscopic structure upon changes of composition
In all the above discussion concerned with geometric and energetic aspect of mixing of water and DMSO species from molecular dynamics computer simulations of different models, we have intended to validate the theoretical predictions w.r.t. experimental results in every detail. One of well established observations is in the existence of the minimum of the excess partial molar volume or of the apparent molar volume of DMSO species at a low . This specific interval of composition is the subject of search of anomalous behavior of different properties in references [20,22,34]. Some insights in those works were obtained from the analysis of the microscopic structure.
Without referring directly to experimental results, it is difficult to prove that the computer simulation predictions for the microscopic structure of the mixtures in question are reasonable. Recently, the contributions of our collaborators [56,57] have scrutinized this problem for pure water and water-methanol mixtures. Unfortunately, we are not aware of comprehensive experimental measurements of the total and partial structure factors in the entire composition interval for water-DMSO mixtures. Thus, the discussion below should be considered as mirror interpretation of thermodynamic results given above, rather than as independent and confident source to capture possible anomalous behaviors.
Consequently, we would like to perform mere analyses of the structural motifs by considering how the partial radial distribution functions change upon adding the DMSO molecules to water. The results are presented in the following figures. We restrict principally to the results for P2-TIP4P-2005 model because it provides the best description of the density dependence on composition and of other properties of water-rich mixtures, as documented above. We observe that some of the functions exhibit anomalous changes whereas the other functions change monotonously. Specifically, the radial distribution function for water oxygens in given in figure 5. The first maximum of this function grows in magnitude with an increasing DMSO concentration. This growth is accompanied by a decrease of the height of the second maximum. Moreover, the location of the first maximum along -axis does not change with , whereas the following oscillations slightly shift to larger distances. These changes can be characterized by the first (inset to panel a) and the second (panel b) coordination numbers evaluated by integration to the respective minimum of the ( ) as common, where is the density of species . This kind of evolution of ( ) for water species in water-DMSO mixtures was discussed already in reference [10] (figure 7 of [10]) and in reference [33]. In addition to previous observations, we would like to note that the first coordination number exhibits a weak nonlinearity as shown in the inset to figure 5. Besides, we observed that the second coordination number changes its shape on at ≈ 0.1667 as a reflection of changes of the second maximum of ( ) for water oxygens. A narrow range of compositions around this value is analyzed in [20] suspecting an anomalous behavior of the number of contacts between methyl groups of DMSO molecules in water-rich mixtures, cf. figures 8 and 9 of reference [20]. In summary to figure 5, we note that the behavior of ( ) for OW-OW can be interpreted as the enhancement of the short-range structure of water upon adding the DMSO species. On the other hand, trends for a weaker structure at larger separations between water molecules witness that on a larger scale, hydrogen bonded network of water becomes dismantled. Nevertheless, the entire interval of composition considered in figures 5 and 6 refers to systems with a rather high degree of hydrogen bonding. The average number of hydrogen bonds (per water molecule) between species of the mixture is illustrated in the inset to figure 6. We used GROMACS software with default parameters to obtain the data. The average number of H-bonds between water molecule drops down from 3.7 to 2.5 within the interval considered. On the other hand, the bonds are formed upon DMSO insertion into water medium, such that their sum, or say bonding state of water molecules, remains almost as high as in pure water. We have not elaborated the H-bonding structure of the systems in question more in detail, however. Certain additional details can be deduced by using a previous observation for our laboratory, see e.g., figure 8 and 9 of reference [33]. It is worth to extend the study of these issues by involving the energetic definition of hydrogen bonding, see e.g., [58], and by searching for the elements of microscopic structure due to cooperativity of hydrogen bonds in the spirit of reference [59]. The hydrogen bonds lifetime should be explored as well.
Some aspects of the structural arrangement of different species in water-DMSO mixtures upon increasing the DMSO concentrations are illustrated in figure 6. The principal effect is shown in panel a of this figure. Namely, the probability of finding water hydrogens close to the oxygens of DMSO molecules substantially increases with an increasing . The second maximum increases as well. The first maximum is observed at the same OD-HW distance upon increasing . Thus, this structure reflects the formation of hydrogen bonds between DMSO and water molecules. The second maximum, presumably corresponds to the second shell of water molecules bonded to molecules of the first shell around the "reference" oxygen of DMSO molecule. At larger distances, the distribution function remains almost flat. A quite similar shape of this function was documented previously at = 0.35, cf. figure 3 of reference [10]. Because of bonding discussed just above, the distribution function of water oxygens w.r.t. methyl groups exhibits similar trends, panel b of figure 6. Its first maximum increases with increasing . However, the second maximum remains constant and then decreases in magnitude when the DMSO concentration increases. Note that both maxima are observed at the distances larger than the first and the second maximum of OD-HW distribution in panel a, as expected.
The pair distribution functions describing the changes of distribution of methyl groups, C-C (the inter-molecular part) and of OD-OD behave differently, in comparison to our discussion above. They are plotted in figure 7. Namely, the first maximum of the radial distributions shown in panels a and b decreases in magnitude when increases from 0.0333 up to 0.1333. Next, the first maximum of both functions grows, if increases further, up to 0.2333. The magnitude of the second maximum exhibits similar trends. The shape of the functions from our calculations favourably compares to the ones obtained in [10] by using the P1 and P2 models combined with the SPC water model, cf. figures 1 and 4 of that reference.
However, the radial distribution functions of methyl groups require additional comments because it should reflect possible aggregation of hydrophobic methyl groups. In order to get more detailed insights into the behavior of this function, we  [20] from the analyses of the total fraction of contacts between methyl groups on composition and of the average number of methyl groups around a given group (see figures 9 and 10 in reference [20]). However, these authors claimed to observe a discontinuous transition termed as weak transition between two regimes, in contrast to our observations. A similar curve coming from the calculations for P2-TIP4P-2005 is given in figure 8 b as well. In this case, there is a very weakly pronounced change of the inclination of the curve between = 0.1667 and 0.2. However, it can hardly be interpreted as a pronounced discontinuity.

33202-11
We would like to finish this part by recalling that the predictions of the OPLS-TIP4P-2005 model for density, for the excess mixing volume and related partials as well as for energetic aspects of mixing of species are not satisfactory. Therefore, it does not seem consistent to discuss anomalous behaviors by using solely certain aspects of microscopic structure from a particular model force field that leads to a quite inaccurate thermodynamic and other properties. Evidently, our comments concerning the microscopic structure are not comprehensive. It would be desirable to validate the structure by comparison with the experimental predictions for the structure factors in the entire composition range. Unfortunately, as it was mentioned at the beginning of this subsection, such data are not available in the literature, up to our best knowledge.
In contrast to microscopic structure, there are different possibilities to validate the force field models against the experimental data. One of the possibilities is to discuss the dynamic properties. Most common, computer simulations predictions are confronted with the experimental results for self-diffusion coefficients of species upon the composition of the solutions in question.

Self-diffusion coefficients of DMSO and water molecules
We proceed now to the description of the behavior of the self-diffusion coefficients of species, ( = , DMSO), upon changing composition of the mixture. In figure 9 we show the results obtained via the Einstein relation, where denotes the time origin. Default settings of GROMACS were used for the separation of the time origins and for the fitting interval. The experimental data were taken from references [60,61]. According to the experiments, the self-diffusion coefficient of water species decreases in magnitude starting from pure water value (at = 0) and reaches minimum at ≈ 0.3. Next, for higher , does not change much, but a certain pecuiliarity is seen in the interval between 60% and 80%. Apparently, the behavior of ( ) is determined by the evolution of density of the mixture and by the hydrogen bonding between species. Theoretical predictions for lead to a bit inaccurate, but still satisfactory agreement with experimental trends for water-rich mixtures. The OPLS-TIP4P-2005 model looks the best, because it underestimates the density in the entire composition interval starting from water-rich and ending up at DMSO-rich mixtures, cf. figure 1. This model yields a good value for close to = 1, as well.

33202-12
Both, the P1-TIP4P-2005 and P2-TIP4P-2005 model yield qualitatively correct trends but essentially underestimate in a wide interval of compositions. Besides two factors mentioned just above, we can attribute such low values for to the pronounced inaccuracy of description of the energetic factors of mixing of species in the framework of these two models, cf. figure 4 b. Better performance of the OPLS-TIP4P-2005 model for the excess mixing enthalpy apparently contributes to better description of . According to the experimental data by Packer and Tomlinson [60], the self-diffusion coefficient of DMSO species has a minimum at ≈ 0.7. Besides, there is a peculiarity of DMSO in the interval of compositions corresponding to water-rich mixtures. The model of the present study yields satisfactory predictions for DMSO close to ≈ 0. Each of the models in question exhibits a weakly pronounced peculiarity around ≈ 0.1. However, this peculiarity is much weaker compared to what one observes from experimental points. On the other hand, close to ≈ 1, the P1-TIP4P-2005 model describes DMSO pretty well. However, the minimum predicted by the models occurs at intermediate compositions around 0.5, rather than at 0.7 from the experimental points. To summarize, the models used in simulation perform better at "wings", i.e., close to ≈ 0 and ≈ 1, rather than at intermediate compositions. Still, the dynamic aspects of mixing are apparently better described for water-rich mixtures. We have not performed a detailed investigation of the selfdiffusion coefficient of species for DMSO models combined with TIP4P/ model. It is known that the diffusion coefficient of pure water in the framework of TIP4P-2005 and TIP4P/ model is very similar [45]. Thus, we do not expect to obtain an improvement of the observed trends by changes of the water models of this type.
The missing elements worth to explore are various. Namely, it is necessary to confirm the present results by using alternative calculations of the self-diffusion coefficients via the velocity auto-correlation functions. In addition, one should attempt to calculate the relaxation times and possibly the power spectra, similarly to e.g., [62], because of the availability of experimental data. It would provide ampler insights into the adequacy of the force fields.

Static dielectric constant of water-DMSO mixtures
Now, proceed to the results for the static dielectric constant. It is one of the most difficult properties to deal with. Recently, we explored the dielectric constant for water-DMSO mixtures by using the SPC-E water model [33]. Therefore, the present discussion will be brief and solely novel findings will be underlined.

33202-13
The long-range, asymptotic behavior of correlations between the molecules possessing a dipole moment is described by the dielectric constant, , which is calculated from the time-average of the fluctuations of the total dipole moment of the system [63] as follows, where B is the Boltzmann constant and is the simulation cell volume. In this work we took care that each of the runs is of sufficient time extension (more than 100 ns in all cases), in order to obtain reasonable statistics. The final value was recorded. In contrast to reference [33], the averaging over blocks was not performed. Besides the stability of the dielectric constant value along the run, we verified the temporal behavior of the average total polarization of the system to ensure that each of the components ( , and ) is close to zero. Evidently, the number of molecules, type of thermostat and barostat, precision of the long-range interactions summation contribute into the accuracy of the evaluation of the dielectric constant. The experimental data are taken from reference [64]. The dielectric constant decreases from a high value for water to a much lower value for pure DMSO. According to the experiment, the inclination of the curve for ( ) changes around ≈ 0.2. Moreover, the dielectric constant is almost linear for > 0.8. It is known that the TIP4P-2005 does not provide a reasonable description of the dielectric constant for pure water. Thus, if one changes the DMSO model, the results do not improve on the water-rich side of the composition, panel a of figure 10. On the DMSO-rich side, the OPLS and P2 models, if combined with TIP4P-2005 water model, yield not very bad values for . However, general trends of decay of with increasing are not satisfactory. In spite of long time of simulations, the values for strongly fluctuate along axis, such that it is impossible to capture any confiable peculiarities. Substantial improvement of the results is reached, if one substitutes TIP4P-2005 water model by the TIP4P/ model. Then, the values for dielectric constant essentially improve for water-rich mixtures. Moreover, the simulation data fluctuate less upon changing the composition . Two regions of behavior can be identified. Namely, the simulation data predict a fast decay of for water-rich mixtures and a slower decay in the opposite interval of compositions, say at > 0.3. The experimental points, by contrast, show a slower decay for a water-rich mixture and a slightly weaker decay at higher values of . After long runs (150 ns), we obtained with confidence a peculiarity in the values of in the interval of between 0.2 and 0.3. It is similar by shape to the behavior of mean square deviation of total dipole moment on composition plotted in figure 3 of reference [22] at between 0.1 and 0.15 with another force field. These trends,

33202-14
however, result from the modelling rather than from the experimental observations. In view of the present findings, the interpretation of data in terms of structural phase transition in laboratory systems seems to be questionable. Apparently, the inaccuracy of the description of energetic trends of mixing the water and DMSO species in the framework of P1 and P2 DMSO models leads to unsatisfactory dependence of the dielectric constant on . In qualitative terms, the conclusion is similar to what was mentioned above concerning the self-diffusion coefficients of species. However, it is worth to complement the calculations of the dielectric constant by the exploration of the behavior of the Kirkwood factor. There is room for intent improvement of the data. One can consider larger systems, exclude the fluctuations of pressure by switching to the canonical, NVT, ensemble, take larger values for cut-off distance. Application of these technical issues does not guarantee better results. Seemingly, present non-polarizable DMSO models with the present schemes of taking into account the long-range electrostatic interactions are incapable of reproducing the experimental data. Definitely, the situation improves if one applies even simple polarizable DMSO model, see e.g. [31]. Then, one would expect that the combination of the TIP4P/ water with even simplified polarizable DMSO model would be profitable. These issues guarantee future studies.

Surface tension of water-DMSO mixtures on composition
Our final remarks concern the behavior of the surface tension of water-DMSO mixtures. The surface tension calculations at each composition were performed by taking the final configuration of particles from the isobaric run. Next, the box edge along -axis was extended by a factor of 3, generating a rectangular box with a liquid slab and two liquid-mixture-vacuum interfaces in the − plane, in close similarity to the procedure applied in [65]. The total number of molecules is sufficient to yield sufficiently big area of the − face of the liquid slab. The elongation of the liquid slab along -axis is satisfactory as well. The executable file was modified by deleting a fixed pressure condition preserving the V-rescale thermostatting with the same parameters as in the NPT runs. Other corrections were not employed. The values for the surface tension, , follow from the combination of the time averages for the components of the pressure tensor, where are the components of the pressure tensor along , axes, and ... denotes the time average. We performed five runs at a constant volume, each piece of 10 ns, and obtained the result for by  taking the block average. The experimental results were taken from reference [66]. The data show that the surface tension rapidly decreases from the pure water value at = 0 till ≈ 0.25. Next, at higher values of , the values for decrease more slowly, the curve behaves almost linearly in that interval of compositions, panel a of figure 11. The excess surface tension from experiment exhibits minimum at ≈ 0.25. All the models in question behave similarly for water-rich mixtures. The decay of is better reproduced by P1-  figure 11. The shape of the surface tension on is quite well reproduced by these two models, just the absolute values of are underestimated. Thus, a reasonable description of and its excess simultaneously by the models of this study restricts to water-rich mixtures. Still, it seems that a better agreement with experimental data may result from the changes of fine details of simulation procedure as documented in reference [67].

Summary and conclusions
Recently we studied the properties of water-DMSO liquid mixtures in the entire composition interval by molecular dynamics computer simulations of the DMSO models combined with the SPC-E water model [33]. In the present work, we explore the same systems with the same methodology but using the non-polarizable, P1, P2 and OPLS united atom DMSO models in conjunction with the TIP4P-2005 and TIP4P/ water models. Concerning the DMSO models, we observed that the flexibility of angles, keeping the bond constrained, requires to include the dihedral angle as recommended in reference [42]. The systems in question are studied at room temperature, 298.15 K and atmospheric pressure, 0.1 MPa, as previously.
A set of novel findings was obtained and discussed in every detail. Namely, in the present work we obtained the excess partial molar volumes and apparent molar volumes. Analyses of these properties permitted us to capture an anomalous behavior observed at low DMSO molar fraction in agreement with experimental results. Besides, the excess partial molar enthalpy of each species was obtained and compared with experimental findings. Particular emphasis was put on the elucidation of the behavior of various radial distribution functions in water-rich mixtures. The coordination numbers and the fractions of molecules participating in hydrogen bonds were involved in the interpretation of the microscopic structure. The self-diffusion coefficients of species were obtained and compared with experimental data. Next, the static dielectric constant was calculated and compared with experimental results. Our final focus was on elucidating the behavior of the surface tension on the composition of water-DMSO mixtures and of the excess surface tension. Their trends were compared with experimental data. A very detailed validation of the predictions of the properties resulting from a set of models for water-DMSO solutions permit to make conclusions w.r.t. their applicability.
Among the variety of our results, solely the discussion of changes of the microscopic structure lacks comparison with experiments. In order to do that, one needs to have data concerning the total and desirably the partial structure factors in the entire composition interval. Then, it would be possible to perform the analyses similar to what was done in this laboratory for water-methanol mixtures [62]. We are not aware of such type of results for water-DMSO mixtures at present. This part of study of water-DMSO mixtures might benefit from the reverse Monte Carlo modelling in order to elucidate fine features of the hydrogen bonded network, see e.g., [56]. An illuminating finding about the existence and lifetime of frequently discussed complexes or associates may be expected.
According to our findings, one can distinguish two types of properties. Namely, for one type, the role of fluctuations is less important. Evidently, at a fixed pressure the systems' volume fluctuates, but density as an average over a set of computer simulation runs is well defined. The density dependence on composition compares well with experiment. Consequently, one is able to describe well the excess partial molar volumes and apparent molar volumes and definitely prove the existence of an anomalous behavior  for water-rich water-DMSO mixtures. Previously, this kind of behavior was obtained and discussed for water-methanol mixtures [68,69].
On the other hand, under the same conditions to obtain the isothermal compressibility, for example, requires time average of fluctuations. Consequently, this property is more difficult to obtain. Similarly, a disappointing situation is observed in the calculations of the static dielectric constant. It results from the time evolution of the total dipole moment fluctuations. If is calculated from the NPT runs, the volume fluctuations contribute to the inaccuracy of the results. In addition, at a constant density, the values for are affected by the choice of the number of particles in the simulation box, by the assumed cut-off of the interactions and by the scheme of summation of long-range interaction effects.
One should have these difficulties in mind while interpreting the changes of various properties on composition. Specifically, we observed that peculiarities of the behavior of different descriptors hardly permit to interpret them as discontinuous structural transitions in laboratory systems, see e.g., figures 3 and 5 of reference [22] and figures 8 and 9 of reference [20], unless a precise description of the underlying properties in terms of microscopic structure and/or the static dielectric constant is obtained by computer simulations.
To summarize, we would like to mention that each of the explored models was designed, or say parametrized, using a restricted number of different criteria at certain thermodynamic conditions. Thus, it is difficult to expect that modelling can be successful in a wide range of temperature and pressure. Similarly, we observed "non-universal" performance of the models along the composition axes. It appears that the P2-TIP4P-2005 model provides a reasonable description of a wide set of properties for water-rich mixtures, except the static dielectric constant. However, this behavior can be mitigated by using the P2 -TIP4P/ model. In addition, we would like to note that the P1-TIP4P-2005 model is better than the P2-TIP4P-2005 model for some properties of DMSO-rich mixtures and describes the surface tension and the excess surface tension reasonably well. Hence, one should choose the working model dependent on how concentrated DMSO solution is necessary to deal with, besides the pressure and temperature values of the given setup. The OPLS-TIP4P-2005 model is successful solely for a few properties. Hence, possible modifications of the parameters and combination rules should be attempted.
Finally, we would like to enlist a few missing elements to extend our knowledge of the properties of water-DMSO mixtures. It would be profitable to explore various velocity auto-correlation functions to enhance the understanding of dynamic properties, as well as to investigate Kirkwood factors and possibly the frequency dependent dielectric constant to understand the dielectric properties better. Another missing issue is to explore the cooperativity of hydrogen bond network in terms of certain complex elements, see e.g., [59], and hydrogen bonds life times.