MeOH-4P - American Chemical Society

Article Cite This: J. Chem. Theory Comput. XXXX, XXX, XXX−XXX

pubs.acs.org/JCTC

A Four-Site Molecular Model for Simulations of Liquid Methanol and Water−Methanol Mixtures: MeOH-4P Manuel Martínez-Jiménez and Humberto Saint-Martin* Instituto de Ciencias Físicas, Universidad Nacional Autónoma de México, Apartado Postal 48-3, Cuernavaca 25510, México S Supporting Information *

ABSTRACT: In this work, we present a new four-site potential for methanol, MeOH-4P, fitted to reproduce the dielectric constant ε, the surface tension γs, and the liquid density ρ of the pure liquid at T = 298.15 K and p = 1 bar. The partial charges on each site were taken from the OPLS/2016 model with the only difference of putting the negative charge on the fourth site (M) instead of on the O atom, as done in four-site water models. The original Lennard-Jones (LJ) parameters of OPLS/2016 for the methyl moiety (Me) were modified for the fitting of ρ and γs, whereas the parameters of the TIP4P-FB water model were used for the O atom without change. Taking into account the energetic cost of the enhanced dipole relative to the isolated molecule, the results from simulations with this model showed good agreement with experiments for ρ, αp, κT, Cp, and ΔHv−l. Also, the temperature dependence of γs and ε is satisfactory in the interval between 260 and 360 K, and the critical point description is similar to that of OPLS/2016. It is shown that orientational correlations, described by the Kirkwood factor Gk, play a prominent role in the appropriate description of dielectric constants in existing models; unfortunately, the enhancement of the dipole moment produced a low diffusion coefficient DMeOH; thus, a compromise was required between a good reproduction of ε and an acceptable DMeOH. The use of a fourth site resulted in a significant improvement for water−methanol mixtures described with TIP4P-FB and MeOH-4P, respectively, but required the modification of the LJ geometric combination rule to allow a good description of the methanol molar-fraction dependence of ρ, ε, and methanol (water) diffusion coefficients DMeOH (DH2O) and excess volume of mixing ΔVmix in the entire range of composition. The resulting free energy of hydration ΔGhyd shows excellent agreement with experiments in the interval between 280 and 360 K.

1. INTRODUCTION One of the challenges of (bio)molecular simulation is the formulation of realistic potential energy functions describing molecular interactions in the condensed phase with accurate force-field parameters. The dielectric constant has particular relevance in solubility processes and for the appropriate description of the separation of liquid phases in mixtures having polar components. Unfortunately, nonpolarizable force fields are known to have difficulties in reproducing this property in the liquid phase. Caleman et al.1 calculated several properties for 146 organic liquids using OPLS/AA and GAFF molecular models, including the density ρ, enthalpy of vaporization ΔHv−l, surface tension γs, heat capacity at constant volume Cv, and dielectric constant ε under ambient conditions of pressure and temperature. For both force fields, most of the analyzed properties showed a reasonable agreement with experiments, but serious issues concerning γs and ε were also evident. A subsequent computation2 of γs values for a smaller set of 61 liquids, this time taking into account proper treatment of long-range interactions, yielded a better agreement with experimental results, but the systematic underestimation of ε, for most liquids, still remained. In the case of water3 and methanol,4 it has been shown that most used models systematically underestimate the value of ε. A plausible rationale for nonsatisfactory values of this property, as given by rigid nonpolarizable models, relies on the fact that point© XXXX American Chemical Society

charge distributions could have intrinsic limitations to reproduce, simultaneously, quantum-mechanical potential energy and dipole moment surfaces (at least two- and threebody). Most models have been obtained from fitting procedures based on energetic criteria. In spite of these considerations, it has been shown that it is possible to reproduce the experimental value of ε for water with nonpolarizable potentials,5−7 being that these models are also able to reproduce several thermodynamic and transport properties such as the temperature of maximum density and the self-diffusion constant. Recently, Salas et al.8 used the same parametrization strategy as in ref 6 with ε, γs, and liquid ρ as target properties; this procedure was applied to pyridine, dichloromethane, methanol, and 1-ethyl-3-methylimidazolium tetrafluoroborate (EMIM-BF4) molecules from the all-atom to coarse-grained level. A three-site model was used for methanol; the temperature dependence of ε, γs, and ρ and the critical point predictions were in good agreement with experiments, but the model yielded a low self-diffusion coefficient and a highly structured liquid. Interestingly, to our present knowledge, none of these potentials fitted to reproduce ε have been used in numerical simulations of water−methanol mixtures, being that these are a stringent test of the transferability of Received: December 18, 2017 Published: March 22, 2018 A

