Prediction of Transport Properties by Molecular Simulation: Methanol and Ethanol and their mixture
Gabriela Guevara-Carrion, Carlos Nieto-Draghi, Jadran Vrabec, Hans Hasse
aa r X i v : . [ phy s i c s . c h e m - ph ] J un Prediction of Transport Properties by Molecular Simulation:Methanol and Ethanol and their mixture
Gabriela Guevara-Carrion † , Carlos Nieto-Draghi ‡ , Jadran Vrabec § ∗ , and HansHasse †† Laboratory for Engineering Thermodynamics, University Kaiserslautern, 67663Kaiserslautern, Germany ‡ IFP, 1-4 Avenue de Bois Préau, 92852 Rueil-Malmaison Cedex, France § Institute of Thermodynamics and Thermal Process Engineering, University Stuttgart,70550 Stuttgart, Germany
Number of pages: 69Number of tables: 6Number of figures: 23 ∗ To whom correspondence should be addressed, tel.: +49-711/685-66107, fax: +49-711/685-66140, email: [email protected] bstract Transport properties of liquid methanol and ethanol are predicted by mole-cular dynamics simulation. The molecular models for the alcohols are rigid,non-polarizable and of united-atom type. They were developed in preced-ing work using experimental vapor-liquid equilibrium data only. Self- andMaxwell-Stefan diffusion coefficients as well as the shear viscosity of metha-nol, ethanol and their binary mixture are determined using equilibrium molec-ular dynamics and the Green-Kubo formalism. Non-equilibrium moleculardynamics is used for predicting the thermal conductivity of the two puresubstances. The transport properties of the fluids are calculated over a widetemperature range at ambient pressure and compared with experimental andsimulation data from the literature. Overall, a very good agreement with theexperiment is found. For instance, the self-diffusion coefficient and the shearviscosity are predicted with average deviations of less 8% for the pure alco-hols and 12% for the mixture. The predicted thermal conductivity agreesin average within 5% with the experimental data. Additionally, some ve-locity and shear viscosity autocorrelation functions are presented and dis-cussed. Radial distribution functions for ethanol are also presented. Thepredicted excess volume, excess enthalpy and the vapor-liquid equilibriumof the binary mixture methanol + ethanol are assessed and the vapor-liquidequilibrium agree well with experimental data.
Keywords:
Green-Kubo; reverse NEMD; self-diffusion coefficient; Maxwell-Stefan diffusion coefficient; shear viscosity; thermal conductivity.2
Introduction
There has been significant progress in recent years to study and predict both staticand dynamic thermophysical properties of fluids by molecular simulation. This isparticularly interesting for transport properties of liquids, where, due to the com-plexity of the involved physical mechanisms, other theoretical approaches are of-ten unsatisfactory. Hence, molecular dynamics is a promising method for predict-ing transport properties of liquids. However, the number of publications regardingtransport properties for technically relevant liquid systems that achieve quantita-tively accurate results is still low. The majority of these publications deals with theself-diffusion coefficient, leaving aside the more demanding Maxwell-Stefan dif-fusion coefficient, shear viscosity and thermal conductivity. Furthermore, molec-ular simulation studies on transport properties for complex liquids, e.g. mixturescontaining highly polar hydrogen-bonding molecules, are rare. This is also thecase for methanol, ethanol and their binary mixture. It should be pointed outthat there is a series of articles on transport properties of pure liquid ethanol byPetravi´c et al. [1, 2, 3]. There are several works on transport properties of associ-ating mixtures: Typically, aqueous solutions of methanol [4, 5, 6, 7, 8] or ethanol[7, 8, 9, 10, 11] have been studied using molecular dynamics simulation. How-ever, to our knowledge the mixture methanol + ethanol has not yet been regardedin that sense.Two fundamentally different methods for calculating transport properties bymolecular dynamics simulation are available. The equilibrium methods (EMD),using either the Green-Kubo formalism or the Einstein relations, analyze the time3ependent response of a fluid system to spontaneous fluctuations. With non-equilibrium molecular dynamics (NEMD), on the other hand, the system’s re-sponse to an externally applied perturbation is analyzed. NEMD was developedin order to increase the signal to noise ratio and to improve statistics and con-vergence. Both methods have different advantages and disadvantages, but arecomparable in efficiency, as shown e.g. by Dysthe et al. [12].The success of molecular dynamics simulation to predict thermodynamic prop-erties is highly determined by the potential model used to describe the intermolec-ular interactions. Numerous molecular models for methanol and ethanol havebeen proposed in the literature [9, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25].However, the ability of these models for predicting transport properties has beenprobed predominantly regarding the self-diffusion coefficient, cf. Table 1.For methanol, mainly the rigid OPLS molecular model from Jorgensen [13]and its modifications proposed by Haughney et al. [14] as well as by van Leeuwenand Smit [16], have been used to predict the self-diffusion coefficient. It wasshown by comparison of the models by Jorgensen [13] and by Haughney et al.[14] that the rigid OPLS model and its modifications perform well for the self-diffusion coefficient, deviations to experimental data range from 7 to 20% [14].Van de Ven-Lucassen et al. [5] showed that the model from van Leeuwen andSmit [16] predicts the self-diffusion coefficient of methanol at ambient tempera-ture more accurately than the polarizable model of Caldwell and Kollman [18],and than the models by Haughney et al. [14]. Wheeler et al. [26] predictedthe shear viscosity of methanol within 10% of experimental values. The thermal4onductivity of methanol at ambient conditions was reported by Nieto-Draghi [7]with a deviation of + − Molecular models
Throughout this work, rigid, non-polarizable molecular models of united-atomtype from earlier work of our group [32, 33] were used. Both models account forthe intermolecular interactions, including hydrogen bonding, by a set of united-atom Lennard-Jones (LJ) sites and superimposed point charges. Hence, the po-tential energy u i j between two alcohol molecules i and j can be written as u i j (cid:0) r i jab (cid:1) = n (cid:229) a = m (cid:229) b = e ab "(cid:18) s ab r i jab (cid:19) − (cid:18) s ab r i jab (cid:19) + q ia q jb pe r i jab , (1)where a is the site index of molecule i , b the site index of molecule j , and n and m are the numbers of interaction sites of molecules i and j , respectively. r i jab represents the site-site distance between molecules i and j . The LJ size and energyparameters are s ab and e ab . q ia and q jb are the point charges located at the sites a and b on the molecules i and j , and e is the permittivity of vacuum.The methanol model has two LJ sites, i.e. one for the methyl group and onefor the hydroxyl group. Additionally, three point charges are superimposed, i.e.one at each of the LJ site centers and one at the nucleus position of the hydroxylhydrogen. Ethanol was modeled by three LJ sites, i.e. one for the methyl, themethylene and the hydroxyl group each. Three point charges are located here aswell on the methylene and hydroxyl LJ site centers and on the nucleus position ofthe hydroxyl hydrogen.The parameter set of the methanol model was obtained using the optimiza-tion procedure of Stoll [34] on the basis of the geometry of the model from van7eeuwen and Smit [16]. The geometry of the ethanol model was derived fromab initio quantum chemistry calculations. The AUA4 parameters of Ungerer etal. [35] were applied for the methyl and methylene LJ sites. The point charge,LJ size and energy parameters as well as the offset of the hydroxyl LJ site werefitted to experimental vapor pressure and saturated liquid density. It should bepointed out that no information on transport or mixture properties was includedin the optimization procedures of any molecular model employed here so that thedata presented below is strictly predictive. The model parameters were taken from[32, 33] and are summarized in Table 2.To define a molecular model for a binary mixture on the basis of pairwise ad-ditive pure substance models, only the unlike interactions have to be specified. Incase of polar interaction sites, i.e. point charges here, this can straightforwardly bedone using the laws of electrostatics. However, for the unlike LJ parameters thereis no physically sound approach [36], and for predictions combination rules haveto be employed. Here, the interaction between unlike LJ sites of two moleculesare calculated applying the Lorentz-Berthelot combining rules: s ab = s aa + s bb , (2)and e ab = √ e aa e bb . (3)8 Methodology
Dynamic properties can be obtained from EMD simulations by means of theGreen-Kubo formalism [37, 38]. These equations are a direct relationship be-tween a transport coefficient and the time integral of an autocorrelation functionof a particular microscopic flux in a system in equilibrium. This method was usedin this work to calculate self- and Maxwell-Stefan diffusion coefficients as well asthe shear viscosity. The thermal conductivity was calculated using a modified ver-sion of the boundary driven (BD) reverse NEMD algorithm from Müller-Plathe[39], also referred to as PeX algorithm.
The self-diffusion coefficient D i is related to the mass flux of single moleculeswithin a fluid. Therefore, the relevant Green-Kubo expression is based on theindividual molecule velocity autocorrelation function as follows D i = N Z ¥ d t (cid:10) v k ( t ) · v k ( ) (cid:11) , (4)where v k ( t ) is the center of mass velocity vector of molecule k at some time t , and < ... > denotes the ensemble average. Eq. (4) is an average over all N moleculesin a simulation, since all contribute to the self-diffusion coefficient. In a binarymixture, the Maxwell-Stefan diffusion coefficient − D i j can also be written in termsof the center of mass velocity 9 D i j = x j N i (cid:18) x i M i + x j M j x j M j (cid:19) Z ¥ d t (cid:10) N i (cid:229) k = v ki ( t ) · N i (cid:229) k = v ki ( ) (cid:11) , (5)where M i and x i are the molar mass and mole fraction of species i . From theexpression above, the collective character of the Maxwell-Stefan diffusion coeffi-cient is evident.Since the present simulations provide self- as well as the Maxwell-Stefan dif-fusion coefficients simultaneously, a comparison to the simple predictive approachsuggested by Darken [40] is possible. The Maxwell-Stefan diffusion coefficient − D i j is estimated from the self-diffusion coefficients of the two components D i and D j in the binary mixture [40] − D i j = x i · D i + x j · D j · (6)This approach was found to be superior to the correlations by Vignes [41] andby Caldwell and Babb [42] in prior work on three binary mixtures [43]. How-ever, it is of little practical use due to the fact that it relies on the self-diffusioncoefficients in the mixture, which are rarely available. The shear viscosity h , as defined by Newton’s "law" of viscosity, is a measure ofthe resistance of a fluid to a shearing force [44]. It is associated with momentumtransport under the influence of velocity gradients. Hence, the shear viscosity can10e related to the time autocorrelation function of the off-diagonal elements of thestress tensor J xyp [45] h = V k B T Z ¥ d t (cid:10) J xyp ( t ) · J xyp ( ) (cid:11) , (7)where V stands for the molar volume, k B is the Boltzmann constant, and T denotesthe temperature. Averaging over all three independent elements of the stress ten-sor, i.e. J xyp , J xzp ,and J yzp , improves the statistics of the simulation. The component J xyp of the microscopic stress tensor J p is given in terms of the molecule positionsand velocities by [44] J xyp = N (cid:229) i = mv xi v yi − N (cid:229) i = N (cid:229) j = i n (cid:229) k = n (cid:229) l = r xi j ¶ u i j ¶ r ykl . (8)Here, the lower indices l and k count the interaction sites, and the upper indices x and y denote the spatial vector components, e.g. for velocity v xi or site-site dis-tance r xi j . The first term of Eq. (8) is the kinetic energy contribution and thesecond term is the potential energy contribution to the shear viscosity. Conse-quently, the Green-Kubo integral (7) can be decomposed into three parts, i.e.the kinetic energy contribution, the potential energy contribution and the mixedkinetic-potential energy contribution [44]. Eqs. (7) and (8) may directly be ap-plied to mixtures. 11 .3 Thermal conductivity The thermal conductivity l , as defined by Fourier’s "law" of heat conduction,characterizes the capability of a substance for molecular energy transport drivenby temperature gradients. In BD-NEMD simulation used in this work, two slabs ofa given thickness, large enough to contain a sufficiently large number of moleculeson average but much smaller than the total length of the simulation volume in z direction, are defined. In order to create a heat flux in z direction, the moleculewith the highest kinetic energy in the cold slab is selected with a frequency u to perform a hypothetical elastic collision with the molecule having the lowestkinetic energy in the hot slab [46]. After such a virtual collision the two moleculeshave exchanged their momentum without changing their respective position inspace. The new velocity of the molecule in the cold region is then v newc = − v oldc + m c v oldc + m h v oldh m c + m h , (9)and, for the molecule in the hot region v newh = − v oldh + m c v oldc + m h v oldh m c + m h , (10)where m c and m h are the respective masses of the selected molecules in the coldand hot slabs. v oldh and v newh stand, respectively, for the velocities of the moleculesbefore and after the collision.During the elastic collision the molecules exchange their momentum. Thus,12he total momentum and the total energy of the system remain constant. The ki-netic energy exchanged by means of these elastic collisions determines the energyflux in the system. The heat flow density in the steady state is then given by h J c i t = At (cid:229) transfers m h (cid:16) ( v newh ) − ( v oldh ) (cid:17) . (11)In this expression, A is the cross sectional area in x and y direction, t is the timeinterval of measure in the steady state during which a given number of momentumtransfers have occurred. h J c i t is the heat flux in z direction, which can indistinctlybe evaluated in the cold or the hot slab. It is required that the number of moleculesin both slabs is sufficiently large so that the intermolecular interactions within theslabs will properly thermalize the velocity distribution before a new collision takesplace. The thermal gradient (cid:209) T z is obtained from a linear fit of the temperatureprofiles from the simulation. In the steady state, the thermal conductivity can beobtained by l = − h J c i t (cid:209) T z · (12) The self-diffusion coefficient of pure methanol and pure ethanol was predicted atambient pressure in the temperature range from 210 to 340 K. Present numerical13ata are given in Table 3. Figures 1 and 2 show the self-diffusion coefficient formethanol and ethanol as a function of temperature, respectively. The statisticalerror of the simulation data is estimated to be in the order of 1%.The self-diffusion coefficient of methanol shows a very good agreement withthe experimental values. At temperatures below 280 K, the self-diffusion coef-ficient was underpredicted by approximately − + − +
5% on average. Furthermore, the predicteddata reproduce the composition dependence well.The Maxwell-Stefan diffusion coefficient of the mixture methanol + ethanolcalculated by molecular simulation is reported in Table 4. In Figure 4, the pre-dicted values are plotted as a function of the composition. The statistical uncer-tainties are on average 15%. These large error bars are attributed to the collectivenature of the Maxwell-Stefan diffusion coefficient. Due to the lack of experi-mental data, present simulation results are only compared to the predictions ofDarken’s model [40], cf. Eq. (6). This model can be considered to be a goodapproximation here, in agreement to previous findings for three other mixtures[43].
The shear viscosity of pure methanol and pure ethanol was predicted at ambientpressure for a temperature range of about 100 K around ambient conditions and isreported in Table 3. The temperature dependence of the shear viscosity is shownin Figures 5 and 6 for methanol and ethanol, respectively. Statistical errors are in15he range from 10 to 17%.In the case of methanol, the predicted shear viscosity agrees very well withthe experimental data over the whole regarded temperature range with an averagedeviation of 8%. At higher temperatures, the predicted viscosities are closer to theexperimental ones than at lower temperatures, cf. Figure 5. It should be pointedout that the highest deviation, being 20%, was found at the lowest regarded tem-perature 220 K. Here, the high density in combination with strongly interactingmolecules causes large simulation uncertainties due to the long-time behavior ofthe shear viscosity autocorrelation function. Present shear viscosity predictionsare better than the EMD data of Wheeler et al. [26] and are comparable with theNEMD Edwald data of the same author.The calculated shear viscosity for ethanol shows a good agreement to exper-imental data with an average deviation of 8%. It can be noticed in Figure 6 thatthere is a tendency to underestimate the shear viscosity. Petravi´c and Delhom-melle [2] also underestimated the shear viscosity using the OPLS model fromJorgensen [13]. At high densities, however, present data better predict the shearviscosity. Zhao et al. [31] recently reported shear viscosity data using the OPLS-AA model [19] and NEMD which is in very good agreement with experimentaldata.The three contributions to the shear viscosity were also calculated. In agree-ment with the results of Meier et al. [48], the potential energy contribution isdominant for both studied alcohols, h pp > h kk < Present results for the thermal conductivity of pure methanol and pure ethanol atambient pressure are summarized in Table 3. The agreement with the experimen-tal data is remarkable with an average deviation of 5% for methanol and of 2%for ethanol, cf. Figures 8 and 9. The thermal conductivity for methanol is slightlyunderestimated by the present predictions, while for ethanol it is slightly overesti-mated. Previous results for the thermal conductivity of methanol [7] and ethanol[3], using the OPLS models, show deviations of about 11 and 10%, respectively.In contrast to present results, the calculated thermal conductivity by Petravi´c [3]lies below the experimental data. 17 .4 Time dependent behavior of the Green-Kubo integrands
Velocity autocorrelation function
The velocity autocorrelation function provides an interesting insight into the liq-uid state. Selected normalized velocity autocorrelation functions (VACF) of puremethanol and ethanol are shown in Figure 10. In agreement with findings in theliterature [49], VACF may have two minima at short times ( t < . T m = .
37 K) [51].18igure 11 shows the behavior of the integrals of the VACF given by Eq. (4).It can be seen that the self-diffusion coefficient converges in all regarded cases toits final value after less than 4 ps. Thus, after this time interval the long-time tailof the VACF does not contribute significantly anymore.The behavior of the VACF of methanol and ethanol was also considered in themixture. From Figures 12 and 13 it can be seen that the VACF of methanol andethanol in the mixture differ from those of the pure alcohols at the same tempera-ture and pressure. Methanol VACF exhibit a significantly increased backscatteringdue to the presence of ethanol, cf. Figure 12. This indicates that at short times,ethanol molecules further restrain the mobility of methanol molecules. EthanolVACF also show a composition dependent behavior. In this case, the backscat-tering is reduced as the concentration of methanol increases, cf. Figure 13. Thissuggests that methanol molecules prevent the formation of cage-like structureshere, allowing ethanol molecules to move more freely in the mixture.
Shear viscosity autocorrelation function
Alder et al. [52] were the first to suggest the existence of a long-time tail inthe shear viscosity autocorrelation function (SACF) near the phase boundary tosolidification due to the low compressibility of such liquid states. Luo and Ho-heisel [53, 54] showed that the long-time behavior of the SACF is determinedprimarily by long-lived orientational correlations, occurring in liquids contain-ing anisotropic molecules. This behavior causes convergence problems in theGreen-Kubo integral, since the SACF must be calculated over a relatively longtime interval. 19igure 14 shows selected normalized SACF of pure methanol. The presenceof a significant long-time tail can clearly be noticed in Figure 15 where the inte-grals of the SACF of methanol are plotted for three temperatures. As expected,the behavior of the long-time tail of the SACF is highly dependent on moleculespecies and state point. For pure methanol at a low temperature of 260 K, approx-imately 60% of the contribution to the total shear viscosity comes from the first2 ps. At a higher temperature of 340 K the first 2 ps contribute 90% to the totalshear viscosity. In the case of ethanol (not shown here), the SACF are significantfor even longer times, i.e. at 263 K the first 2 ps of the SACF contribute only 45%to the final value of the shear viscosity. At 338 K this contribution reaches 75%.This difference can be explained by the higher molecular anisotropy of ethanol.The kinetic, potential and mixed kinetic-potential energy contributions to theSACF were also considered in this work. Over the whole regarded range of statepoints, the behavior of the kinetic contribution is analogous to the monotonicchair-like decay of autocorrelation functions for non-attractive fluids. The cor-responding time integral exhibits no important contribution from the long-timebehavior, cf. Figure 16. The potential energy contribution shows at short timesoscillations similar to those of the VACF, cf. Figures 10 and 17. These oscillationshave been also related to the presence of "bound states" by Michels and Trappe-niers [50] and by Meier et al. [48]. The long-time tail is significant and approxi-mates the behavior of the total SACF. The mixed kinetic-potential energy contri-bution shows a completely different behavior at short times, cf. Figure 18, wheredata is presented for pure methanol at four temperatures. At very short times ( t < . ) , two roughly mirrored patterns can be seen, the SACF decay/raise20harply to reach a minimum/maximum and oscillate strongly during their decayto zero. Meier et al. [48], who investigated simple LJ fluids, only observed initialpositive peaks in the mixed kinetic-potential energy contribution. To our knowl-edge, the presence of negative initial peaks has not been reported before. Thelong-time tail behavior of the mixed contribution is difficult to analyze because ofthe presence of strong oscillations. Petravi´c et al. [1] have shown that internal degrees of freedom play a crucialrole on the relaxation of the hydrogen bonding dynamics of liquid ethanol. Forinstance, they found that "freezing" bending and torsion in the OPLS model ofethanol produces a slowdown of the dynamics of the hydrogen bonds at 273 K(they observed an increase of the hydrogen bond life time of about 27%). For themethanol model used in this work, an excellent agreement between simulated andexperimental site-site radial distribution functions and hydrogen bond statisticshas been reported in previous work [32]. In the present work, the hydrogen bonddynamics of the rigid ethanol model was compared with the results using a flexi-ble model version including bending and torsion. For this porpouse, a geometricalhydrogen bonding criterion with r OH below 2.6 Å, r OO below 3.5 Å and q ( OH ... O ) below 30 ◦ together with a life time probability accounting for intermittent hy-drogen bonds [55] was used. In the flexible version, bending and torsion angleswere allowed to evolve using a standard harmonic and a Fourier series potential,respectively. The potential constants were taken from the TraPPE intermolecu-lar potential for alcohols [22], keeping the equilibrium values of the bending and21orsion angles of the rigid model and all the remaining inter-molecular model pa-rameters as well as the bond distances between all sites. Simulations were thenperformed at a specified temperature of 273 K (using the Nosé-Hoover thermo-stat) and a density of 17 .
62 mol / l to compute the radial distribution functions, thehydrogen bond dynamics and self-diffusion coefficient for both, rigid and flexibleversion of the present ethanol model.The comparison of the radial distribution functions from the flexible and therigid model can be seen in Figure 19. It was found that there are only minor differ-ences in the liquid structure as revealed by the three radial distribution functions g OO , g OH and g HH . The position and height of the first peaks are almost identi-cal for both models, with the exception of g HH , where a more pronounced peakis observed for the rigid model, cf. Figure 19. Only small differences are notedaround the second peak of g OO and g OH , mainly related to the molecules hydrogenbonded in the second solvation shell. In the latter case, a more pronounced secondpeak was found for the rigid model. In general, it can be said that the inclusionof the flexibility does not strongly affect the structure of liquid ethanol around thefirst solvation shell.In the next step, differences in the hydrogen bond dynamics for the two mod-els were analyzed. Average life times of hydrogen bonds of 137 and 80 ps wereobtained for the rigid and flexible versions, respectively. The flexibility signifi-cantly reduces the hydrogen bond life time by a factor of 1 .
7. Petravi´c et al. foundin a similar investigation for the OPLS ethanol model a reduction in the hydrogenbond dynamics of about 1 .
3. Hence, the effect is more pronounced for the present22odel. The hydrogen bond dynamics has a significant impact on the transportproperties. Eg. at 273 K a self-diffusion coefficient of 0 . × − m / s wasobtained for the flexible version, which is in perfect agreement with the values re-ported by Petravi´c et al. [1]. This value is higher by a factor of 1 . . × − m / s) whichis in excellent agreement with the experimental data (0 . × − m / s) [56]. Thiscomparison shows the impact of the internal degrees of freedom on the hydrogenbond dynamics and the self-diffusion coefficient. Similar effects can be expectedalso for other transport properties. Flexible and rigid models of H-Bonding sys-tems obviously show significant differences regarding transport properties. Thisshould also be kept in mind when assessing the comparisons of different modelsshown throughout this work. Vapor-liquid equilibria (VLE) often serve as a benchmark to assess molecularmodels. This is a useful choice, due to the fact that VLE represent the dominantfeature in the region of fluid states. Indeed, they are nowadays preferably beingused in the development of molecular models by optimizing the model parametersto VLE properties. This approach was also followed in the development of bothregarded pure fluid models. In the next step, molecular models can be assessedon how they perform in mixtures. Based on experience for numerous binaries andternaries, very good success can be expected, if the pure substance models aregood. 23alculated VLE data of pure methanol and ethanol using present potentialmodels have been reported previously for a temperature range from 270 to 490 K[32, 33]. The average deviations found for methanol in vapor pressure, saturatedliquid density and heat of vaporization are of 1.1, 0.6, and 5.5%, respectively[32]. The corresponding reported average deviations for ethanol are of 3.7, 0.3and 0.9%, respectively [33].Present simulation results for the binary VLE are given in Table 5 and com-pared in Figures 20 and 21 to experimental data by Butcher and Robinson [57] andto the Peng-Robinson equation of state [58]. As in the molecular simulations noadjustment to binary data was carried out, also the Peng-Robinson equation wasused without any adjustment to such data ( k i j = 0). Figure 20 shows the almost"ideal" behavior of methanol + ethanol for two temperatures. In fact, these werethe highest two isotherms for which experimental data were available to us. At thehigher one (443.15 K), however, it can be seen from the experimental data that thebubble line is not fully linear. The comparison of simulation and experiment is ex-cellent throughout, both the pressure slope and the composition relation betweenvapor and liquid match good, only minor deviations are found in the ethanol-richcomposition range where slightly too high vapor compositions are obtained. Excess properties are basic mixing properties which contain real mixing behavior.In case of methanol + ethanol non-ideal mixing effects are very weak. At equimo-lar composition, where the mixing influence is strongest, experimental data sug-24est a volume increase of 0.017% only. Also for enthalpy, a tiny positive excesscontribution of 0.011 % was found by experiment. This resolution is difficult toachieve by molecular simulation.To determine excess properties, a straightforward approach was used here.Three simulations at specified temperature and pressure were performed, two forthe pure substances and one of the mixture at a given composition. A binaryexcess property given by y E = y − x · y − x · y , (13)where y E represents the excess volume v E or the excess enthalpy h E . Enthalpiesor volumes of the pure liquids methanol and ethanol as well as of their binarymixture are denoted by y , y and y , respectively.Table 6 contains present simulation data for ambient conditions which arecompared to experimental results [59, 60, 61, 62] in Figures 22 and 23. Theexcess volume, cf. Figure 22, agrees with the experiment in all cases within itsseemingly large error bars. However, it should be noted that the simulation uncer-tainty is only around 0.1% of the total volume of the mixture. A slight tendency ofthe simulation data towards higher excess volumes might be assumed. In Figure23, the present excess enthalpy is shown together with experimental data. Themaximal excess contribution from experiment is only 0.011%, which is well be-low the simulation uncertainty of around 0.1%. By inspection of the simulationdata set as a whole, an even better agreement with experiment can be assumed.25 Conclusion
Transport and mixture properties for methanol and ethanol were predicted bymolecular dynamics simulation on the basis of rigid, non-polarizable molecularmodels from preceding work. Self- and Maxwell-Stefan diffusion coefficients aswell as shear viscosity were calculated by EMD and the Green-Kubo formalism.Comparing to experimental data for the two pure alcohols, present self-diffusioncoefficient and shear viscosity data are very good to excellent over a temperaturerange of about 100 K. Slightly higher deviations were found for the shear vis-cosity at low temperatures, which are mainly due to the long-time behavior ofthe shear viscosity autocorrelation function. Regarding the mixture methanol +ethanol at ambient conditions over the full composition range, both self-diffusioncoefficients also excellently agree with experimental data. The same holds forshear viscosity, especially in the methanol-rich region. Maxwell-Stefan diffusioncoefficient data is, to our knowledge, not available in the literature. Present pre-dictions correspond very well with the model of Darken. The thermal conductivitywas determined by BD reverse NEMD. Over a temperature range of about 100 K,an excellent agreement with the experiment was found as well for both pure fluids.Furthermore, predicted VLE data for the mixture methanol + ethanol along twoisotherms correspond with experimental data practically throughout within thesimulation uncertainty. The same holds for excess volume and excess enthalpy atambient conditions in the full composition range. It has been demonstrated thatpresent rigid, non-polarizable models for ethanol and methanol show a better per-formance to reproduce transport properties when compared to standard rigid and26exible models available in the literature.27 ppendix Simulation details
In this work, EMD simulations for transport properties were made in two steps.In the first step, a short simulation in the isobaric-isothermal (
N pT ) ensemble wasperformed at ambient pressure to calculate the density at the desired temperature.In the second step, a canonic (
NV T ) ensemble simulation was performed at thistemperature and density, to determine the transport properties. The simulationswere carried out in a cubic box with periodic boundary conditions containing2048 molecules. Newton’s equations of motion were solved using a fifth-orderGear predictor-corrector numerical integrator. The temperature was controlledby velocity scaling. In all EMD simulations, the integration time step was 0 . r c =
15 Å. Electrostatic long-range correctionswere made using the reaction field technique with conducting boundary conditions ( e RF = ¥ ) . On the basis of a center of mass cut-off scheme, the LJ long-rangeinteractions were corrected using angle averaging [63]. The simulations wereequilibrated in the NV T ensemble over 10 time steps, followed by productionruns of 0 . × time steps. Diffusion coefficients and shear viscosity werecalculated using Eqs. (4) to (7) with up to 9 900 independent time origins of theautocorrelation functions. The sampling length of the autocorrelation functionsvaried between 9 and 14 ps. The separation between time origins was chosen sothat all autocorrelation functions have decayed at least to 1 / e of their normalizedvalue to guarantee their time independence [64]. The uncertainties of the predictedvalues were estimated using a block averaging method [65].In the NEMD simulations with the PeX algorithm for the thermal conductiv-28ty, 800 molecules were placed in a parallelepiped box with L z = L x = L y ,where L i is the length of the simulation box in the different spatial directions. Pe-riodic boundary conditions were applied in all directions. The system was thenequilibrated over 1 ns at the desired temperature and pressure by N pT simulationusing a weak coupling bath [66] with long-range corrections [65] for pressure andenergy. The resulting density of the equilibration run was then employed to gen-erate a new set of simulations in order to develop the thermal gradient in a run of2 ns using the NEMD scheme. In this case an exchange frequency u = . − (i.e. every 150 time steps) for the energy flux was used. Electrostatic interactionswere treated with the reaction field with conducting boundary conditions. Thethermal conductivity was then computed as the average of 4 different runs over 1ns. The integration of the equations of motion was performed with a time step of2 fs with a cut-off radius of 12 Å. A Verlet neighbor list was employed to improvethe performance of the simulations.The VLE calculations of the mixture methanol + ethanol were performed us-ing the Grand Equilibrium method [67] in combination with the gradual insertiontechnique [68]. The latter, a Monte Carlo based expanded ensemble method, wasapplied to determine the chemical potential in the liquid phase. In comparison toWidom’s test particle method [69], where a full molecule is inserted into the fluid,the gradual insertion method introduces a fluctuating molecule. This molecule un-dergoes changes in a set of discrete states of coupling with all other real moleculesof the fluid. Preferential sampling was performed in the vicinity of the fluctuatingmolecule to improve accuracy. The gradual insertion simulations were done with500 molecules in the liquid phase. Starting from a face-centered cubic lattice ar-29angement, the simulation runs were equilibrated over 2 000 cycles in the N pT ensemble. The production phase was performed over 30 000 Monte Carlo cycles.Further simulation parameters of these runs were taken from Vrabec et al. [70].To evaluate the excess properties, EMD simulations in the
N pT ensemblewere performed. The system contained 1372 molecules in a cubic box with pe-riodic boundary conditions. The molecular trajectories were sampled with thealgorithms and parameters as described above. The simulations were equilibratedin the
N pT ensemble over 80 000 time steps, followed by a production run over600 000 time steps. 30 eferences [1] Petravic, J.; Delhommele, J.
J. Chem. Phys. , 303.[2] Petravic, J.; Delhommele, J.
J. Chem. Phys. , 234509.[3] Petravic, J.
J. Chem. Phys. , 174503.[4] Wheeler, D. R.; Rowley, R. L.
Mol. Phys. , 555.[5] van de Ven-Lucassen, I. M. J. J.; Vlugt, T. J. H.; van der Zanden, A. J. J.;Kerkhof, P. J. A. M. Mol. Simul. , 79.[6] Hawlicka, E.; Swiatla-Wojcik, D. Phys. Chem. Chem. Phys. , 3175.[7] Nieto-Draghi, C. Ph.D. Thesis, Universitat Rovira i Virgili, Terragona, 2003.[8] Wensink, E. J. W.; Hoffmann, A.; van Mareen, P. J.; van der Spoel, D. J.Chem. Phys. , 37308.[9] Müller-Plathe, F.
Mol. Simul. , 133.[10] Zhang, L.; Wang, Q.; Liu, Y.-C.; Zhang, L.-Z. J. Chem. Phys. ,104502.[11] Noskov, S. Y.; Lamourex, G.; Roux, B.
J. Phys. Chem. B , 6705.[12] Dysthe, D. K.; Fuchs, A. H.; Rousseau, B.
J. Chem. Phys. , 4047.[13] Jorgensen, W. L.
J. Phys. Chem. , 1276.[14] Haughney, M.; Ferrario, M.; McDonald, I. R. J. Phys. Chem. , 4934.3115] Mayo, S. L.; Olafson, B. D.; Goddard III, W. A. J. Phys. Chem. ,8897.[16] van Leeuwen, M. E.; Smit, B. J. Phys. Chem. , 1831.[17] Gao, J.; Habibollazadeh, D.; Shao, L. J. Phys. Chem. , 16460.[18] Caldwell, J. W.; Kollman, P. A. J. Phys. Chem. , 6208.[19] Jorgensen, W. L.; Maxwell, D. S.; Tirado-Rives, J. J. Am. Chem. Soc. , 11225.[20] Potoff, J. J.; Errington, J. R.; Panagiotopoulos, A. Z.
Mol. Phys. ,1073.[21] Stubbs, J. M.; Chen, B.; Potoff, J. J.; Siepmann, J. I. Fluid. Phase Equil. , 301.[22] Chen, B.; Potoff, J. J.; Siepmann, J. I.
J. Phys. Chem. B , 3093.[23] Kahre, R.; Sum, A. K.; Nath, S. K.; de Pablo, J. J.
J. Phys. Chem. B , 10071.[24] Liem X. Dang, T. M. C.
J. Chem. Phys. , 9851.[25] Wang, S.; Cann, N.
J. Chem. Phys. , 214502.[26] Wheeler, D. R.; Fuller, N. G.; Rowley, R. L.
Mol. Phys. , 55.[27] Taylor, R. S.; Shields, R. L. J. Chem. Phys. , 12569.[28] Saiz, L.; Padró, J. A.; Guàrdia, E.
J. Phys. Chem. B , 78.3229] Patel, S.; Brooks III, C. L.
J. Chem. Phys. , 164502.[30] González, M. A.; Encisco, E.; Bermejo, F. J.; Bée, M.
J. Chem. Phys. , 8045.[31] Zhao, L.; Wang, X.; Wang, L.; Sun, H.
Fluid. Phase Equil. , 212.[32] Schnabel, T.; Srivastava, A.; Vrabec, J.; Hasse, H.
J. Phys. Chem. B , 9871.[33] Schnabel, T.; Vrabec, J.; Hasse, H.
Fluid Phase Equil. , 134.[34] Stoll, J.
Molecular Models for the Prediction of Thermophysical Propertiesof Pure Fluids and Mixtures ; Fortschritt-Berichte VDI, Reihe 3, 836; VDIVerlag: Düsseldorf, 2005.[35] Ungerer, P.; Beauais, C.; Delhommelle, J.; Boutin, A.; Fuchs, A. H.
J. Chem.Phys. , 5499.[36] Haslam, A. J.; Galindo A.; Jackson G.
Fluid Phase Equil. , 105.[37] Green, M. S.
J. Chem. Phys. , 398.[38] Kubo, R. J. Phys. Soc. Jpn. , 570.[39] Müller-Plathe, F. J. Chem. Phys. , 6082.[40] Darken, L.
AIME , 184.[41] Vignes, A.
Ind. Eng. Chem. Fund. , 189.[42] Caldwell, C. S.; Babb, A. L. Phys. Rev. A , 51.3343] Fernández, G. A.; Vrabec, J.; Hasse, H. Int. J. Thermophys. , 1389.[44] Hoheisel, C. Phys. Rep. , 111.[45] Gubbins, K. E.
Statistical Mechanics,
Vol. 1; The Chemical Society Burling-ton house: London, 1972.[46] Nieto-Draghi, C.; Avalos, J. B.
Mol. Phys. , 2303.[47] Pereira, J. C. G.; Catlow, C. R. A.; Price, G. D.
J. Phys. Chem. A ,1909.[48] Meier, K.; Laesecke, A.; Kabelac, S.
J. Chem. Phys. , 3671.[49] Alonso, J.; Bermejo, F. J.; García-Hernández, M.; Martínez, J. L.; Howells,W. S.
J. Mol. Struct. , 147.[50] Michels, P. J.; Trappeniers, N. J.
Chem. Phys. Lett. , 195.[51] Sindzingre, P.; Klein, M. L. J. Chem. Phys. , 4681.[52] Alder, B. J.; Gass, D. M.; Wainwright, T. J. Chem. Phys. , 3813.[53] Luo, H.; Hoheisel, C. J. Chem. Phys. , 8378.[54] Luo, H.; Hoheisel, C. Phys. Rev. A , 6421.[55] Nieto-Draghi, C.; Avalos, J. B; Rosseau, B. J. Chem. Phys. , 4782.[56] Karger, N.; Vardag, T.; Lüdemann, H. D.
J. Chem. Phys. , 3437.[57] Butcher, K. L.; Robinson, W. I. J. Appl. Chem. , 289.3458] Peng, D.-Y.; Robinson, D. B. Ind. Eng. Chem. Fund. , 59.[59] Benson, G. C.; Pflug, H. D. J. Chem. Eng. Data , 382.[60] Ogawa, H.; Murakami, S. J. Solution Chem. , 315.[61] Zarei, H. A.; Jalili, F.; Assadi, S. J. Chem. Eng. Data , 2517.[62] Pflug, H. D.; Pope, A. E.; Benson, G. C. J. Chem. Eng. Data , 408.[63] Lustig, R. Mol. Phys. , 175.[64] Schoen, M.; Hoheisel, C. Mol. Phys. , 33.[65] Allen, M. P.; Tidesley, D. J. Computer Simulation of Liquids,
J. Chem. Phys. , 3684.[67] Vrabec, J.; Hasse, H. Mol. Phys. , 3375.[68] Nezbeda, I.; Kolafa, J.
Mol. Simul. , 391.[69] Widom, B. J. Chem. Phys. , 2808.[70] Vrabec, J.; Kettler, M.; Hasse, H. Chem. Phys. Lett. , 431.[71] Casulleras, J.; Guardia, E.
Mol. Simul. , 273.[72] Asahi, N.; Nakamura, Y. J. Chem. Phys. , 9879.[73] Johnson, P. A.; Babb, A. L.
Chem. Rev. , 587.3574] Dullien, F. A. L. AIChE J. , 62.[75] Hurle, R. L.; Eastel, A. J.; Woolf, L. A. J. Chem. Soc. Faraday Trans. 1 , 769.[76] Meckl, S.; Zeidler, M. D. Mol. Phys. , 85.[77] Johnson, P. A.; Babb, A. L. J. Chem. Phys. , 14.[78] Rauf, M. A.; Steward, G. H.; Farhataziz. J. Chem. Eng. Data , 324.[79] Vargaftik, N. B.; Vinogradov, Y. K.; Yargin, V. S. Handbook of physicalProperties of liquids and Gases. Pure Substances and Mixtures,
J. Chem. Eng. Data , 46.[81] Poling, B. E.; Prausnitz, J. M.; O ′ Connell, J. P
The Properties of Gases andLiquids,
Thermal Conductivity. Nonmetallic Liq-uids and Gases;
Thermophysical Properties of Matter Vol.3; Plenum: NewYork, 1970.[83] Lide, D. R., Ed.
CRC Handbook of Chemistry and Physics, − − − k B and theelectronic charge is e .Site s aa e aa / k B q ia Å K e MethanolS
CH3 + OH − H − − + CH3 − S CH2 + OH − H − − + . T r D i h l K mol l − − m − s − − Pa s Wm − K − Methanol213 27.02 0.318 (5) 3.1 (4)220 26.80 0.433 (6) 2.0 (2) 0.22 (3)240 26.22 0.750 (8) 1.7 (2)260 25.63 1.21 (1) 1.0 (2) 0.207 (9)278.15 25.14 1.72 (1) 0.64 (8)288 24.81 2.09 (2) 0.62 (7)298.15 24.51 2.41 (2) 0.57 (6)300 24.50 0.191 (9)308.15 24.19 2.97 (2) 0.47 (5)318.15 23.91 3.54 (2) 0.41 (6)320 23.88 0.188 (8)328.15 23.58 4.15 (3) 0.42 (5)340.15 23.21 4.91 (3) 0.29 (5)Ethanol239 18.11 0.192 (6) 3.1 (3)253.15 17.99 0.250 (5) 2.4 (2) 0.18 (1)263 17.83 0.363 (7) 2.2 (3)273.15 17.62 0.493 (8) 1.7 (2) 0.178 (9)280 17.46 0.64 (1) 1.4 (2)298.15 17.08 1.07 (1) 1.1 (1) 0.171 (9)308.15 16.88 1.33 (2) 1.0 (2)320 16.66 1.68 (1) 0.8 (1)328.15 16.45 2.01 (2) 0.61 (8) 0.161 (8)333 16.31 2.20 (2) 0.60 (7)338.15 16.26 2.40 (2) 0.47 (7)39able 4: Density, self-diffusion and Maxwell-Stefan diffusion coefficient as wellas shear viscosity of the mixture methanol (1) + ethanol (2) at 298.15 K and 0 . x r D D − D h mol mol − mol l − − m − s − − m − s − − m − s − − Pa s0 27.02 1.07 (1)0.143 26.22 1.37 (2) 1.12 (1) 1.4 (3) 1.17 (9)0.305 25.63 1.52 (1) 1.25 (1) 1.8 (3) 0.76 (9)0.500 25.14 1.73 (1) 1.44 (1) 1.7 (2) 0.62 (8)0.570 24.81 1.82 (1) 1.53 (1) 1.6 (3) 0.64 (9)0.754 24.51 2.07 (1) 1.72 (2) 2.3 (3) 0.70 (9)0.796 23.91 2.14 (2) 1.82 (2) 2.0 (3) 0.59 (7)0.890 23.58 2.28 (1) 1.90 (2) 2.6 (3) 0.60 (8)1 23.21 4.95 (3) 40able 5: Vapor-liquid equilibrium data for the mixture methanol (1) + ethanol(2) at 393.15 K and 433.15 K from present molecular dynamics simulation. Thenumber in parenthesis indicates the statistical uncertainty in the last digit. x p y r ′ r ′′ D h v mol mol − MPa mol mol − mol l − mol l − kJ mol − T = 393.15 K0 0.551 (6) 0 14.84 (2) 0.200 (1) 33.24 (7)0.150 0.448 (9) 0.207 (5) 15.62 (2) 0.154 (3) 33.86 (6)0.282 0.482 (8) 0.369 (2) 16.35 (2) 0.168 (3) 33.32 (7)0.492 0.538 (9) 0.589 (2) 17.54 (3) 0.185 (3) 32.67 (8)0.678 0.588 (9) 0.752 (5) 18.87 (3) 0.212 (3) 31.92 (9)0.838 0.572 (9) 0.884 (3) 19.94 (3) 0.205 (3) 31.43 (8)1.0 0.767 (5) 1.0 21.30 (3) 0.294 (1) 29.74 (8) T = 433.15 K0 1.32 (3) 0 13.58 (4) 0.469 (2) 28.41 (7)0.150 1.29 (2) 0.192 (4) 14.34 (4) 0.453 (7) 28.52 (9)0.282 1.38 (2) 0.339 (5) 14.96 (4) 0.486 (7) 28.08 (9)0.488 1.48 (1) 0.549 (5) 16.09 (3) 0.529 (5) 27.50 (8)0.588 1.54 (2) 0.659 (5) 18.87 (3) 0.548 (7) 27.28 (8)0.674 1.58 (2) 0.734 (5) 17.22 (4) 0.560 (7) 27.07 (9)0.836 1.64 (2) 0.872 (2) 18.31 (4) 0.594 (7) 26.55 (8)0.942 1.72 (2) 0.957 (1) 19.15 (5) 0.633 (7) 26.12 (9)1.0 1.89 (3) 1.0 19.42 (4) 0.711 (2) 26.12 (9)41able 6: Volume, excess volume, enthalpy and excess enthalpy of the mixturemethanol (1) + ethanol (2) at 298 .
15 K and 0 . x v v E h h E mol mol − cm − mol − cm − mol − kJ mol − J mol − − .
46 (4)0.070 57.15 (5) 0.00 (7) − .
20 (4) 3 (52)0.200 54.94 (5) 0.01 (6) − .
75 (4) 2 (49)0.270 53.72 (5) 0.03 (6) − .
50 (4) − − .
22 (4) − − .
95 (4) 20 (45)0.500 49.65 (4) 0.02 (5) − .
69 (4) − − .
32 (4) 5 (44)0.697 46.19 (4) 0.04 (5) − .
96 (4) 20 (44)0.800 44.34 (4) 0.00 (5) − .
61 (4) 1 (45)0.905 42.48 (4) 0.00 (5) − .
24 (4) − − .
90 (3)42 ist of Figures . • ) are comparedto experimental data [56, 72, 73, 74, 75] ( + ) and to other simu-lation results [5] ( △ ), [14] ( ), [72] ( ◦ ), and [47] ( ). Presentsimulation error bars are within symbol size. . . . . . . . . . . . . 472 Temperature dependence of the self-diffusion coefficient of pureethanol at 0 . • ) are comparedto experimental data [56, 73, 74, 75, 76] ( + ) and to other simula-tion results [1] ( ), [25] ( ◦ ), [27] ( △ ), [28] (star), [29] ( ▽ ), [30](hexagon), and [47] ( ). Present simulation error bars are withinsymbol size. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 483 Composition dependence of the self-diffusion coefficients for themixture methanol + ethanol at 298 .
15 K and 0 . • ) and ethanol ( N ) are comparedto experimental data [77] ( ◦ ) and ( △ ). . . . . . . . . . . . . . . . 494 Composition dependence of the Maxwell-Stefan diffusion coef-ficient for the mixture methanol + ethanol at 298 .
15 K and 0 . • ) are compared to Darken’smodel [40] ( ◦ ), cf. Eq. (6) . . . . . . . . . . . . . . . . . . . . . 5043 Temperature dependence of the shear viscosity of pure methanolat 0 . • ) andEdwald sumation ( N ) are compared to experimental data [78, 79]( + ) as well as to other EMD ( ) and NEMD ( ◦ ) simulation results[26]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 516 Temperature dependence of the shear viscosity of pure ethanol at0 . • ) andEdwald summation ( N ) are compared to experimental data [78]( + ) and to other simulation results [2] ( ) and [31] ( ◦ ). . . . . . . 527 Composition dependence of the shear viscosity for the mixturemethanol + ethanol at 298 .
15 K and 0 . • ) are compared to experimental data [80] ( + ). . . . . . . 538 Temperature dependence of the thermal conductivity of pure methanolat 0 . • ) are compared to exper-imental data [79, 81, 82, 83] ( + ) and to other simulation results[7] ( ◦ ). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 549 Temperature dependence of the thermal conductivity of pure ethanolat 0 . • ) are compared to exper-imental data [79, 81, 82, 83] ( + ) and to other simulation results[3] ( ) and [7] ( ◦ ). . . . . . . . . . . . . . . . . . . . . . . . . . 5510 Velocity autocorrelation function at 0.1 MPa of pure methanol at180 K ( − ) and 340 K ( −− ) as well as of pure ethanol at 173 K( − · − ) and 333 .
15 K ( · · · ). . . . . . . . . . . . . . . . . . . . . . 56441 Integral of the velocity autocorrelation function at 0.1 MPa of puremethanol at 180 K ( − ) and 340 K ( −− ) as well as of ethanol at173 K ( − · − ) and 333 .
15 K ( · · · ). . . . . . . . . . . . . . . . . . 5712 Velocity autocorrelation function of methanol in methanol (1) +ethanol (2) at 298.15 K and 0.1 MPa for different compositions: x = . − ( · · · ), x = . − ( − · − ), x = . − ( −− ) and x = − ( − ). . . . . . . . . . . 5813 Velocity autocorrelation function of ethanol in methanol (1) +ethanol (2) at 298.15 K and 0.1 MPa for different compositions: x = . − ( · · · ), x = . − ( − · − ), x = . − ( −− ) and x = − ( − ). . . . . . . . . . . 5914 Shear viscosity autocorrelation function of pure methanol at 0.1MPa and 260 K ( − ), 278 .
15 K ( − · − ) and 340 K ( −− ). . . . . . . 6015 Integral of the shear viscosity autocorrelation function of puremethanol at 0.1 MPa and 260 K ( − ), 278 .
15 K ( − · − ) and 340K ( −− ). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6116 Kinetic energy contribution of the shear viscosity autocorrelationfunction and its integral (inset) for pure methanol at 308 K and0.1 MPa. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6217 Potential energy contribution to the shear viscosity autocorrela-tion function and its integral (inset) for pure methanol at 308 Kand 0.1 MPa. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 63458 Mixed kinetic-potential energy contribution to the shear viscosityautocorrelation function of pure methanol at 0.1 MPa and 260 K( − ), 278 K ( −− ), 298 K ( − · − ) and 308 K ( · · · ). . . . . . . . . . 6419 Comparison of the radial distribution functions obtained with therigid ( − ) and flexible ( −− ) model versions for ethanol at 273 Kand 17 .
62 mol/l. . . . . . . . . . . . . . . . . . . . . . . . . . . . 6520 Vapor-liquid equilibria of the mixture methanol + ethanol at 393 . .
15 K. Present simulation results ( • ) are compared to ex-perimental data [57] ( + ). The lines represent the Peng-Robinsonequation of state. . . . . . . . . . . . . . . . . . . . . . . . . . . 6621 Vapor to liquid mole fraction ratio of the mixture methanol +ethanol at 393 .
15 and 433 .
15 K. Present simulation results ( • ) arecompared to experimental data [57] ( + ). The lines represent thePeng-Robinson equation of state. . . . . . . . . . . . . . . . . . . 6722 Composition dependence of the excess volume for the mixturemethanol + ethanol at 298 .
15 K and 0.1 MPa. Present simulationresults ( • ) are compared to experimental data [59, 60, 61] ( + ). . . 6823 Composition dependence of the excess enthalpy for the mixturemethanol + ethanol at 298 .
15 K and 0.1 MPa. Present simulationresults ( • ) are compared to experimental data [62] ( + ). . . . . . . 6946igure 1:47igure 2:48igure 3:49igure 4:50igure 5:51igure 6:52igure 7:53igure 8:54igure 9:55igure 10:56igure 11:57igure 12:58igure 13:59igure 14:60igure 15:61igure 16:62igure 17:63igure 18:64igure 19:65igure 20:66igure 21:67igure 22:68igure 23:69 r X i v : . [ phy s i c s . c h e m - ph ] J un Prediction of Transport Properties by Molecular Simulation:Methanol and Ethanol and their mixture
Gabriela Guevara-Carrion † , Carlos Nieto-Draghi ‡ , Jadran Vrabec § ∗ , and HansHasse †† Laboratory for Engineering Thermodynamics, University Kaiserslautern, 67663Kaiserslautern, Germany ‡ IFP, 1-4 Avenue de Bois Préau, 92852 Rueil-Malmaison Cedex, France § Institute of Thermodynamics and Thermal Process Engineering, University Stuttgart,70550 Stuttgart, Germany
Number of pages: 69Number of tables: 6Number of figures: 23 ∗ To whom correspondence should be addressed, tel.: +49-711/685-66107, fax: +49-711/685-66140, email: [email protected] bstract Transport properties of liquid methanol and ethanol are predicted by mole-cular dynamics simulation. The molecular models for the alcohols are rigid,non-polarizable and of united-atom type. They were developed in preced-ing work using experimental vapor-liquid equilibrium data only. Self- andMaxwell-Stefan diffusion coefficients as well as the shear viscosity of metha-nol, ethanol and their binary mixture are determined using equilibrium molec-ular dynamics and the Green-Kubo formalism. Non-equilibrium moleculardynamics is used for predicting the thermal conductivity of the two puresubstances. The transport properties of the fluids are calculated over a widetemperature range at ambient pressure and compared with experimental andsimulation data from the literature. Overall, a very good agreement with theexperiment is found. For instance, the self-diffusion coefficient and the shearviscosity are predicted with average deviations of less 8% for the pure alco-hols and 12% for the mixture. The predicted thermal conductivity agreesin average within 5% with the experimental data. Additionally, some ve-locity and shear viscosity autocorrelation functions are presented and dis-cussed. Radial distribution functions for ethanol are also presented. Thepredicted excess volume, excess enthalpy and the vapor-liquid equilibriumof the binary mixture methanol + ethanol are assessed and the vapor-liquidequilibrium agree well with experimental data.
Keywords:
Green-Kubo; reverse NEMD; self-diffusion coefficient; Maxwell-Stefan diffusion coefficient; shear viscosity; thermal conductivity.2
Introduction
There has been significant progress in recent years to study and predict both staticand dynamic thermophysical properties of fluids by molecular simulation. This isparticularly interesting for transport properties of liquids, where, due to the com-plexity of the involved physical mechanisms, other theoretical approaches are of-ten unsatisfactory. Hence, molecular dynamics is a promising method for predict-ing transport properties of liquids. However, the number of publications regardingtransport properties for technically relevant liquid systems that achieve quantita-tively accurate results is still low. The majority of these publications deals with theself-diffusion coefficient, leaving aside the more demanding Maxwell-Stefan dif-fusion coefficient, shear viscosity and thermal conductivity. Furthermore, molec-ular simulation studies on transport properties for complex liquids, e.g. mixturescontaining highly polar hydrogen-bonding molecules, are rare. This is also thecase for methanol, ethanol and their binary mixture. It should be pointed outthat there is a series of articles on transport properties of pure liquid ethanol byPetravi´c et al. [ ? , ? , ? ]. There are several works on transport properties of associ-ating mixtures: Typically, aqueous solutions of methanol [ ? , ? , ? , ? , ? ] or ethanol[ ? , ? , ? , ? , ? ] have been studied using molecular dynamics simulation. However,to our knowledge the mixture methanol + ethanol has not yet been regarded in thatsense.Two fundamentally different methods for calculating transport properties bymolecular dynamics simulation are available. The equilibrium methods (EMD),using either the Green-Kubo formalism or the Einstein relations, analyze the time3ependent response of a fluid system to spontaneous fluctuations. With non-equilibrium molecular dynamics (NEMD), on the other hand, the system’s re-sponse to an externally applied perturbation is analyzed. NEMD was developedin order to increase the signal to noise ratio and to improve statistics and con-vergence. Both methods have different advantages and disadvantages, but arecomparable in efficiency, as shown e.g. by Dysthe et al. [ ? ].The success of molecular dynamics simulation to predict thermodynamic prop-erties is highly determined by the potential model used to describe the intermolec-ular interactions. Numerous molecular models for methanol and ethanol havebeen proposed in the literature [ ? , ? , ? , ? , ? , ? , ? , ? , ? , ? , ? , ? , ? , ? ]. However,the ability of these models for predicting transport properties has been probedpredominantly regarding the self-diffusion coefficient, cf. Table ?? .For methanol, mainly the rigid OPLS molecular model from Jorgensen [ ? ] andits modifications proposed by Haughney et al. [ ? ] as well as by van Leeuwen andSmit [ ? ], have been used to predict the self-diffusion coefficient. It was shownby comparison of the models by Jorgensen [ ? ] and by Haughney et al. [ ? ] thatthe rigid OPLS model and its modifications perform well for the self-diffusioncoefficient, deviations to experimental data range from 7 to 20% [ ? ]. Van deVen-Lucassen et al. [ ? ] showed that the model from van Leeuwen and Smit [ ? ]predicts the self-diffusion coefficient of methanol at ambient temperature moreaccurately than the polarizable model of Caldwell and Kollman [ ? ], and than themodels by Haughney et al. [ ? ]. Wheeler et al. [ ? ] predicted the shear viscosityof methanol within 10% of experimental values. The thermal conductivity of4ethanol at ambient conditions was reported by Nieto-Draghi [ ? ] with a deviationof + ? , ? ], semi-flexible [ ? ], flexible [ ? ] and polarizable [ ? , ? ]molecular models for ethanol have been used to predict the self-diffusion coef-ficient. Thereby, it has been found that rigid, united-atom molecular models pre-dict the self-diffusion coefficient more reliably than polarizable, flexible models.Most ethanol models assessed by Wang and Cann [ ? ] as well as the Jorgensenmodel [ ? ] investigated by Saiz et al. [ ? ] overestimate the self-diffusion coefficientrather significantly with deviations ranging from +8 to +50%. Allowing a po-tential model for bending results in a further increase in the predicted values [ ? ].González et al. [ ? ] found that rigid, non-polarizable models perform better thanthe polarizable ones for predictions of the self-diffusion coefficient. The shear vis-cosity of ethanol has been predicted using EMD [ ? ] and NEMD [ ? ]. Both worksachieved good results (average deviations are 11 and 7%, respectively) using theunited-atom model and the all-atom model of Jorgensen et al. [ ? , ? ], respectively.Simulations performed by Petravi´c [ ? ] underestimate the thermal conductivity ofliquid ethanol near ambient temperatures by approximately − ? , ? ]. Thesemodels have been optimized on the basis of experimental data of vapor pressureand saturated liquid density in the whole temperature range from the triple pointto the critical point. No experimental data for transport properties was thereby5aken into account. The most significant transport properties, i.e. self-diffusioncoefficient, shear viscosity and thermal conductivity, were regarded here for thepure liquids using EMD and NEMD. Furthermore, the two molecular models weretested regarding the prediction of self- and Maxwell-Stefan diffusion coefficientsand shear viscosity of the binary mixture. Additionally, binary vapor-liquid equi-libria as well as excess volume and excess enthalpy were predicted. The goal ofthis study is to demonstrate the ability of rigid, non-polarizable molecular models,adjusted to pure substance vapor-liquid equilibria only,to accurately describe thetransport properties of strongly polar and H-bonding liquids.This work is organized as follows: Firstly, the employed molecular modelsare briefly described. Then, the simulation techniques are explained, the resultsfor self-diffusion coefficient, shear viscosity and thermal conductivity of the pureliquids are presented and compared to experimental and simulation data from theliterature. The self- and Maxwell-Stefan diffusion coefficients as well as the shearviscosity for the mixture methanol + ethanol are presented and compared withexperimental data where possible. Furthermore, the velocity and shear viscosityautocorrelation functions are described and discussed. Simulated radial distribu-tion functions of ethanol are also presented. Thereafter, the simulation results forvapor-liquid equilibria, excess volume and excess enthalpy for the binary mixtureare given in comparison to experimental values. Finally, conclusions are drawn.Simulation details are summarized in the Appendix.6 Molecular models
Throughout this work, rigid, non-polarizable molecular models of united-atomtype from earlier work of our group [ ? , ? ] were used. Both models account for theintermolecular interactions, including hydrogen bonding, by a set of united-atomLennard-Jones (LJ) sites and superimposed point charges. Hence, the potentialenergy u i j between two alcohol molecules i and j can be written as u i j (cid:0) r i jab (cid:1) = n (cid:229) a = m (cid:229) b = e ab "(cid:18) s ab r i jab (cid:19) − (cid:18) s ab r i jab (cid:19) + q ia q jb pe r i jab , (1)where a is the site index of molecule i , b the site index of molecule j , and n and m are the numbers of interaction sites of molecules i and j , respectively. r i jab represents the site-site distance between molecules i and j . The LJ size and energyparameters are s ab and e ab . q ia and q jb are the point charges located at the sites a and b on the molecules i and j , and e is the permittivity of vacuum.The methanol model has two LJ sites, i.e. one for the methyl group and onefor the hydroxyl group. Additionally, three point charges are superimposed, i.e.one at each of the LJ site centers and one at the nucleus position of the hydroxylhydrogen. Ethanol was modeled by three LJ sites, i.e. one for the methyl, themethylene and the hydroxyl group each. Three point charges are located here aswell on the methylene and hydroxyl LJ site centers and on the nucleus position ofthe hydroxyl hydrogen.The parameter set of the methanol model was obtained using the optimizationprocedure of Stoll [ ? ] on the basis of the geometry of the model from van Leeuwen7nd Smit [ ? ]. The geometry of the ethanol model was derived from ab initioquantum chemistry calculations. The AUA4 parameters of Ungerer et al. [ ? ]were applied for the methyl and methylene LJ sites. The point charge, LJ sizeand energy parameters as well as the offset of the hydroxyl LJ site were fittedto experimental vapor pressure and saturated liquid density. It should be pointedout that no information on transport or mixture properties was included in theoptimization procedures of any molecular model employed here so that the datapresented below is strictly predictive. The model parameters were taken from[ ? , ? ] and are summarized in Table ?? .To define a molecular model for a binary mixture on the basis of pairwise ad-ditive pure substance models, only the unlike interactions have to be specified. Incase of polar interaction sites, i.e. point charges here, this can straightforwardly bedone using the laws of electrostatics. However, for the unlike LJ parameters thereis no physically sound approach [ ? ], and for predictions combination rules haveto be employed. Here, the interaction between unlike LJ sites of two moleculesare calculated applying the Lorentz-Berthelot combining rules: s ab = s aa + s bb , (2)and e ab = √ e aa e bb . (3)8 Methodology
Dynamic properties can be obtained from EMD simulations by means of theGreen-Kubo formalism [ ? , ? ]. These equations are a direct relationship between atransport coefficient and the time integral of an autocorrelation function of a par-ticular microscopic flux in a system in equilibrium. This method was used in thiswork to calculate self- and Maxwell-Stefan diffusion coefficients as well as theshear viscosity. The thermal conductivity was calculated using a modified versionof the boundary driven (BD) reverse NEMD algorithm from Müller-Plathe [ ? ],also referred to as PeX algorithm. The self-diffusion coefficient D i is related to the mass flux of single moleculeswithin a fluid. Therefore, the relevant Green-Kubo expression is based on theindividual molecule velocity autocorrelation function as follows D i = N Z ¥ d t (cid:10) v k ( t ) · v k ( ) (cid:11) , (4)where v k ( t ) is the center of mass velocity vector of molecule k at some time t , and < ... > denotes the ensemble average. Eq. ( ?? ) is an average over all N moleculesin a simulation, since all contribute to the self-diffusion coefficient. In a binarymixture, the Maxwell-Stefan diffusion coefficient − D i j can also be written in termsof the center of mass velocity 9 D i j = x j N i (cid:18) x i M i + x j M j x j M j (cid:19) Z ¥ d t (cid:10) N i (cid:229) k = v ki ( t ) · N i (cid:229) k = v ki ( ) (cid:11) , (5)where M i and x i are the molar mass and mole fraction of species i . From theexpression above, the collective character of the Maxwell-Stefan diffusion coeffi-cient is evident.Since the present simulations provide self- as well as the Maxwell-Stefan dif-fusion coefficients simultaneously, a comparison to the simple predictive approachsuggested by Darken [ ? ] is possible. The Maxwell-Stefan diffusion coefficient − D i j is estimated from the self-diffusion coefficients of the two components D i and D j in the binary mixture [ ? ] − D i j = x i · D i + x j · D j · (6)This approach was found to be superior to the correlations by Vignes [ ? ] andby Caldwell and Babb [ ? ] in prior work on three binary mixtures [ ? ]. However, itis of little practical use due to the fact that it relies on the self-diffusion coefficientsin the mixture, which are rarely available. The shear viscosity h , as defined by Newton’s "law" of viscosity, is a measure ofthe resistance of a fluid to a shearing force [ ? ]. It is associated with momentumtransport under the influence of velocity gradients. Hence, the shear viscosity can10e related to the time autocorrelation function of the off-diagonal elements of thestress tensor J xyp [ ? ] h = V k B T Z ¥ d t (cid:10) J xyp ( t ) · J xyp ( ) (cid:11) , (7)where V stands for the molar volume, k B is the Boltzmann constant, and T denotesthe temperature. Averaging over all three independent elements of the stress ten-sor, i.e. J xyp , J xzp ,and J yzp , improves the statistics of the simulation. The component J xyp of the microscopic stress tensor J p is given in terms of the molecule positionsand velocities by [ ? ] J xyp = N (cid:229) i = mv xi v yi − N (cid:229) i = N (cid:229) j = i n (cid:229) k = n (cid:229) l = r xi j ¶ u i j ¶ r ykl . (8)Here, the lower indices l and k count the interaction sites, and the upper indices x and y denote the spatial vector components, e.g. for velocity v xi or site-site dis-tance r xi j . The first term of Eq. ( ?? ) is the kinetic energy contribution and thesecond term is the potential energy contribution to the shear viscosity. Conse-quently, the Green-Kubo integral ( ?? ) can be decomposed into three parts, i.e.the kinetic energy contribution, the potential energy contribution and the mixedkinetic-potential energy contribution [ ? ]. Eqs. ( ?? ) and ( ?? ) may directly be ap-plied to mixtures. 11 .3 Thermal conductivity The thermal conductivity l , as defined by Fourier’s "law" of heat conduction,characterizes the capability of a substance for molecular energy transport drivenby temperature gradients. In BD-NEMD simulation used in this work, two slabs ofa given thickness, large enough to contain a sufficiently large number of moleculeson average but much smaller than the total length of the simulation volume in z direction, are defined. In order to create a heat flux in z direction, the moleculewith the highest kinetic energy in the cold slab is selected with a frequency u to perform a hypothetical elastic collision with the molecule having the lowestkinetic energy in the hot slab [ ? ]. After such a virtual collision the two moleculeshave exchanged their momentum without changing their respective position inspace. The new velocity of the molecule in the cold region is then v newc = − v oldc + m c v oldc + m h v oldh m c + m h , (9)and, for the molecule in the hot region v newh = − v oldh + m c v oldc + m h v oldh m c + m h , (10)where m c and m h are the respective masses of the selected molecules in the coldand hot slabs. v oldh and v newh stand, respectively, for the velocities of the moleculesbefore and after the collision.During the elastic collision the molecules exchange their momentum. Thus,12he total momentum and the total energy of the system remain constant. The ki-netic energy exchanged by means of these elastic collisions determines the energyflux in the system. The heat flow density in the steady state is then given by h J c i t = At (cid:229) transfers m h (cid:16) ( v newh ) − ( v oldh ) (cid:17) . (11)In this expression, A is the cross sectional area in x and y direction, t is the timeinterval of measure in the steady state during which a given number of momentumtransfers have occurred. h J c i t is the heat flux in z direction, which can indistinctlybe evaluated in the cold or the hot slab. It is required that the number of moleculesin both slabs is sufficiently large so that the intermolecular interactions within theslabs will properly thermalize the velocity distribution before a new collision takesplace. The thermal gradient (cid:209) T z is obtained from a linear fit of the temperatureprofiles from the simulation. In the steady state, the thermal conductivity can beobtained by l = − h J c i t (cid:209) T z · (12) The self-diffusion coefficient of pure methanol and pure ethanol was predicted atambient pressure in the temperature range from 210 to 340 K. Present numerical13ata are given in Table ?? . Figures ?? and ?? show the self-diffusion coefficientfor methanol and ethanol as a function of temperature, respectively. The statisticalerror of the simulation data is estimated to be in the order of 1%.The self-diffusion coefficient of methanol shows a very good agreement withthe experimental values. At temperatures below 280 K, the self-diffusion coef-ficient was underpredicted by approximately − + ? ], using amodified Jorgensen model, is comparable with present results near ambient tem-perature where average deviations are 9 %. However, the data from Haughney etal. is scattering and the self-diffusion coefficient is significantly underestimatedabove 300 K.There is an excellent agreement between experimental data and present self-diffusion coefficient for ethanol. The predicted data also correctly reproducesthe temperature dependence over the whole temperature range from 239 to 338 K.Contrary to previous simulation reports [ ? , ? , ? , ? , ? , ? , ? ], where the self-diffusioncoefficient of ethanol has mostly been significantly overestimated, deviations ofthe present data are negative and on average only − ?? and plotted in Figure ?? . The estimatedstatistical uncertainty of the simulation data is also 1%. As shown in Figure ?? ,there is an excellent agreement between predicted self-diffusion coefficients andthe experimental data for both alcohols. Present simulation results slightly over-estimate the self-diffusion of both alcohols by +
5% on average. Furthermore, thepredicted data reproduce the composition dependence well.The Maxwell-Stefan diffusion coefficient of the mixture methanol + ethanolcalculated by molecular simulation is reported in Table ?? . In Figure ?? , the pre-dicted values are plotted as a function of the composition. The statistical uncer-tainties are on average 15%. These large error bars are attributed to the collectivenature of the Maxwell-Stefan diffusion coefficient. Due to the lack of experi-mental data, present simulation results are only compared to the predictions ofDarken’s model [ ? ], cf. Eq. ( ?? ). This model can be considered to be a goodapproximation here, in agreement to previous findings for three other mixtures[ ? ]. The shear viscosity of pure methanol and pure ethanol was predicted at ambientpressure for a temperature range of about 100 K around ambient conditions and isreported in Table ?? . The temperature dependence of the shear viscosity is shownin Figures ?? and ?? for methanol and ethanol, respectively. Statistical errors are15n the range from 10 to 17%.In the case of methanol, the predicted shear viscosity agrees very well withthe experimental data over the whole regarded temperature range with an averagedeviation of 8%. At higher temperatures, the predicted viscosities are closer to theexperimental ones than at lower temperatures, cf. Figure ?? . It should be pointedout that the highest deviation, being 20%, was found at the lowest regarded tem-perature 220 K. Here, the high density in combination with strongly interactingmolecules causes large simulation uncertainties due to the long-time behavior ofthe shear viscosity autocorrelation function. Present shear viscosity predictionsare better than the EMD data of Wheeler et al. [ ? ] and are comparable with theNEMD Edwald data of the same author.The calculated shear viscosity for ethanol shows a good agreement to exper-imental data with an average deviation of 8%. It can be noticed in Figure ?? that there is a tendency to underestimate the shear viscosity. Petravi´c and Del-hommelle [ ? ] also underestimated the shear viscosity using the OPLS model fromJorgensen [ ? ]. At high densities, however, present data better predict the shearviscosity. Zhao et al. [ ? ] recently reported shear viscosity data using the OPLS-AA model [ ? ] and NEMD which is in very good agreement with experimentaldata.The three contributions to the shear viscosity were also calculated. In agree-ment with the results of Meier et al. [ ? ], the potential energy contribution is domi-nant for both studied alcohols, h pp > h kk < ?? and plotted in Figure ?? together with the available experimental data.There is a very good agreement between experimental and predicted shear vis-cosity, being mostly within the simulation uncertainty of around 12%. Deviationsfrom the experimental data are on average 12% as well. Present results for the thermal conductivity of pure methanol and pure ethanol atambient pressure are summarized in Table ?? . The agreement with the experi-mental data is remarkable with an average deviation of 5% for methanol and of2% for ethanol, cf. Figures ?? and ?? . The thermal conductivity for methanol isslightly underestimated by the present predictions, while for ethanol it is slightlyoverestimated. Previous results for the thermal conductivity of methanol [ ? ] andethanol [ ? ], using the OPLS models, show deviations of about 11 and 10%, re-spectively. In contrast to present results, the calculated thermal conductivity byPetravi´c [ ? ] lies below the experimental data.17 .4 Time dependent behavior of the Green-Kubo integrands Velocity autocorrelation function
The velocity autocorrelation function provides an interesting insight into the liq-uid state. Selected normalized velocity autocorrelation functions (VACF) of puremethanol and ethanol are shown in Figure ?? . In agreement with findings in theliterature [ ? ], VACF may have two minima at short times ( t < . ? ] attributed this phe-nomenon to the formation of “bound states” since this effect can only be observedfor molecular models that include attractive forces.For higher temperatures at ambient pressure, where the liquid density is lower,molecular collisions are less frequent and the depth of the minima is graduallyreduced until they almost disappear, cf. Figure ?? . At short times, the VACF decayfaster at lower temperatures. Moreover, the initial decay of the VACF of methanoloccurs earlier than in the corresponding VACF of ethanol. Methanol VACF exhibita more pronounced backscattering oscillations than those of ethanol, particularlyat low temperatures near its melting temperature ( T m = .
37 K) [ ? ].18igure ?? shows the behavior of the integrals of the VACF given by Eq. ( ?? ).It can be seen that the self-diffusion coefficient converges in all regarded cases toits final value after less than 4 ps. Thus, after this time interval the long-time tailof the VACF does not contribute significantly anymore.The behavior of the VACF of methanol and ethanol was also considered in themixture. From Figures ?? and ?? it can be seen that the VACF of methanol andethanol in the mixture differ from those of the pure alcohols at the same tempera-ture and pressure. Methanol VACF exhibit a significantly increased backscatteringdue to the presence of ethanol, cf. Figure ?? . This indicates that at short times,ethanol molecules further restrain the mobility of methanol molecules. EthanolVACF also show a composition dependent behavior. In this case, the backscat-tering is reduced as the concentration of methanol increases, cf. Figure ?? . Thissuggests that methanol molecules prevent the formation of cage-like structureshere, allowing ethanol molecules to move more freely in the mixture. Shear viscosity autocorrelation function
Alder et al. [ ? ] were the first to suggest the existence of a long-time tail in theshear viscosity autocorrelation function (SACF) near the phase boundary to solid-ification due to the low compressibility of such liquid states. Luo and Hoheisel[ ? , ? ] showed that the long-time behavior of the SACF is determined primarily bylong-lived orientational correlations, occurring in liquids containing anisotropicmolecules. This behavior causes convergence problems in the Green-Kubo inte-gral, since the SACF must be calculated over a relatively long time interval.19igure ?? shows selected normalized SACF of pure methanol. The presenceof a significant long-time tail can clearly be noticed in Figure ?? where the inte-grals of the SACF of methanol are plotted for three temperatures. As expected,the behavior of the long-time tail of the SACF is highly dependent on moleculespecies and state point. For pure methanol at a low temperature of 260 K, approx-imately 60% of the contribution to the total shear viscosity comes from the first2 ps. At a higher temperature of 340 K the first 2 ps contribute 90% to the totalshear viscosity. In the case of ethanol (not shown here), the SACF are significantfor even longer times, i.e. at 263 K the first 2 ps of the SACF contribute only 45%to the final value of the shear viscosity. At 338 K this contribution reaches 75%.This difference can be explained by the higher molecular anisotropy of ethanol.The kinetic, potential and mixed kinetic-potential energy contributions to theSACF were also considered in this work. Over the whole regarded range of statepoints, the behavior of the kinetic contribution is analogous to the monotonicchair-like decay of autocorrelation functions for non-attractive fluids. The cor-responding time integral exhibits no important contribution from the long-timebehavior, cf. Figure ?? . The potential energy contribution shows at short timesoscillations similar to those of the VACF, cf. Figures ?? and ?? . These oscillationshave been also related to the presence of "bound states" by Michels and Trappe-niers [ ? ] and by Meier et al. [ ? ]. The long-time tail is significant and approximatesthe behavior of the total SACF. The mixed kinetic-potential energy contributionshows a completely different behavior at short times, cf. Figure ?? , where data ispresented for pure methanol at four temperatures. At very short times ( t < . ) ,two roughly mirrored patterns can be seen, the SACF decay/raise sharply to reach20 minimum/maximum and oscillate strongly during their decay to zero. Meier etal. [ ? ], who investigated simple LJ fluids, only observed initial positive peaks inthe mixed kinetic-potential energy contribution. To our knowledge, the presenceof negative initial peaks has not been reported before. The long-time tail behaviorof the mixed contribution is difficult to analyze because of the presence of strongoscillations. Petravi´c et al. [ ? ] have shown that internal degrees of freedom play a crucialrole on the relaxation of the hydrogen bonding dynamics of liquid ethanol. Forinstance, they found that "freezing" bending and torsion in the OPLS model ofethanol produces a slowdown of the dynamics of the hydrogen bonds at 273 K(they observed an increase of the hydrogen bond life time of about 27%). For themethanol model used in this work, an excellent agreement between simulated andexperimental site-site radial distribution functions and hydrogen bond statisticshas been reported in previous work [ ? ]. In the present work, the hydrogen bonddynamics of the rigid ethanol model was compared with the results using a flexi-ble model version including bending and torsion. For this porpouse, a geometricalhydrogen bonding criterion with r OH below 2.6 Å, r OO below 3.5 Å and q ( OH ... O ) below 30 ◦ together with a life time probability accounting for intermittent hy-drogen bonds [ ? ] was used. In the flexible version, bending and torsion angleswere allowed to evolve using a standard harmonic and a Fourier series potential,respectively. The potential constants were taken from the TraPPE intermolecu-lar potential for alcohols [ ? ], keeping the equilibrium values of the bending and21orsion angles of the rigid model and all the remaining inter-molecular model pa-rameters as well as the bond distances between all sites. Simulations were thenperformed at a specified temperature of 273 K (using the Nosé-Hoover thermo-stat) and a density of 17 .
62 mol / l to compute the radial distribution functions, thehydrogen bond dynamics and self-diffusion coefficient for both, rigid and flexibleversion of the present ethanol model.The comparison of the radial distribution functions from the flexible and therigid model can be seen in Figure ?? . It was found that there are only minor differ-ences in the liquid structure as revealed by the three radial distribution functions g OO , g OH and g HH . The position and height of the first peaks are almost identi-cal for both models, with the exception of g HH , where a more pronounced peakis observed for the rigid model, cf. Figure ?? . Only small differences are notedaround the second peak of g OO and g OH , mainly related to the molecules hydrogenbonded in the second solvation shell. In the latter case, a more pronounced secondpeak was found for the rigid model. In general, it can be said that the inclusionof the flexibility does not strongly affect the structure of liquid ethanol around thefirst solvation shell.In the next step, differences in the hydrogen bond dynamics for the two mod-els were analyzed. Average life times of hydrogen bonds of 137 and 80 ps wereobtained for the rigid and flexible versions, respectively. The flexibility signifi-cantly reduces the hydrogen bond life time by a factor of 1 .
7. Petravi´c et al. foundin a similar investigation for the OPLS ethanol model a reduction in the hydrogenbond dynamics of about 1 .
3. Hence, the effect is more pronounced for the present22odel. The hydrogen bond dynamics has a significant impact on the transportproperties. Eg. at 273 K a self-diffusion coefficient of 0 . × − m / s wasobtained for the flexible version, which is in perfect agreement with the values re-ported by Petravi´c et al. [ ? ]. This value is higher by a factor of 1 . . × − m / s) whichis in excellent agreement with the experimental data (0 . × − m / s) [ ? ]. Thiscomparison shows the impact of the internal degrees of freedom on the hydrogenbond dynamics and the self-diffusion coefficient. Similar effects can be expectedalso for other transport properties. Flexible and rigid models of H-Bonding sys-tems obviously show significant differences regarding transport properties. Thisshould also be kept in mind when assessing the comparisons of different modelsshown throughout this work. Vapor-liquid equilibria (VLE) often serve as a benchmark to assess molecularmodels. This is a useful choice, due to the fact that VLE represent the dominantfeature in the region of fluid states. Indeed, they are nowadays preferably beingused in the development of molecular models by optimizing the model parametersto VLE properties. This approach was also followed in the development of bothregarded pure fluid models. In the next step, molecular models can be assessedon how they perform in mixtures. Based on experience for numerous binaries andternaries, very good success can be expected, if the pure substance models aregood. 23alculated VLE data of pure methanol and ethanol using present potentialmodels have been reported previously for a temperature range from 270 to 490K [ ? , ? ]. The average deviations found for methanol in vapor pressure, saturatedliquid density and heat of vaporization are of 1.1, 0.6, and 5.5%, respectively[ ? ]. The corresponding reported average deviations for ethanol are of 3.7, 0.3 and0.9%, respectively [ ? ].Present simulation results for the binary VLE are given in Table ?? and com-pared in Figures ?? and ?? to experimental data by Butcher and Robinson [ ? ] andto the Peng-Robinson equation of state [ ? ]. As in the molecular simulations noadjustment to binary data was carried out, also the Peng-Robinson equation wasused without any adjustment to such data ( k i j = 0). Figure ?? shows the almost"ideal" behavior of methanol + ethanol for two temperatures. In fact, these werethe highest two isotherms for which experimental data were available to us. At thehigher one (443.15 K), however, it can be seen from the experimental data that thebubble line is not fully linear. The comparison of simulation and experiment is ex-cellent throughout, both the pressure slope and the composition relation betweenvapor and liquid match good, only minor deviations are found in the ethanol-richcomposition range where slightly too high vapor compositions are obtained. Excess properties are basic mixing properties which contain real mixing behavior.In case of methanol + ethanol non-ideal mixing effects are very weak. At equimo-lar composition, where the mixing influence is strongest, experimental data sug-24est a volume increase of 0.017% only. Also for enthalpy, a tiny positive excesscontribution of 0.011 % was found by experiment. This resolution is difficult toachieve by molecular simulation.To determine excess properties, a straightforward approach was used here.Three simulations at specified temperature and pressure were performed, two forthe pure substances and one of the mixture at a given composition. A binaryexcess property given by y E = y − x · y − x · y , (13)where y E represents the excess volume v E or the excess enthalpy h E . Enthalpiesor volumes of the pure liquids methanol and ethanol as well as of their binarymixture are denoted by y , y and y , respectively.Table ?? contains present simulation data for ambient conditions which arecompared to experimental results [ ? , ? , ? , ? ] in Figures ?? and ?? . The excessvolume, cf. Figure ?? , agrees with the experiment in all cases within its seem-ingly large error bars. However, it should be noted that the simulation uncertaintyis only around 0.1% of the total volume of the mixture. A slight tendency of thesimulation data towards higher excess volumes might be assumed. In Figure ?? ,the present excess enthalpy is shown together with experimental data. The maxi-mal excess contribution from experiment is only 0.011%, which is well below thesimulation uncertainty of around 0.1%. By inspection of the simulation data setas a whole, an even better agreement with experiment can be assumed.25 Conclusion
Transport and mixture properties for methanol and ethanol were predicted bymolecular dynamics simulation on the basis of rigid, non-polarizable molecularmodels from preceding work. Self- and Maxwell-Stefan diffusion coefficients aswell as shear viscosity were calculated by EMD and the Green-Kubo formalism.Comparing to experimental data for the two pure alcohols, present self-diffusioncoefficient and shear viscosity data are very good to excellent over a temperaturerange of about 100 K. Slightly higher deviations were found for the shear vis-cosity at low temperatures, which are mainly due to the long-time behavior ofthe shear viscosity autocorrelation function. Regarding the mixture methanol +ethanol at ambient conditions over the full composition range, both self-diffusioncoefficients also excellently agree with experimental data. The same holds forshear viscosity, especially in the methanol-rich region. Maxwell-Stefan diffusioncoefficient data is, to our knowledge, not available in the literature. Present pre-dictions correspond very well with the model of Darken. The thermal conductivitywas determined by BD reverse NEMD. Over a temperature range of about 100 K,an excellent agreement with the experiment was found as well for both pure fluids.Furthermore, predicted VLE data for the mixture methanol + ethanol along twoisotherms correspond with experimental data practically throughout within thesimulation uncertainty. The same holds for excess volume and excess enthalpy atambient conditions in the full composition range. It has been demonstrated thatpresent rigid, non-polarizable models for ethanol and methanol show a better per-formance to reproduce transport properties when compared to standard rigid and26exible models available in the literature.27 ppendix Simulation details
In this work, EMD simulations for transport properties were made in two steps.In the first step, a short simulation in the isobaric-isothermal (
N pT ) ensemble wasperformed at ambient pressure to calculate the density at the desired temperature.In the second step, a canonic (
NV T ) ensemble simulation was performed at thistemperature and density, to determine the transport properties. The simulationswere carried out in a cubic box with periodic boundary conditions containing2048 molecules. Newton’s equations of motion were solved using a fifth-orderGear predictor-corrector numerical integrator. The temperature was controlledby velocity scaling. In all EMD simulations, the integration time step was 0 . r c =
15 Å. Electrostatic long-range correctionswere made using the reaction field technique with conducting boundary conditions ( e RF = ¥ ) . On the basis of a center of mass cut-off scheme, the LJ long-rangeinteractions were corrected using angle averaging [ ? ]. The simulations were equi-librated in the NV T ensemble over 10 time steps, followed by production runsof 0 . × time steps. Diffusion coefficients and shear viscosity were cal-culated using Eqs. ( ?? ) to ( ?? ) with up to 9 900 independent time origins of theautocorrelation functions. The sampling length of the autocorrelation functionsvaried between 9 and 14 ps. The separation between time origins was chosen sothat all autocorrelation functions have decayed at least to 1 / e of their normalizedvalue to guarantee their time independence [ ? ]. The uncertainties of the predictedvalues were estimated using a block averaging method [ ? ].In the NEMD simulations with the PeX algorithm for the thermal conductiv-28ty, 800 molecules were placed in a parallelepiped box with L z = L x = L y ,where L i is the length of the simulation box in the different spatial directions. Pe-riodic boundary conditions were applied in all directions. The system was thenequilibrated over 1 ns at the desired temperature and pressure by N pT simulationusing a weak coupling bath [ ? ] with long-range corrections [ ? ] for pressure andenergy. The resulting density of the equilibration run was then employed to gen-erate a new set of simulations in order to develop the thermal gradient in a run of2 ns using the NEMD scheme. In this case an exchange frequency u = . − (i.e. every 150 time steps) for the energy flux was used. Electrostatic interactionswere treated with the reaction field with conducting boundary conditions. Thethermal conductivity was then computed as the average of 4 different runs over 1ns. The integration of the equations of motion was performed with a time step of2 fs with a cut-off radius of 12 Å. A Verlet neighbor list was employed to improvethe performance of the simulations.The VLE calculations of the mixture methanol + ethanol were performed us-ing the Grand Equilibrium method [ ? ] in combination with the gradual insertiontechnique [ ? ]. The latter, a Monte Carlo based expanded ensemble method, wasapplied to determine the chemical potential in the liquid phase. In comparison toWidom’s test particle method [ ? ], where a full molecule is inserted into the fluid,the gradual insertion method introduces a fluctuating molecule. This molecule un-dergoes changes in a set of discrete states of coupling with all other real moleculesof the fluid. Preferential sampling was performed in the vicinity of the fluctuatingmolecule to improve accuracy. The gradual insertion simulations were done with500 molecules in the liquid phase. Starting from a face-centered cubic lattice ar-29angement, the simulation runs were equilibrated over 2 000 cycles in the N pT ensemble. The production phase was performed over 30 000 Monte Carlo cycles.Further simulation parameters of these runs were taken from Vrabec et al. [ ? ].To evaluate the excess properties, EMD simulations in the N pT ensemblewere performed. The system contained 1372 molecules in a cubic box with pe-riodic boundary conditions. The molecular trajectories were sampled with thealgorithms and parameters as described above. The simulations were equilibratedin the
N pT ensemble over 80 000 time steps, followed by a production run over600 000 time steps. 30 eferences [1] Petravic, J.; Delhommele, J.
J. Chem. Phys. , 303.[2] Petravic, J.; Delhommele, J.
J. Chem. Phys. , 234509.[3] Petravic, J.
J. Chem. Phys. , 174503.[4] Wheeler, D. R.; Rowley, R. L.
Mol. Phys. , 555.[5] van de Ven-Lucassen, I. M. J. J.; Vlugt, T. J. H.; van der Zanden, A. J. J.;Kerkhof, P. J. A. M. Mol. Simul. , 79.[6] Hawlicka, E.; Swiatla-Wojcik, D. Phys. Chem. Chem. Phys. , 3175.[7] Nieto-Draghi, C. Ph.D. Thesis, Universitat Rovira i Virgili, Terragona, 2003.[8] Wensink, E. J. W.; Hoffmann, A.; van Mareen, P. J.; van der Spoel, D. J.Chem. Phys. , 37308.[9] Müller-Plathe, F.
Mol. Simul. , 133.[10] Zhang, L.; Wang, Q.; Liu, Y.-C.; Zhang, L.-Z. J. Chem. Phys. ,104502.[11] Noskov, S. Y.; Lamourex, G.; Roux, B.
J. Phys. Chem. B , 6705.[12] Dysthe, D. K.; Fuchs, A. H.; Rousseau, B.
J. Chem. Phys. , 4047.[13] Jorgensen, W. L.
J. Phys. Chem. , 1276.[14] Haughney, M.; Ferrario, M.; McDonald, I. R. J. Phys. Chem. , 4934.3115] Mayo, S. L.; Olafson, B. D.; Goddard III, W. A. J. Phys. Chem. ,8897.[16] van Leeuwen, M. E.; Smit, B. J. Phys. Chem. , 1831.[17] Gao, J.; Habibollazadeh, D.; Shao, L. J. Phys. Chem. , 16460.[18] Caldwell, J. W.; Kollman, P. A. J. Phys. Chem. , 6208.[19] Jorgensen, W. L.; Maxwell, D. S.; Tirado-Rives, J. J. Am. Chem. Soc. , 11225.[20] Potoff, J. J.; Errington, J. R.; Panagiotopoulos, A. Z.
Mol. Phys. ,1073.[21] Stubbs, J. M.; Chen, B.; Potoff, J. J.; Siepmann, J. I. Fluid. Phase Equil. , 301.[22] Chen, B.; Potoff, J. J.; Siepmann, J. I.
J. Phys. Chem. B , 3093.[23] Kahre, R.; Sum, A. K.; Nath, S. K.; de Pablo, J. J.
J. Phys. Chem. B , 10071.[24] Liem X. Dang, T. M. C.
J. Chem. Phys. , 9851.[25] Wang, S.; Cann, N.
J. Chem. Phys. , 214502.[26] Wheeler, D. R.; Fuller, N. G.; Rowley, R. L.
Mol. Phys. , 55.[27] Taylor, R. S.; Shields, R. L. J. Chem. Phys. , 12569.[28] Saiz, L.; Padró, J. A.; Guàrdia, E.
J. Phys. Chem. B , 78.3229] Patel, S.; Brooks III, C. L.
J. Chem. Phys. , 164502.[30] González, M. A.; Encisco, E.; Bermejo, F. J.; Bée, M.
J. Chem. Phys. , 8045.[31] Zhao, L.; Wang, X.; Wang, L.; Sun, H.
Fluid. Phase Equil. , 212.[32] Schnabel, T.; Srivastava, A.; Vrabec, J.; Hasse, H.
J. Phys. Chem. B , 9871.[33] Schnabel, T.; Vrabec, J.; Hasse, H.
Fluid Phase Equil. , 134.[34] Stoll, J.
Molecular Models for the Prediction of Thermophysical Propertiesof Pure Fluids and Mixtures ; Fortschritt-Berichte VDI, Reihe 3, 836; VDIVerlag: Düsseldorf, 2005.[35] Ungerer, P.; Beauais, C.; Delhommelle, J.; Boutin, A.; Fuchs, A. H.
J. Chem.Phys. , 5499.[36] Haslam, A. J.; Galindo A.; Jackson G.
Fluid Phase Equil. , 105.[37] Green, M. S.
J. Chem. Phys. , 398.[38] Kubo, R. J. Phys. Soc. Jpn. , 570.[39] Müller-Plathe, F. J. Chem. Phys. , 6082.[40] Darken, L.
AIME , 184.[41] Vignes, A.
Ind. Eng. Chem. Fund. , 189.[42] Caldwell, C. S.; Babb, A. L. Phys. Rev. A , 51.3343] Fernández, G. A.; Vrabec, J.; Hasse, H. Int. J. Thermophys. , 1389.[44] Hoheisel, C. Phys. Rep. , 111.[45] Gubbins, K. E.
Statistical Mechanics,
Vol. 1; The Chemical Society Burling-ton house: London, 1972.[46] Nieto-Draghi, C.; Avalos, J. B.
Mol. Phys. , 2303.[47] Pereira, J. C. G.; Catlow, C. R. A.; Price, G. D.
J. Phys. Chem. A ,1909.[48] Meier, K.; Laesecke, A.; Kabelac, S.
J. Chem. Phys. , 3671.[49] Alonso, J.; Bermejo, F. J.; García-Hernández, M.; Martínez, J. L.; Howells,W. S.
J. Mol. Struct. , 147.[50] Michels, P. J.; Trappeniers, N. J.
Chem. Phys. Lett. , 195.[51] Sindzingre, P.; Klein, M. L. J. Chem. Phys. , 4681.[52] Alder, B. J.; Gass, D. M.; Wainwright, T. J. Chem. Phys. , 3813.[53] Luo, H.; Hoheisel, C. J. Chem. Phys. , 8378.[54] Luo, H.; Hoheisel, C. Phys. Rev. A , 6421.[55] Nieto-Draghi, C.; Avalos, J. B; Rosseau, B. J. Chem. Phys. , 4782.[56] Karger, N.; Vardag, T.; Lüdemann, H. D.
J. Chem. Phys. , 3437.[57] Butcher, K. L.; Robinson, W. I. J. Appl. Chem. , 289.3458] Peng, D.-Y.; Robinson, D. B. Ind. Eng. Chem. Fund. , 59.[59] Benson, G. C.; Pflug, H. D. J. Chem. Eng. Data , 382.[60] Ogawa, H.; Murakami, S. J. Solution Chem. , 315.[61] Zarei, H. A.; Jalili, F.; Assadi, S. J. Chem. Eng. Data , 2517.[62] Pflug, H. D.; Pope, A. E.; Benson, G. C. J. Chem. Eng. Data , 408.[63] Lustig, R. Mol. Phys. , 175.[64] Schoen, M.; Hoheisel, C. Mol. Phys. , 33.[65] Allen, M. P.; Tidesley, D. J. Computer Simulation of Liquids,
J. Chem. Phys. , 3684.[67] Vrabec, J.; Hasse, H. Mol. Phys. , 3375.[68] Nezbeda, I.; Kolafa, J.
Mol. Simul. , 391.[69] Widom, B. J. Chem. Phys. , 2808.[70] Vrabec, J.; Kettler, M.; Hasse, H. Chem. Phys. Lett. , 431.[71] Casulleras, J.; Guardia, E.
Mol. Simul. , 273.[72] Asahi, N.; Nakamura, Y. J. Chem. Phys. , 9879.[73] Johnson, P. A.; Babb, A. L.
Chem. Rev. , 587.3574] Dullien, F. A. L. AIChE J. , 62.[75] Hurle, R. L.; Eastel, A. J.; Woolf, L. A. J. Chem. Soc. Faraday Trans. 1 , 769.[76] Meckl, S.; Zeidler, M. D. Mol. Phys. , 85.[77] Johnson, P. A.; Babb, A. L. J. Chem. Phys. , 14.[78] Rauf, M. A.; Steward, G. H.; Farhataziz. J. Chem. Eng. Data , 324.[79] Vargaftik, N. B.; Vinogradov, Y. K.; Yargin, V. S. Handbook of physicalProperties of liquids and Gases. Pure Substances and Mixtures,
J. Chem. Eng. Data , 46.[81] Poling, B. E.; Prausnitz, J. M.; O ′ Connell, J. P
The Properties of Gases andLiquids,
Thermal Conductivity. Nonmetallic Liq-uids and Gases;
Thermophysical Properties of Matter Vol.3; Plenum: NewYork, 1970.[83] Lide, D. R., Ed.
CRC Handbook of Chemistry and Physics, ? , ? , ? , ? , ? ] [ ? ] [ ? ]Ethanol [ ? , ? , ? , ? , ? , ? , ? ] [ ? , ? ] [ ? , ? ]Methanol + Ethanol − − − ?? ). Boltzmanns’s constant is k B andthe electronic charge is e .Site s aa e aa / k B q ia Å K e MethanolS
CH3 + OH − H − − + CH3 − S CH2 + OH − H − − + . T r D i h l K mol l − − m − s − − Pa s Wm − K − Methanol213 27.02 0.318 (5) 3.1 (4)220 26.80 0.433 (6) 2.0 (2) 0.22 (3)240 26.22 0.750 (8) 1.7 (2)260 25.63 1.21 (1) 1.0 (2) 0.207 (9)278.15 25.14 1.72 (1) 0.64 (8)288 24.81 2.09 (2) 0.62 (7)298.15 24.51 2.41 (2) 0.57 (6)300 24.50 0.191 (9)308.15 24.19 2.97 (2) 0.47 (5)318.15 23.91 3.54 (2) 0.41 (6)320 23.88 0.188 (8)328.15 23.58 4.15 (3) 0.42 (5)340.15 23.21 4.91 (3) 0.29 (5)Ethanol239 18.11 0.192 (6) 3.1 (3)253.15 17.99 0.250 (5) 2.4 (2) 0.18 (1)263 17.83 0.363 (7) 2.2 (3)273.15 17.62 0.493 (8) 1.7 (2) 0.178 (9)280 17.46 0.64 (1) 1.4 (2)298.15 17.08 1.07 (1) 1.1 (1) 0.171 (9)308.15 16.88 1.33 (2) 1.0 (2)320 16.66 1.68 (1) 0.8 (1)328.15 16.45 2.01 (2) 0.61 (8) 0.161 (8)333 16.31 2.20 (2) 0.60 (7)338.15 16.26 2.40 (2) 0.47 (7)39able 4: Density, self-diffusion and Maxwell-Stefan diffusion coefficient as wellas shear viscosity of the mixture methanol (1) + ethanol (2) at 298.15 K and 0 . x r D D − D h mol mol − mol l − − m − s − − m − s − − m − s − − Pa s0 27.02 1.07 (1)0.143 26.22 1.37 (2) 1.12 (1) 1.4 (3) 1.17 (9)0.305 25.63 1.52 (1) 1.25 (1) 1.8 (3) 0.76 (9)0.500 25.14 1.73 (1) 1.44 (1) 1.7 (2) 0.62 (8)0.570 24.81 1.82 (1) 1.53 (1) 1.6 (3) 0.64 (9)0.754 24.51 2.07 (1) 1.72 (2) 2.3 (3) 0.70 (9)0.796 23.91 2.14 (2) 1.82 (2) 2.0 (3) 0.59 (7)0.890 23.58 2.28 (1) 1.90 (2) 2.6 (3) 0.60 (8)1 23.21 4.95 (3) 40able 5: Vapor-liquid equilibrium data for the mixture methanol (1) + ethanol(2) at 393.15 K and 433.15 K from present molecular dynamics simulation. Thenumber in parenthesis indicates the statistical uncertainty in the last digit. x p y r ′ r ′′ D h v mol mol − MPa mol mol − mol l − mol l − kJ mol − T = 393.15 K0 0.551 (6) 0 14.84 (2) 0.200 (1) 33.24 (7)0.150 0.448 (9) 0.207 (5) 15.62 (2) 0.154 (3) 33.86 (6)0.282 0.482 (8) 0.369 (2) 16.35 (2) 0.168 (3) 33.32 (7)0.492 0.538 (9) 0.589 (2) 17.54 (3) 0.185 (3) 32.67 (8)0.678 0.588 (9) 0.752 (5) 18.87 (3) 0.212 (3) 31.92 (9)0.838 0.572 (9) 0.884 (3) 19.94 (3) 0.205 (3) 31.43 (8)1.0 0.767 (5) 1.0 21.30 (3) 0.294 (1) 29.74 (8) T = 433.15 K0 1.32 (3) 0 13.58 (4) 0.469 (2) 28.41 (7)0.150 1.29 (2) 0.192 (4) 14.34 (4) 0.453 (7) 28.52 (9)0.282 1.38 (2) 0.339 (5) 14.96 (4) 0.486 (7) 28.08 (9)0.488 1.48 (1) 0.549 (5) 16.09 (3) 0.529 (5) 27.50 (8)0.588 1.54 (2) 0.659 (5) 18.87 (3) 0.548 (7) 27.28 (8)0.674 1.58 (2) 0.734 (5) 17.22 (4) 0.560 (7) 27.07 (9)0.836 1.64 (2) 0.872 (2) 18.31 (4) 0.594 (7) 26.55 (8)0.942 1.72 (2) 0.957 (1) 19.15 (5) 0.633 (7) 26.12 (9)1.0 1.89 (3) 1.0 19.42 (4) 0.711 (2) 26.12 (9)41able 6: Volume, excess volume, enthalpy and excess enthalpy of the mixturemethanol (1) + ethanol (2) at 298 .
15 K and 0 . x v v E h h E mol mol − cm − mol − cm − mol − kJ mol − J mol − − .
46 (4)0.070 57.15 (5) 0.00 (7) − .
20 (4) 3 (52)0.200 54.94 (5) 0.01 (6) − .
75 (4) 2 (49)0.270 53.72 (5) 0.03 (6) − .
50 (4) − − .
22 (4) − − .
95 (4) 20 (45)0.500 49.65 (4) 0.02 (5) − .
69 (4) − − .
32 (4) 5 (44)0.697 46.19 (4) 0.04 (5) − .
96 (4) 20 (44)0.800 44.34 (4) 0.00 (5) − .
61 (4) 1 (45)0.905 42.48 (4) 0.00 (5) − .
24 (4) − − .
90 (3)42 ist of Figuresist of Figures