DOI: 10.1021/acs.jctc.7b01265 J. Chem. Theory Comput. XXXX, XXX, XXX−XXX

Article

Journal of Chemical Theory and Computation

= 298.15 K; this was done by a linear scaling of the geometric combination rule between the oxygen in water (OH2O) and the oxygen and the methyl group in methanol (OMeOH, Me, respectively)

analytical potentials which have been optimized, individually, for the liquid and/or solid phase of each component. Because of their nonideal physical behavior and practical applications, water−methanol liquid mixtures have received considerable attention in literature from pioneering works9−11 considering single water molecules in methanol to extensive studies12−18 of structural, dynamic, and thermodynamic properties in the entire range of composition. Recently, Galicia-Andrés et al.19 reported a comparative analysis of OPLS-AA20 and OPLS-UA21 methanol models, coupled with SPC-E22 and TIP4P-Ew23 water models, for several properties including liquid densities ρ, dielectric constant ε, excess volume of mixing ΔVmix, excess enthalpy of mixing ΔHmix, and water (methanol) diffusion coefficients DH2O (DMeOH) in the entire range of methanol composition. Their results showed reasonable agreement of liquid ρ and ΔVmix with experiments and a qualitative reproduction of DH2O (DMeOH). However, these models were clearly unsuccessful in reproducing ε and ΔHmix as a function of the methanol molar fraction, and the authors pointed out the need for a refined parametrization of charges with ad hoc LJ combination rules for the interactions between water and methanol. About this, Moučka and Nezbeda24 implemented small deviations from the Lorentz− Berthelot rule for these cross-interactions and analyzed their effect on the properties of the mixtures using the nonpolarizable potentials OPLS-UA (methanol) and TIP4P25 (water), keeping the geometric LJ rule intact for the pure fluid of methanol. This approach provided a qualitative improvement of ΔHmix and a better agreement of ΔVmix with experiments, being also clear that this strategy could not lead, by itself, to a dramatic improvement of mixture properties. A natural choice in trying to overcome present limitations is the use of polarizable models combined with optimized LJ combination rules. While further work in this direction is clearly imperative, it still remains of interest to understand to what extent the performance of nonpolarizable models can be improved over current benchmarks, particularly when they are applied to the study of multicomponent mixtures. For this reason, we decided to take advantage of recent advances in nonpolarizable models to engineer a new analytical potential for methanol, hereafter denoted MeOH-4P, using ε, γs, and ρ of liquid methanol as target quantities under ambient conditions of pressure and temperature. The key differences with respect to existing potentials are the addition of a fourth interaction site in methanol monomer-geometry and the transferability of LJ hydroxyl-group parameters from the TIP4P-FB water model.7 Several thermodynamic properties, namely, ρ, αp, κt, Cp, ΔHv−l, the temperature-dependence of ε, and the critical point were calculated to assess the performance of the new model. A comparative analysis of the orientational correlations, measured by the finite Kirkwood factor Gk, was done to better understand the issues for the adequate reproduction of ε in liquid phase and get further insight into new routes of improvement. Finally, we used the TIP4P/ε and TIP4P-FB water models combined with MeOH-4P to investigate water−methanol mixtures, incorporating a modified combination rule for LJ crossinteractions between water and methanol and observing their effect in mixing properties such as ρ, ΔVmix, ΔHmix, ε, and DH2O (DMeOH) diffusion coefficients, as functions of the methanol mix molar fraction. The new εmix a−b and σa−b LJ parameters for crossinteractions between water and methanol molecules were optimized using ρ and ΔHmix as target quantities at p = 1 and T

εamix − b = λ1 εaεb ;

σamix − b = λ 2 σaσb ;

a ∈ {OH2O}, b ∈ {OMeOH , Me}

(1)

while the LJ parameters between molecules of the same type were kept identical as in pure water (methanol).

2. SIMULATION DETAILS The GROMACS 5.0.426 package was used for all the simulations in this work, and the equations of motion were integrated using the leapfrog algorithm, with a time step of 2 fs. Liquid phase simulations were performed in the isotropic NPT ensemble, under ambient temperature (T = 298.15 K), pressure (p = 1 bar), and a fixed number of molecules N = 435. The LJ and the real part of the electrostatic interactions were truncated at 1 nm, and the PME method was used for the long-range electrostatic part with the following settings: tolerance of 10−5 for the real space contribution, grid spacing of 0.12 nm, and 4order interpolation. The energy and pressure corrections27 implemented in GROMACS 5.0.4 were also applied due to the use of a finite cutoff for LJ interactions; dielectric constant values were computed from the system dipole-moment fluctuations. Coexistence properties were obtained from NVT simulations, using an explicit liquid−vapor interface containing N = 2250 molecules with T in the interval between 250 and 500 K. This interface was generated by setting up a liquid slab surrounded by vacuum in a simulation box with periodic boundary conditions in the three spatial directions; the dimensions of the simulation cell were Lx = Ly = 6 nm with Lz = 3Lx, with z the normal direction to the liquid−vapor interface. The PME method was used for an electrostatic longrange contribution with a reciprocal grid of 0.12 nm (x and y directions) and 0.3 nm (z direction). A 2.6 nm cutoff was used for the LJ and the real electrostatic contributions, with a tolerance of 10−5, and long-range dispersion corrections for energy and pressure were turned off. In liquid mixtures, the LJ cross-interactions between water and methanol molecules were computed using modified combination rules (the parameters are given in Models section). Temperature coupling was implemented using the Nosé−Hoover thermostat (τT = 1.0 ps), and the isotropic Berendsen barostat was used for pressure control (τP = 0.6 ps). The Bennett acceptance ratio (BAR) method that was used for free energy calculations with the GROMACS 5.0.4 code required the stochastic dynamics integrator instead of the leapfrog scheme. Bond distances were kept rigid with the LINCS algorithm. Quantum mechanical calculations were performed using the Gaussian0928 package: methanol geometric optimizations, population analysis, and partial charge determinations were obtained using the basis sets aug-cc-PVxZ (x = D, T, Q, and 5) with the MP2 electronic-correlation treatment. The dielectric constant with periodic boundary conditions was calculated from ε=1+

4π 4π ⟨M2⟩ = 1 + Nμ2 Gk 3kTV 3kTV

(2)

where V and M denote the volume and the total dipole moment of the simulation box. Gk is the finite system Kirkwood B


Article


Table 1. Partial Charges, Lennard-Jones, and Structural Parameters of Methanol Models MD1, MD2, TraPPE-UA, OPLS-UA, OPLS/2016, and MeOH-4Pa Model

εOO /K

σOO /Å

εMeMe /K

σMeMe /Å

qO /e

qH /e

qMe /e

qM /e

lOM /Å

MD1 MD2 TraPPE-UA OPLS-UA OPLS/2016 MeOH-4P TIP4P-FB TIP4P/ε TIP4P/2005 TIP4P-Ew

83.710 88.345 92.995 85.547 97.775 90.118 90.118 93.000 93.200 81.954

3.0200 3.0260 3.0200 3.0700 3.1659 3.1655 3.1655 3.1650 3.1589 3.16435

88.208 93.098 97.998 104.166 110.450 106.032 − − − −

3.7500 3.7580 3.7500 3.7750 3.6499 3.6370 − − − −

−0.7490 0.0000 −0.7000 −0.7000 −0.6544 0.0000 0.0000 0.0000 0.0000 0.0000

0.4650 0.3844 0.4350 0.4350 0.4998 0.4998 0.52587 0.5270 0.5564 0.52422

0.2840 0.2474 0.2650 0.2650 0.1546 0.1546 − − − −

− −0.6318 − − − −0.6544 −1.05174 −1.0540 −1.1128 −1.04844

− −0.1500 − − − −0.037 0.10527 0.1050 0.1546 0.1250

Water models are TIP4P-FB,7 TIP4P/ε, TIP4P/2005,50 and TIP4P-Ew.23 Values of the ϵ are given in units of kB. All methanol models in this table use rO−H = 0.945 Å, rO−Me = 1.43 Å, and θ = 108.5°. TraPPE-UA and MD1 flexural constant is κ = 460.67 kJ mol−1 rad2. All water models in this table use rO−H = 0.9572 Å and θ = 104.52°.

a

obtained with the pressure normal to the liquid−vapor interface, and data was fitted to the Clausius−Clapeyron equation

factor related to the relative orientation of individual molecular dipoles Gk = 1 +

2 N

N

∑ ui·uj

log(P /P0) = a − b/T (3)

i>j

and P0 was chosen to be 1 bar. The excess enthalpy of mixing ΔHmix and excess volume of mixing ΔVmix were obtained from equations

We must emphasize the fact that long simulations, at least 100 ns in this work, are required to get reliable statistics from models having enhanced dipole moments because of the large fluctuations in the dynamics. The surface tension γs was computed using an explicit interface through the diagonal elements Pxx, Pyy, and Pzz of the pressure tensor and the length Lz of the box which is perpendicular to the interface γs =

⟨Pxx⟩ + ⟨Pyy⟩ ⎤ Lz ⎡ ⎢⟨Pzz⟩ − ⎥ 2⎣ 2 ⎦

ΔVmix = V − (1 − X m)Vw − X mVm ΔHmix = H − (1 − X m)Hw − X mHm

(Self-)diffusion coefficients were estimated in pure methanol and water−methanol mixtures by using the Einstein equation

αp =

(5)

t →∞

(ρl − ρg )/ρc = A 0|τ | + A1|τ |

βc +Δ

+ A 2 |τ |

βc + 2Δ

+ ···

⎛ ρ[290K] − ρ[306K] ⎞ 1 ⎟ ⎜ ⎠ 16 ρ[298.15K] ⎝

⎛ ∂ log ρ ⎞ ⎟ κT = ⎜ ⎝ ∂P ⎠T

Cp =

(11)

⎛ ∂H ⎞ ⎜ ⎟ ⎝ ∂T ⎠ P

(12)

(13)

3. MODELS Inspired by the well-known improved performance of four-site water models, we investigated the consequences of incorporating an additional massless site M into the charge distribution of the methanol monomer, using the dielectric constant, surface tension, and density at 298.15 K and p = 1 bar as target quantities. The new potential is denoted as MeOH-4P, and its parametrization is based on those of the TIP4P/ε6 and TIP4PFB7 water models and on that of the recent OPLS/20164 model for methanol, which was fitted to reproduce liquid properties and solid−liquid and vapor−liquid equilibria. To carry out a comparative analysis with several existing potentials, we also designed a second model, denoted as MD2, as it was parametrized following the strategy used for the three-site MD1 in ref 8; thus, MD2 differs from MeOH-4P in the

(6)

and the critical density was deduced from the expansion of the diameter31 (ρl + ρg )/2ρc = 1 + D1 − α |τ |1 − α + D1|τ | + D2|τ |1 − α +Δ + D3|τ |1 − α + 2Δ ···

(10)

while isothermal compressibility κT and specific heat Cp are from the thermodynamics definitions

where ri(t) denotes the position of the i particle at time t. The highest dilute mixtures were simulated placing a single molecule of water (methanol) in the methanol (water) solvent box, and the diffusivity of each molecule was calculated as a function of the molar fraction of methanol. The system-size effects on the values of diffusion coefficients were taken into account following the method reported by Yeh and Hummer.29 The critical temperature was determined by the Wegner expansion30 of the order parameter βc

(9)

where the enthalpy of the system H was estimated by Upot + PV. Following previous work,34 the isobaric thermal expansivity αp at 298.15 K was obtained by using the incremental method through the equation

(4)

6Dt = lim ⟨ri(t ) − ri(0)⟩2

(8)

(7)

where ρl,ρg are the liquid and vapor densities, respectively, and T τ is given by τ = 1 − T with βc = 0.325 and Δ = 0.5 being c

scaling parameters according to the 3D Ising universality class,32,33 and α = 0.11 is the critical exponent for the weak singularity of the heat capacity. The pressure of the system was C


Article


4. RESULTS 4.1. Liquid Methanol. Figure 1 shows the satisfactory reproduction of the experimental dielectric constant as a

determination method for the partial charges, as described in the following paragraphs. The inclusion of a fourth site in water models has been relatively simple, along the bisector of the angle H−O−H, but the charge-asymmetry between methyl and hydroxyl groups in methanol requires a different approach. In this work, we propose a simple yet reasonably realistic modeling by obtaining the induced dipole moment μ⃗ PCM from the wave function of the methanol monomer in a polarizable continuum, using the integral equation formalism (IEFPCM) in Gaussian09. The orientation of μ⃗PCM resulted in a shift of 1.78° with respect to the H−O−CH3 angle bisector, and we used the ray along this direction to place the M site, with the distance lOM chosen to match the dielectric constant, together with the partial charges. Unlike water models, the position of the M site was found to be in the negative direction of the OM ray, where the positive direction is defined by the bond angle CH3−O−H. In both models, the hydroxyl atoms were explicitly modeled, and the united-atom approach was used for the methyl group. For MD2, geometry optimizations were performed in Gaussian09 with the MP2 level of theory and the aug-ccPVxZ (x = D, T, Q, and 5) basis sets to ensure convergence. Next, we describe the fitting procedure for each model, and their parameters are shown in Table 1: • MeOH-4P: Atomic partial charges were the same as those in OPLS/2016, only turning the oxygen atom chargeless and assigning its charge value to the new massless M virtual site. Performing several simulations, the location of M along μ⃗ PCM was varied until the reproduction of the dielectric constant was satisfactory. Finally, methyl-group LJ values ε and σ were adjusted to reproduce the surface tension and the liquid density, keeping the same parameters of the TIP4P-FB7 water model for the oxygen atom. • MD2: Atomic partial charges for hydroxyl atoms and the methyl group were calculated with ChelpG analysis from single-methanol optimization in a polarizable continuum, using the integral equation formalism (IEFPCM); in this way, polarization effects in liquid phase are reasonably taken into account. These partial charges were scaled up, and the location of M along μ⃗PCM was simultaneously varied until the reproduction of the dielectric constant was satisfactory. The Lennard-Jones parameters of both the methyl group and oxygen were fitted to reproduce the surface tension and the liquid density. The values of the parameters of MeOH-4P, from eq 1, are the following: λ1 = 0.915;

λ 2 = 0.996

Figure 1. Dielectric constant as a function of temperature calculated for MeOH-4P and existent methanol models. Experimental data are from refs 53 and 54.

function of temperature in the interval between 260 and 360 K for the two new models presented in this work: MeOH-4P and MD2. Compared to the predictions of OPLS-UA, OPLS/2016, and MD18 three-site models, the curve slope also shows a better agreement with experiments in temperatures other than 298.15 K. With ε a target quantity in this work, this fact should not be surprising, as there are more interesting insights about the role of the effective dipole moment μmod and the Kirkwood factor Gk in the appropriate reproduction of ε that can be extracted from a simple comparative analysis considering several existing models (Table 2). In Figure 2, three regions are distinguished according to the values of the dielectric constant. In the first group, models L1,35 L2,36,37 H1,38,39 and OPLS-UA are agglutinated around ε = 22, and although the TraPPE-UA40 model has a slightly higher epsilon, it can also be included in this group. The parameters of all these models were adjusted to reproduce liquid phase and/ or liquid−vapor coexistence properties, and the absolute difference between Gk and μmod is very similar for them. The second region, about ε = 27, includes only the model OPLS/ 2016 whose absolute difference between Gk and μmod is notably larger than the same value calculated for all models in the first group, and the parametrization of this model included liquid, liquid−vapor, and solid−liquid equilibrium properties. Finally, the third group comprises MD1, MD2, and MeOH-4P models, which were intentionally adjusted to reproduce the experimental ε and also show large Gk values. Each of them can be considered a reparameterization of some model belonging to the previous two groups: MeOH-4P comes from OPLS/2016, while MD1 and MD2 are based on TraPPE-UA. For these potentials, the two procedures presented in this work to match experimental ε were found not to affect significantly the correlation between μmod and Gk observed in their parent models. When dealing with rigid nonpolarizable models having fixed point charges, it is clear that their ability for reproducing the experimental dielectric constant relies on their effective dipole moment μmod and the Kirkwood factor Gk (see eq 2). The first

(14)

They were obtained from simulations of water−methanol liquid mixtures using the TIP4P-FB and TIP4P/ε water models. While the fitting of mixing density was straightforward, the reproduction of experimental excess enthalpy of mixing, in the entire range of composition, turned out to be challenging. In this case, the main criterion for optimization was not to exceed the minimum experimental value for this property. The final parameters reflect a compromise between the accurate reproduction of the selected target quantities and the sensitivity of each one to the adjustment of parameters in the fitting; for this reason, a maximum deviation of 5% from the experimental value was considered acceptable for each target property. D


Article


Table 2. Comparative of Kirkwood Factor Gk, Dipole μmod, and Quadrupole Principal Components Qxx, Qyy, and Qzz of Models Analyzed in This Worka Model

μmod /D

Gk

ε

Qxx /DÅ

Qyy /DÅ

Qzz /DÅ

QT

|1 − (μmod/QT)|

H1 L1 L2 TraPPE-UA MD1 MD2 OPLS-UA OPLS/2016 MeOH-4P Gas(QMb) Gas(QMc)

2.336 2.221 2.147 2.221 2.432 2.481 2.221 2.178 2.302 1.699 1.701

2.837 3.039 3.021 3.064 3.044 3.283 2.939 3.620 3.931 − −

21.5 22.8 21.2 24.2 30.5 31.6 21.6 26.8 31.4 − −

2.671 2.687 2.660 2.687 2.873 2.631 2.687 3.050 3.110 3.080 3.081

−1.442 −1.437 −1.417 −1.437 −1.537 −1.733 −1.437 −1.575 −1.691 −1.863 −1.874

−1.229 −1.250 −1.243 −1.250 −1.336 −0.858 −1.250 −1.475 −1.419 −1.216 −1.207

2.056 2.062 2.039 2.062 2.205 2.202 2.062 2.312 2.400 2.478 2.472

0.136 0.077 0.053 0.077 0.103 0.127 0.077 0.058 0.041 − −

a

Quantum-mechanics (QM) calculations were taken from the NIST-CCCB database. (Pseudo)atoms O, H, and Me lie in the plane xy, and the z axis is perpendicular to this plane following a right-handed coordinate system. Experimental values of ε have been reported to be in the interval from 31.5 to 33.0.51,52 Quadrupole magnitude is QT = (Qxx − Qyy)/2, as in previous analysis for water models reported in the literature.3,43. bMP2 = Full/ aug-cc-pVTZ. cMP2 = Full/aug-cc-pVQZ.

However, the most widely used methods to assign partial charges in empirical models are based on reproducing the electrostatic potential produced by the molecule on its vicinity, the CHELPG-type42 methods. The second quantity, Gk, takes into account the short-range orientation correlations between a given dipole and its surrounding neighbors, being expected to play a prominent role when simulating polar liquids. However, the difficulty of nonpolarizable models to reproduce the experimental dielectric constant has been attributed, mainly, to the absence of explicit polarizability, without examining in more detail the differences about orientation correlations provided by existing models. From Figure 2, it is evident that the improved performance of OPLS/2016 about ε is due to the fact that its Kirkwood factor is significantly higher than its counterparts in previous models, even though its effective dipole moment μmod is one of the smallest among them. This fact is clearly a consequence of its robust empirical fitting, which led to an optimal set of partial charges and LJ parameters able to improve the description of orientation correlations over all existent models, including those whose parameters were obtained from high-quality ab initio calculations. It is worth noticing that the partial charges of OPLS/2016 are not a simple rescaling of those from OPLS-UA (Table 1); thus, they do not correspond to a CHELPG-type procedure. Figure 2 also shows that the same considerations also apply to MeOH-4P with respect to other models with an enhanced μmod. Following the analysis of Abascal and Vega43 with four-site water models, we calculated the principal components of the quadrupole tensor of single methanol for each model and from ab initio calculations in gas phase (Table 2). It was found that molecular

Figure 2. Comparison of the values of the effective dipole moment μmod (triangle) and Kirkwood factors Gk (circles) as given by the models analyzed in this work. For each model, the magnitude of their relative difference can be appreciated using the horizontal line joining them.

value, μmod, is completely determined by the charge distribution (geometry and partial charges) fixed during the simulation. The main problem of assigning partial charges to the atoms in a molecule is that they are not observables of the system (the molecule); that is, they can neither be measured nor computed unequivocally. Several different population analyses have been designed, based on charge densities obtained from quantum calculations (e.g., the “atoms-in-molecules” approach, AIM41).

Table 3. Comparison of Values for Several Properties of Liquid Methanol As Predicted by Models Analyzed in This Work with Experimental Data at 298.15 K (and p = 1 bar when appropriate)a Model

ρ /kg m−3

αp /(kK)−1

κT /(TPa)−1

Cp J /(mol K)−1

ε

ΔHv−l /kJ mol−1

γs /mN m−1

DMeOH /cm2 s−1

MD1 MD2 OPLS-UA OPLS/2016 MeOH-4P Exp.

786 786 762 785 786 78655

1.30 1.30 1.33 1.28 1.27 1.2055

1180 1210 1302 1236 1231 124856

83.0 87.8 86.0 86.0 86.5 81.557

30.5 31.6 21.6 26.8 31.4 31.5−33.051,52

39.72 39.72 34.91 37.55 39.73 37.4358

22.5 22.5 22.7 24.7 23.3 22.5159

0.60 0.90 2.11 2.72 2.00 2.4460

Here, ρ denotes density, αp isobaric thermal expansivity, κT isothermal compressibility, Cp isobaric heat capacity, ε dielectric constant, ΔHv−l vaporization enthalpy (including polarization correction), γs surface tension, and DMeOH self-diffusion coefficient. a

E


Article

Journal of Chemical Theory and Computation quadrupoles of the best models, MeOH-4P and OPLS/2016, are very close to that of the isolated molecule obtained with quantum-mechanics calculations; this feature was also observed in water models having top performance. In the mentioned study, a quotient μmod/QT ≈ 1 was shown to be necessary for the relative stability between different ice phases; the respective values of MeOH-4P and OPLS/2016 are also very close to unity. This not being fortuitous as their partial charges came from fitting involving solid−liquid equilibria and the density ρ of the β solid at 170 K. The following thermodynamics properties were calculated for the new model: ρ, αp, κT, Cp, ΔHv−l, and γs at ambient conditions. It is observed, from Table 3, that OPLS/2016 original predictions are not significantly deteriorated. The enthalpy of vaporization ΔHv−l was corrected considering the energetic cost of molecular polarization22 from the gas phase to the liquid phase Epol =

Table 4. Critical Temperature, Critical Pressure, and Critical Density for Models Analyzed in This Work Compared to Experimental Measurements Tc/K

pc/bar

ρc/kg m−3

L1 L2 MD1 MD2 OPLS-UA OPLS/2016 MeOH-4P Exp.

505.0 510.4 533.0 525.5 490.0 521.5 521.9 512.661

77.0 80.0 76.0 74.8 70.0 75.0 88.2 81.061

269 273 251.2 252 265 272 255 27258/26961

specifically optimized for the vapor−liquid coexistence). A common feature of the new models is the appreciable underestimation of vapor density and, to a lesser extent, the overestimation of liquid density along the coexistence curve. This fact should be expected as their dipole moments were increased for matching ε, enhancing attractive interactions in the liquid, and consequently, the number of molecules in the vapor region decreased, shifting their critical points to higher Tc (lower ρc) values: 521.9 K (0.255 g cm−3) for MeOH-4P and 525.5 K (0.252 g cm−3) for MD2, with pc are 88.2 and 74.8 bar, respectively. When compared to experiments Tc = 512.6 K, ρc = 0.272/0.269 g cm−3, and pc = 81 bar, it is clear that the gasphase densities are underestimated, although a definitive conclusion should bear in mind that experimental values of critical density are affected by large uncertainties. MeOH-4P and MD2 performed rather well for the surface tension values, as a function of the temperature, in the interval between 250 and 500 K (Figure 4). We conclude that the vapor-pressure

(μmod − μgas )2 2αpol

Model

(15)

where μmod = 2.302 D is the effective dipole moment of MeOH-4P, and μgas = 1.70 D and αpol = 3.29 Å3 are the gas dipole moment and the molecular polarizability of methanol, respectively.44 This term is necessary for proper comparison with other models because of its physical basis and its large value: Epol = 3.06 kJ/mol for MeOH-4P and Epol = 7 kJ/mol for TIP4P/ε and TIP4P-FB water models, which also reproduce the dielectric constant at p = 1 bar and T = 298.15 K. The resulting value of ΔHv−l = 39.73 kJ/mol for MeOH-4P is higher than the experimental value ΔHv−l = 37.43 kJ/mol, although it does not overestimate it largely, a latent risk when the model has an enhanced dipole moment. The methanol self-diffusion coefficient DMeOH = 2.0 × 10−5 of MeOH-4P is better than the values DMeOH = 0.90 × 10−5 of MD1 and DMeOH = 0.60 × 10−5 of MD2, but still lower than the experimental value DMeOH = 2.44 × 10−5 cm2/s and the OPLS-2016 prediction DMeOH = 2.72 × 10−5 cm2/s. The predictions of MeOH-4P and MD2 models for liquid− vapor equilibria are also of interest; their orthobaric curves and critical points are shown in Figure 3 and Table 4 (L1 and L2 models are included, for comparison purposes, as they were

Figure 4. Temperature dependence of the surface tension of MeOH4P and existent methanol models. Experimental data are from ref 59.

(Figure 5) curve and the deviation of the critical point from experiments are similar to OPLS/2016, suggesting that the description of liquid−vapor equilibria by these models is reasonable, keeping in mind they were not particularly optimized for liquid−vapor equilibria and that the molecules evaporating from the liquid should gain the energy Epol. 4.2. Water−Methanol Liquid Mixtures. Thermodynamic excess properties measure the nonideal behavior of multicomponent systems, being a relevant test for analytic potentials in the simulation of liquid mixtures. For nonpolarizable models, the absence of local electrostatic response will lead to inherent

Figure 3. Orthobaric densities and critical points of MeOH-4P and existent methanol models. Experimental data are from ref 61. F


Article


XMeOH < 0.9, and a slight overestimation is observed for waterrich mixtures (XMeOH < 0.2). The dielectric constant showed close agreement with experiments, also yielding a correct slope of the predicted curve. Interestingly, the behavior of ΔVmix and ΔHmix, for this pair of models, is poorly described by the geometric rule, and in the case of the excess enthalpy, an endothermic mixing is predicted, contrary to what is observed in experiment. After implementing the modified rule, the excess volume curve was appreciably closer to experiments, and the negative trend of the excess enthalpy was recovered. However, quantitative agreement for this property was still unsatisfactory. An identical analysis applied to the mixtures described by TIP4P-FB and OPLS/2016 yielded similar results (Supporting Information) with λ1 = 0.980 and λ2 = 0.993. A comparative analysis was made bringing together water− methanol simulations19,45 which use the geometric combination rule TIP4P-Ew − OPLS-UA and MD2 − TIP4P/ε and those with the modified combination rule of eq 1 TIP4P-FB − MeOH-4P, TIP4Pε − MeOH-4P, and TIP4P-FB − OPLS/ 2016. We found that ρ, as shown in Figure 7, is best described by MeOH-4P, OPLS-UA, and OPLS/2016. They share slightly overestimated values for low methanol concentrations (XMeOH

MeOH-4P - American Chemical Society

Recommend Documents