Force Field Benchmark of the TraPPE UA for Polar Liquids: Density

Subscriber access provided by CORNELL UNIVERSITY LIBRARY

Article

Force Field Benchmark of the TraPPE UA for Polar Liquids: Density, Heat of Vaporization, Dielectric Constant, Surface Tension, Volumetric Expansion Coefficient and Isothermal Compressibility Edgar Núñez-Rojas, Jorge Alberto Aguilar-Pineda, Alexander Pérez de la Luz, Edith Nadir De Jesús González, and Jose Alejandre J. Phys. Chem. B, Just Accepted Manuscript • Publication Date (Web): 10 Jan 2018 Downloaded from http://pubs.acs.org on January 10, 2018

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a free service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are accessible to all readers and citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

The Journal of Physical Chemistry B is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 32 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

Force Field Benchmark of the TraPPE UA for Polar Liquids: Density, Heat of Vaporization, Dielectric Constant, Surface Tension, Volumetric Expansion Coefficient and Isothermal Compressibility December 15, 2017 ˜ ´ Edgar Nu´ nez-Rojas, Jorge Alberto Aguilar-Pineda, Alexander Pérez de la Luz, Edith Nadir de Jesus González and José Alejandrea∗ a ´ Departamento de Qu´ımica, Universidad Autonoma Metropolitana-Iztapalapa. Av. San Rafael Atlixco 186, Col. Vicentina, 09340, México Distrito Federal, México.

∗

Author to whom correspondence should be addressed. Electronic mail: [email protected]

1

ACS Paragon Plus Environment

The JournalAbstract of Physical Chemistry The transferable potential for phase equilibria force field in its united atom version, TraPPE UA, is

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

evaluated for 41 polar liquids that include alcohols, thiols, ethers, sulfides, aldehydes, ketones and esters to determine its ability to reproduce experimental properties that were not included in the parametrization procedure. The intermolecular force field parameters for pure components were fit to reproduce experimental boiling temperature, vapor-liquid coexisting densities and critical point (temperature, density and pressure) using Monte Carlo simulations in different ensembles. The properties calculated in this work are liquid density, heat of vaporization, dielectric constant, surface tension, volumetric expansion coefficient and isothermal compressibility. Molecular dynamics simulations were performed in the gas and liquid phases, and also at the liquid-vapor interface. We found that relative error between calculated and experimental data is 1.2% for density, 6% for heat of vaporization and 6.2% for surface tension, in good agreement with experimental data. The dielectric constant is systematically underestimated, the relative error is 37%. Evaluating the performance of the force field to reproduce the volumetric expansion coefficient and isothermal compressibility requires more experimental data.

2


Page 2 of 32

Page 3 of 32 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

1

Introduction


Evaluating force fields give insights to improve future parameterizations. In 2011, Vega and Abascal1 evaluated the performance of several non-polarizable models of water over 17 properties in the vapor, liquid and solid phases. The main conclusions were that all the models fail to reproduce, at the same time, the experimental value of melting point, vapor pressure and dielectric constant. The TIP4P/2005 model was the best ranked. The TIP4Q2 and TIP4P/3 non-polarizable models were parametrized in 2011 and 2014, respectively, using as target properties the dielectric constant and temperature of maximum density. The models also reproduce several of the density anomalies of water. The polarizable iAMOEBA force field of water4 was parametrized in 2013 and improved the results of the nonpolarizable models. Different approaches have been used to increase the computational time and sampling properly the phase and configurational spaces. The first approach is to include groups of atoms in a single sphere or bead to develop coarse grain models. The dissipative particle dynamics is one of such methods where it is possible to simulate mesoscopic systems for several microseconds using short ranged soft potential and electrostatics.5 The second is to eliminate the vibrational motion between carbon and hydrogen atoms where the methyl and methylene groups are treated as united atoms, as it is done in the TraPPE UA force field, which decreases the number of atoms involved in the simulations. The third is to include all the atoms where the interaction parameters are given for every atom and the charge distribution can be obtained from quantum mechanical calculations. The computational cost in this latter case is several times larger than the united atom and coarse grain approaches. In 2012, Caleman et al.6 performed a benchmark for OPLS/AA,7 GAFF8 and CGenFF9 force fields to evaluate their ability to reproduce experimental results of density, heat of vaporization, heat capacity, surface tension, isothermal compressibility, volumetric expansion coefficient and dielectric constant for 146 organic liquids. The parametrization procedure of those force fields involves ab initio calculations to obtain the intramolecular parameters and the partial charges of every atom on isolated molecules. The Lennard-Jones (LJ) parameters were obtained with Monte Carlo (MC) or molecular dynamics (MD) simulations to reproduce the liquid density and heat of vaporization7–9 at room conditions. The parameters were obtained for selected atoms on small molecules of the main functional groups of oganic chemistry and then they were made transferable to build other molecules and to study other thermodynamic states. Those force fields are widely used to study biological systems. Caleman et al. found that the calculated dielectric constant was almost half of the experimental value for most of the liquids and that the surface tension was systematically smaller. The other properties were in reasonable good agreement with experimental data. The surface tension was calculated using a slab of liquid surrounded by a vapor with a cut-off distance of 1.1 nm on the LJ potential. In that set up, the surface tension is cut-off dependent.10 Zubillaga et al.11 calculated the surface tension of 61 liquids from the Caleman’s list using the OPLS/AA model with a truncation distance of 2.3 nm. They found that surface tension increased and the agreement with experimental data was improved. In 2015, Fisher et al.12 perfomed simulations for all the organic liquids reported by Caleman using the particle mesh Ewald method on the LJ and Coulomb potentials. With that method is not needed to add long-range corrections to energy and pressure, the calculated bulk and interfacial properties are not cut-off dependent.13 Fisher et al. found that the surface tension increased in around 10 mN/m for 3


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Theother Journal of Physical Chemistry the three force fields. The results for the properties were almost unchanged when the Ewald sum was applied in both potentials. Those works show the importance of calculating properly the molecular interactions in order to parametrize force fields. With the development of the Gibbs ensemble Monte Carlo (GEMC) method14 it was possible to study efficiently the vapor-liquid phase equilibria in a single simulation using a small number of molecules in homogenenous systems where the long-range corrections on the LJ potential are applied to internal energy and pressure. In that way the properties are independent of the truncation distance. In addition, with the introduction of the configurational-bias Monte Carlo method,15, 16 the TraPPE UA17 and NERD18 united atom force fields on hydrocarbons were developed in 1998 to determine orthobaric densities, over a wide range of temperatures and the critical point, of small and long alkanes. The surface tension of linear hydrocarbons containing up to 100 carbons in their molecular structure was obtained using molecular dynamics in the canonical ensemble for the NERD and TraPPE UA force fields.19 The TraPPE UA gives the lowest deviation respect the experimental data for all hydrocarbons. It was found that the OPLS united atom model failed to reproduce the experimental behavior as a function of temperature and pressure.17, 18, 20 The group of Professor Siepmann have developed several transferable potentials for phase equilibria21 that include united atoms (TraPPE UA), explicit hydrogens (TraPPE EH), polarizability (TrAPPE polarizable) and TraPPE-coarse grain. All those force fields were parametrized using MC methods in different ensembles such as GEMC, isothermal-isobaric (NPT) and grand canonical (GCMC) combined mainly with the configurational bias and histogram reweighting techniques. Molecular parameters for all those force field are available in a web page21 and, in particular for the TraPPE UA model of polar liquids, in the original sources.22–28 The idea behind the development of TraPPE force field it was not only to predict several thermophysical properties accurately but also to minimize the number of atom types to build new molecules.28 The intramolecular parameters (bond distances, bending and torsional angles) for the TraPPE UA were taken from other force fields such as OPLS. For polar liquids, the partial charges for linear alcohols, ethers, glycols, ketones and aldehydes were taken from OPLS UA model.29, 30 For the amines, amides, nitriles, pyridine and pyrimidine the charges were taken from previous TraPPE parametrization and the carbon and hydrogen atoms bonded to a nitrogen atom were fixed by considering the molecule was electroneutral. In the case of aromatic compounds they added three additional charges to explicitly represent the out of plane quadrupole moment. The charges for thiols, sulfides, disulfides and thiophene were calculated using the CHELPG (charges from electrostatic potentials using a grid based method) analysis of electrostatic potential energy surfaces derived from ab initio calculation at the HF/6-31g+(d,p) level on isolated molecules. The Lennard-Jones parameters for the TraPPE UA model were obtained by matching the experimental vapor-liquid densities at different temperatures along the coexistence line, the critical point (density, pressure and temperature) and the boiling temperature. The parameters were also used to simulate some binary mixtures in the vapor-liquid equilibrium. The TraPPE UA force field improves several thermodynamic properties of the OPLS, GAFF and CHARMM force fields as a function of temperature and pressure. Developing efficient force fields requires experimental data for several properties. The TraPPE UA model uses coexisting densities at the vapor-liquid equilibrium and critical properties, however, there is a lack of information for several polar liquids such as thiols26 and acrylates.28 Work in that direction is desirable to develop force field with a better performance. In this work we performed molecular dynamics simulations on 41 organic liquids to calculate the

4


Page 4 of 32

Page 5 of 32 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journalconstant, of Physical Chemistry liquid density, heat of vaporization, The dielectric surface tension, volumetric expansion coefficient and isothermal compressibility. The systems include alcohols, thiols, ethers, sulfides, aldehydes, ketones and esters. The main goal is to evaluate the abbility of TraPPE UA force field for polar liquids to reproduce experimental data that were not used as target properties to fix the interaction parameters. The results will also provide useful information to understand the limits of the parametrization procedure applied on non-polarizable force fields of polar liquids. The manuscript is organized as follows: In Section 2 the force field is defined. The simulation details are given in Section 3. The MD results are discussed in Section 4. Finally, the main conclusions are hightlighted together with a list of references. The equivalence between results from OPLS and Ryckaert-Bellemans (RB) functions in the dihedral angle is discussed in Appendix A. The analysis of choosing optimum thermostat and barostat parameters for a given time step is given in Appendix B.

2

TraPPE UA force field

The functional form of the non-polarizable force field contains bonded and non-bonded interactions where the bond distances are considered to be rigid. The bending of an angle, θ, is described by a harmonic potential, kθ (θ − θ0 )2 (1) 2 where kθ is the spring constant and θ0 is the equilibrium angle. The torsional potential of angle φ is of the OPLS type with two functional forms. The first type is, Ubend (θ) =

Utors1 (φ) = F0 + F1 [1 + cos(φ)] + F2 [1 − cos(2φ)] + F3 [1 + cos(3φ)] where Fi are the Fourier coefficients. The convention φtrans =180 is assumed. The second type is, 6 X Gi cos(iφ) Vtors2 (φijkl ) = G0 +

(2)

(3)

i=1

The Gromacs molecular dynamics package has implemented the function given by equation (2) with parameters Fig = 2Fi and F0g =0. The variable Fig denotes the Fourier coefficients used in the topology file within a Gromacs simulation, see the manual.31 The second type of dihedral potential defined by equation (3) is not implemented in Gromacs. In this work, the functional form of the torsional potential given by equations (2) and (3) were transformed to a Ryckaert-Bellemans32 function to perform all the MD simulations, URB (ψ) =

5 X

Cn [cos(ψ)]n =

n=0

5 X

(−1)n Cn [cos(φ)]n

(4)

n=0

where ψ = φ − 180 and ψtrans =0. The RB function is implemented in the Gromacs program. The relations between Fi and Gi with the Cn values are given in Appendix A. In order to check consistency between the results from this work with those reported for TraPPE UA using MC methods, we performed NPT MD simulations on methyl tert-butyl ether at 298.15 K and 1

5


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal ofangles Physical Chemistry bar, using the Gromacs package, with The the dihedral described by equations (2) and (4). The results of angle distributions and several physical properties are given in Appendix A. The non-bonded interactions, atoms in different molecules or within the same molecule but separated by four or more bonds, are described by a pairwise-additive potential that contains Lennard-Jones and Coulombic contributions, " V (rij ) = 4εij

σij rij

12

−

σij rij

6 # +

qi qj 4πε0 rij

(5)

where rij is the distance between atoms i and j, σij and εij are the size and energy LJ parameters, respectively. The permitivity of vacuum is ε0 and the charge of atom i is denoted by qi . The intramolecular Coulombic interactions separated by three bonds in esters28 groups were reduced by a factor of 1/2, as in the original works. The Lorentz-Berthelot mixing rules are used to calculate the interactions between atoms of different type, σij = (σii + σjj )/2 and ij = (ii ij )1/2 . We do not simulate systems that include a repulsive term to calculate intramolecular interactions between atoms separated for four bonds as in glycol ethers24 and acrylates.28 Systems with benzene rings, where additional charges out of the plane were included in the parametrization procedure, are also not simulated. The TraPPE UA parameters used to simulate polar liquids in this work are obtained from a web page21 and from the original sources: Alcohols,23 ethers, glycols, ketones and aldehydes,24 thioles, sulfides, disulfides26 and esters.27 We note that for 2-ethylhexyl acrylate the LJ parameters for the united atom CH, denoted as (CHx)2-[CH]-CHx, are reported in the web page but they are not given in Table 3 of Ref.28 The values used to perform the MD simulations are σ=0.468 nm and =0.083 kJ/mol. On the other hand, the parameter C1 of the dihedral angle CH3-O-C=O in Table 3 of Ref.27 has a wrong sign, the correct value is -2194 K (private communication), which is used to perform the simulations of acetates. The coordinates, force field parameters and molecular geometry for all the systems studied in this work are provided as Supporting Information in the format of the Gromacs package.

3 Simulation details We made a careful analysis to choose appropriately the time step, thermostat and barostat parameters involved in the solution of equations of motion in molecular dynamics using the Gromacs package.33–36 This is important to warrant the results are independent of such parameters and to discard that differences between MC methods, used to parametrize the TraPPE UA force field, and MD do not come from reasons related with the methodology. The stability of Berendsen and Nosé-Hoover/ParrinelloRahman thermostats and barostats37–40 is evaluated for ethyl acetate where we have observed that the liquid density is sensitive to the simulation parameters. The analysis made for the parameters involved in the solution of the equations of motion shows that Parrinello-Rahman barostat is better than the Berendesen barostat to produce correct volume fluctuations. The results show that for rigid bonds, the time step of 0.002 ps is adequate if τT , the thermostat parameter, and τP , the barostat parameter, are 0.5 ps and 1.0 ps, respectively, see Appendix B. The results are the same, within error in the simulation, if τP > 1. When the values of τT = τP =0.5, the 6


Page 6 of 32

Page 7 of 32 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Thefor Journal of Physical density and isothermal compressibility the ethyl acetateChemistry are 1.25 and 5.6 times larger, respectively, compared with those obtained with τP = 1.0 ps. A bad choice of thermostat and barostat parameters can give wrong results of physical properties affecting the parametrization of force fields. The volume fluctuations are used to calculate the compressibility factor, it is therefore important to evaluate their proper sampling during the simulations. In addition, the heat of vaporization, dielectric constant and compressibility factor are also obtained. The simulations were perfomed using the version 4.5.4 of Gromacs where the equations of motion were solved using the leap-frog algorithm with a time step of 2 fs. Periodic boundary conditions were used in all directions. The bond distances were kept rigid by using the LINCS algorithm.41, 42 The electrostatic interactions were treated with the Ewald sums with a tolerance of 1 × 10−6 for the real space contribution. The real part of the Ewald summation and the LJ interactions were truncated at 1.4 nm including long range corrections to the energy and pressure in the case of homogeneous systems. The long range contribution of electrostatic interactions was determined through the particle mesh Ewald method43 with a grid of reciprocal vectors of 1.2 nm and spline of order 4.

3.1 NPT liquid simulations They were carried out in the isothermal-isobaric ensemble (constant number of molecules, pressure and temperature). NPT MD simulations at 298.15 K and 1 bar are performed on methyl tert-butyl ether to analyze the effect that number of molecules has on several properties, see Appendix B. The number of molecules used in this work was 1000. The internal temperature was coupled to a Nosé-Hoover thermostat with a time parameter τT =0.5 ps while the pressure was coupled to a Parrinello-Rahman barostat with a time parameter τP =1.0 ps. The calculated properties from these simulations were the density, dielectric constant, isothermal compressibility and volumetric expansion coefficient. The potential energy was also obtained and it was used to calculate the heat of vaporization. The average properties were obtained for 50 ns after an equilibration period of 20 ns.

3.2 NVT vapor-liquid simulations They were carried out in the canonical ensemble (constant number of molecules, volume and temperature). The number of molecules was 4000. The internal temperature was coupled to a Nosé-Hoover thermostat with a time parameter τT =0.5 ps. The calculated property from these simulations was the surface tension. The initial set up was a slab of a liquid surrounded by two empty regions in a parallelepiped cell. The system evolved to reach equilibrium. The final configuration contained liquid and vapor phases in contact. The dimensions of the simulation cell were at least Lx = Ly =5.1 nm, large enough to avoid finite size effects of surface area on the surface tension.44 The average value of Lz was 45 nm. The smallest and largest values were 30 nm for acetaldehyde and 85 nm for octanethiol, respectively. The size of the liquid region in all the systems varied from 15 nm to 30 nm. The truncation distance was 2.5 nm to calculate properly the long range interactions on the LJ potential.11 The average components of the pressure tensor were obtained for 20 ns after an equilibration period of 5 ns.

7


3.3 NVT gas simulations 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Page 8 of 32

They were performed using a stochastic dynamics integrator6 in the canonical ensemble as implemented in the Gromacs software. The temperature was the same as that used to calculate the liquid density and the time parameter was 2.0 ps. The system contains one molecule and the simulations were run for 2 ns after an equilibation period of additional 1 ns. The potential energy was obtained and it was used to calculate the heat of vaporization.

4

Results

The differences between MD results from this work and those reported using MC methods can not be associated to methodological aspects neither to the use of the Ryckaert-Bellemans function for the dihedral angle. We have checked that the distribution of dihedral angles are the same as that obtained with the functional forms used in the parameterization of TraPPE UA models, see Appendix B. The differences in liquid density and heat of vaporization come from the use of different number of molecules and the pressure where the density was calculated, we use 1 bar in NPT simulations and in the MC simulations is used the vapor pressure at the vapor-liquid equilibrium. The MC pressures are close to 1 bar but in general they are slightly different. In order to evaluate the force field, the percentage of the absolute relative error on the calculated properties, XM D , respect to experimental values, Xexp , is obtained as, ∆X = |(Xexp − XM D )| ∗ 100/Xexp

(6)

The relative error was also calculated between the values of liquid density and heat of vaporization from MD and those reported with MC. In that case Xexp is replaced by XM C in equation (6). We made an analysis of the systems that have a larger value of the relative error in every property. The results are given in bold face letter type in the Tables S1-S6 as Supporting Information.

4.1 Liquid density The liquid density was obtained using, M where M is the mass of the system and < V > the average volumen of the simulation cell. ρ=

(7)

The MD calculations are performed at temperatures where experimental liquid densities are available at 1 bar. The results for 41 polar liquids from NPT MD simulations are shown in Fig. 1. Only 32 MC results reported in the web page21 were used in the comparison with MD calculations because the temperature was the same or because it was possible to estimate the density from a linear fit using the results of two temperatures close to the experimental data. We recall that the density reported in MC calculations corresponds to the liquid phase in the vapor-liquid equilibrium. In that conditions, the MC pressure is slightly different from 1 bar. The average relative error between MD and MC for 32 results is 1.0%. The same quantity respect to experimental data is 1.2% and 1.6% for MD and MC results, respectively. The agreement between MD and MC results is excellent. The effect of using different number

8


1100

1000

MC [21]

3

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry of molecules, pressure and dihedral functions seems to be small for the molecules studied in this work. The numerical values are given in Table S1 as Supporting Information. The far right data point corresponds to dimethyl disulfide. It is the system with the largest density. We checked carefully the input force field parameters and simulations details. The simulation was run for longer time than the other systems. The MC calculations were performed in the grand canonical ensemble. The insertion of molecules is difficult at high densities. That might be one of the possible sources of error, even though in the original work were used special techniques which sample properly the configurational space.

ρCalc(kg/m )

Page 9 of 32

900

800

700 700

800

900 3 ρExp(kg/m )

1000

1100

Figure 1: Correlation between calculated and experimental?, 6, 45, 46 liquid density for the TraPPE UA force field at temperatures given in Table S1 (Supporting Information) and 1 bar. The 41 MD results from this work are shown with green circles. The 32 results reported in the web page21 are shown with green squares. The maximum estimated error in the calculated density is 0.4 kg/m3 and it is less than the symbol size. The symbols are as follows: squares for alcohols, circles for thiols, diamonds for ethers, triangles for sulfides, pluses for aldehydes, equis for ketones and stars for esters.

4.2 Heat of vaporization The heat of vaporization, the enthalpy change produced in the transition of one mol of liquid to the vapor phase at equilibrium conditions, is calculated assuming that the gas is ideal. We do not correct the results to take into account the quantum vibrational and polarization effects,47 ∆Hvap = Hg − Hl = RT + Ug − (U + P V )l

(8)

where H is the enthalpi, U is the potential energy and the subscripts g and l denote gas and liquid phases, respectively.

9


90 80

MC [23,28]

70 ∆Hvap,Calc(kJ/mol)

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Physical Chemistry The calculations were performed The at temperatures where experimental data are available at 1 bar, for that reason, the temperature in some systems is not the same as that used to calculate density. The potential energy of the gas phase was obtained for one molecule. The heat of vaporization results are shown in Fig. 2 and the numerical values are given in Table S2 as Supporting Information. Although ∆Hvap is not a target property in the parametrization of the TraPPE UA force field, the agreement with experimental data is good. The relative error between 39 results from this work and experimental data is 6.0%. The ∆Hvap from MC were obtained using the potential energy, saturated pressure and volume in GEMC simulations when the vapor pressure is not too small, otherwise, it is calculated using the NPT method assuming the gas is ideal.23, 28 The relative error between 12 points from MD and MC is 3.2%.

60 50 40 30 20 10 10

20

30

40 60 50 ∆Hvap,Exp(kJ/mol)

70

80

90

Figure 2: Correlation between calculated and experimental heat of vaporization for polar liquids obtained from NPT MD simulations using the TraPPE UA force field. The experimental values were taken from references.6, 28, 30, 45, 46, 48 The MC results for alcohols and acrylates were reported at 300 K23 and 298 K,28 respectively, at the vapor-liquid equilibrium at almost zero pressure. The maximum estimated error in the MD calculated value is 0.2 kJ/mol and it is smaller than the size of the symbol. The symbols have the same meaning as in Fig. 1.

4.3 Dielectric constant The static dielectric constant, , is obtained using the dipole moment, M, fluctuations of the system,49 4π (< M2 > − < M >2 ) (9) 3kB T V where kB is the Boltzmann constant, T the absolute temperature and V the simulation cell volume. =1+

The calculations were performed at temperatures where experimental data are available at 1 bar. Most of the results are obtained from the liquid simulations to obtain the density. The results are shown 10


Page 10 of 32

εCalc

Page 11 of 32 in Figure 3 and the numerical values given The Journal of S3 Physical ChemistryInformation. The calculated values in Table as Supporting are systematically lower than experimental data. The average relative error is 37%. We checked the systems were well equilibrated, the average value of the dipole moment of the system was around cero. 1 Caleman et al.6 also found that the dielectric constant for OPLS AA and GAFF force fields was under2 3 estimated for most of the systems and they argued that deviations from experimental data was due to 4 a lack of polarization of the molecules. Several atomic charges of polar liquids used in the parametriza5 tion procedure of TraPPE UA were taken from OPLS and AMBER force fields23, 24, 29, 30 or calculated 6 using the CHELPG (charges from electrostatic potentials using a grid based method) analysis of elec7 8 trostatic potential energy surfaces from ab initio calculation at the HF/6-31g+(d,p) level of theory on 9 isolated molecules.24, 26 The reason why the OPLS/AA, GAFF and TraPPE UA models fail to reproduce 10 the dielectric constant might be related with the use of atomic charges from quantum mechanical cal11 culations on isolated molecules. In a recent work we used the dielectric constant as a target property 12 13 and the charges of different force fields were scaled linearly to reproduce its experimental value.50 The 14 procedure has been used to study acetamide and formamide where the dielectric constants are 60 and 15 110, respectively. We found that is important to reproduce the correct solubility of hexan-2-one in for16 mamide51 and acetamide in water.52 17 18 19 30 20 21 22 25 23 24 25 26 20 27 28 29 15 30 31 32 10 33 34 35 5 36 37 38 0 39 0 10 20 30 5 15 25 40 εExp 41 42 Figure 3: Correlation between calculated and experimental6, 45, 53–59 dielectric constant for the 43 TraPPE UA force field at 1 bar and temperatures were experimental values are available. The aver44 age relative error for 34 values is 37.0%. The maximum estimated error in the calculated value is 0.4 and 45 it is around the same size of the symbol. The symbols have the same meaning as in Fig. 1. 46 47 48 49 4.4 Surface tension 50 51 The surface tension of a planar interface using the mechanical definition of the atomic pressure60, 61 is, 52 53 54 11 55 56 57 58 59 ACS Paragon Plus Environment 60


Page 12 of 32

γ = 0.5Lz [< Pzz > −0.5(< Pxx > + < Pyy >)]

where Lz is the length of the simulation cell in the longest direction and Pαα are the diagonal components of the pressure tensor. The factor 0.5 outside the squared brackets takes into account the two symmetrical interfaces in the system. The results are shown in Fig. 4 and the numerical values given in Table S4 as Supporting Information. The calculations were performed at temperatures where experimental data are available. We checked that the average pressure components Pxx and Pyy in the interface simulations were equal, within simulation error, and that the density of the liquid phase was also equal to that obtained from NPT MD at the same temperature. The agreement with experimental data is good, the average relative error between calculated and experimental values is 6.2%. It is seems that a force field that reproduces the critical point has a good chance to describe quantitatively also the surface tension.

50

40 γCalc(mN/m)

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

(10)

30

20 20

30 γExp(mN/m)

40

50

Figure 4: Correlation between calculated and experimental6, 45, 46, 62, 63 surface tension at the vapor-liquid interface for the TraPPE UA force field. The relative error for 36 data is 6.2%. The symbols have the same meaning as in Fig. 1.

4.5 Volumetric expansion coefficient The thermal expansion coefficient, αP , was obtained using, 1 αP = − ρ

∂ρ ∂T

(11) P

Liquid simulations in the NPT ensemble were performed from 290 K to 310 K in steps of 5 K at 1 bar using the annealing procedure as implemented in the Gromacs package. The calculations are performed at temperatures were experimental data are available. The density was fitted to a linear equation of density as a function of temperature which is used to obtain αP . The results are shown in Fig. 5 and 12


αP,Calc(1/K)

Page 13 of 32 the numerical data are given in TableThe of Physical Chemistry The agreement with experiment S5 Journal as Supporting Information. is good. The average relative error between 19 calculated and experimental data is 7.6%. The choice of coexisting densities along the vapor-liquid equilibrium as target properties to parametrize a force field 1 seems to give also good results for the volumetric expansion coefficient. More experimental information 2 3 is needed to have a better evaluation of this property. 4 5 6 2.5 7 8 9 2 10 11 12 13 1.5 14 15 16 1 17 18 19 20 0.5 21 22 23 0 24 0 1 2 0.5 1.5 2.5 25 αP,Exp(1/K) 26 27 Figure 5: Correlation between calculated and experimental6, 45, 46 volumetric expansion coefficient. The 28 symbols have the same meaning as in Fig. 1. 29 30 31 32 4.6 Isothermal compressibility 33 34 The isothermal compressibility, κT , was obtained from volume fluctuations, 35 36 < V 2 > − < V >2 37 κT = (12) kB T < V > 38 The calculations are performed at temperatures were experimental data are available. The results 39 40 are shown in Fig. 6 and the numerical data are given in Table S6 as Supporting Information. 41 42 We made in Appendix B an analysis to determine the optimum time step, thermostat and barostat 43 parameters involved in the solution of the equations of motion. We are confident that the results of 44 45 κT are well calculated with the use of Nosé-Hoover thermostat and Parrinello-Rahman barostat. More 46 experimental information is needed to evaluate the performance of TraPPE UA model to reproduce this 47 property. 48 49 50 51 52 53 54 55 56 57 58 59 60

13


2.5

Page 14 of 32

2 κT,Calc(1/GPa)

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


1.5 1 0.5 0 0

0.5

1 1.5 κT,Exp(1/GPa)

2

2.5

Figure 6: Correlation between calculated and experimental6, 45, 64–66 isothermal compressibility. The average relative error for 12 results is 13.8%. The symbols have the same meaning as in Fig. 1.

5

Conclusions

Excellent agreement for liquid densities between 32 results obtained in this work using the MD method with those reported for the TraPPE UA from MC simulations. The experimental heat of vaporization and surface tension, which are not target properties in the TraPPE UA parametrization procedure, are also well reproduced. The calculated dielectric contant is systematically lower than experimental values for most of the simulated liquids. The relative error for 34 data is 37%. The same trend was observed by Caleman et al.6 for most of the polar liquids studied for them using different force fields. The reason of that behaviour in both cases might be related with the use of charges obtained from quantum mechanical calculations of isolated molecules where the electric field around a molecule is not taken into account. This is a new physical insight which indicates that better charge distributions are needed to simulate polar liquids. The TraPPE UA force field has to be used with care because it can fail to reproduce the experimental miscibility of one liquid in another when there is a large difference in the values of the dielectric constants.51 The relative error between calculated and experimental results is 7.6% for αP and 14.1% for κT . More experimental data are needed to evaluate the perfomance of the TraPPE UA force field to reproduce both properties but also to develop new parameters and force fields. The analysis of outliers, systems with relative error larger than the average value for every property and highlighted with bold face letter type in Tables of Supporting Information, shows that MD results for ethers, sulfides and esters have a large deviation respect to experimental values for at least 4 properties 14


Page 15 of 32 (density, heat of vaporization, dielectric The Journal and of Physical constant surfaceChemistry tension). It is strange the failure of the force field for ethers because there is experimental data for all the critical properties that would have helped to develop good parameters. For sulfides and esters there is a lack of experimental coexisting densities and 1 critical data. That might be the reason why the parameterization fail to predict the properties obtained 2 3 in this work. 4 The results from this work show that the procedure applied to parametrize the TraPPE UA force 5 field is better than that used for the popular force fields OPLS/AA, GAFF and CGenFF. 6 7 8 Supporting information 9 10 The geometries of all the systems studied are supplied. 11 Tables with experimental and calculated data. 12 13 Topologies of molecular systems studied. 14 15 Acknowledgements 16 17 We would like to thank Conacyt for financial support under the project 257422. APL, JAP and 18 19 NJG thank Conacyt for a graduate scholarship. The support and computer time of Laboratorio de Su20 ´ ´ en Paralelo de la UAM-Iztapalapa is greatly appreciated. percomputo y Visualizacion 21 22 23 24 References 25 26 [1] Vega, C.; Abascal, L. F. Simulating water with rigid non-polarizable models: a general perspective. 27 Phys. Chem. Chem. Phys. 2011, 13, 19663-19688. 28 29 [2] Alejandre, J.; Chapela, G. A.; Saint-Martin, H.; Mendoza, N. Non-polarizable model of water that 30 yields the dielectric constant and the density anomalies of the liquid: TIP4Q. Phys. Chem. Chem. 31 Phys. 2011, 13, 19728-19741. 32 33 [3] Fuentes-Azcatl, R.; Alejandre, J. Non-polarizable force field of water based on the dielectric con34 35 stant: TIP4P/. J. Phys. Chem. B. 2014, 118, 1263-1272. 36 [4] Wang, L-P.; Head-Gordon, T.; Ponder, J. W.; Ren, P.; Chodera, J. D.; Eastman, P. K.; Mart´ınez, T. J.; 37 38 Pande, V. S. Systematic improvement of a classical molecular model of water. J. Phys. Chem. B. 2013, 39 117, 9956-9972. 40 41 [5] González-Melchor, M.; Mayoral E.; Velázquez, M. E.; Alejandre J. Electrostatic interactions in dissi42 pative particle dynamics using the Ewald sums. J. Chem. Phys. 2006, 125, 224107-1 - 224107-8. 43 44 [6] Caleman, C.; van Maaren, P. J.; Hong, M.; Hub, H. S.;Costa, L. T.; van der Spoel, D. Force Field 45 Benchmark of Organic Liquids: Density, Enthalpy of Vaporization, Heat Capacities, Surface Ten46 47 sion, Isothermal Compressibility, Volumetric Expansion Coefficient, and Dielectric Constant. J. 48 Chem. Theory Comput. 2012, 8, 61-74. 49 50 51 52 53 54 15 55 56 57 58 59 ACS Paragon Plus Environment 60

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Physical Chemistryand Testing of the OPLS All-Atom [7] Jorgensen, W. L.; Maxwell, D. S.;The Tirado-Rives, J. Development Force Field on Conformational Energetics and Properties of Organic Liquids. J. Am. Chem. Soc. 1996, 118, 11225-11236. [8] Wang, J.; Wolf, R. M.; Caldwell, J. W.; Kollman, P. A.; Case, D. A. Development and testing of a general amber force field. J. Comput. Chem. 2004, 25, 1157-1174. [9] Vanommeslaeghe, K.; et al. CHARMM General Force Field (CGenFF): A force field for drug-like molecules compatible with the CHARMM all-atom additive biological force fields. J. Comput. Chem. 2010, 31, 671-690. [10] Truckymchuk, A.; Alejandre, J. Computer simulations of liquid/vapor interface in Lennard-Jones fluids: Some questions and answers. J. Chem. Phys. 1999, 111, 8510-8523. [11] Zubillaga, R. A.; Labastida, A.; Cruz, B.; Mart´ınez, J. C.; Sánchez, E.; Alejandre, J. Surface Tension of Organic Liquids Using the OPLS/AA Force Field. J. Chem. Theory Comput. 2013, 9, 1611-1615. [12] Fisher, N, M.; van Maaren, P. J.; Ditz, J. C, Yildirim, A.; van der Spoel, D. Properties of organic liquids when simulated with Long-Range Lennard-Jones interactions. J. Chem. Theory Comput. 2015, 11, 2938-2944. ´ [13] Lopez-Lemus, J.; Alejandre, J. Thermodynamic and transport properties of simple fluids using lattice sums: bulk phases and liquid-vapour interface. Mol. Phys. 2002, 100, 2883-2992. [14] Panagiotopoulos, A. Direct determination of phase coexistence properties of fluids by Monte Carlo simulation in a new ensemble. Molec. Phys. 1987, 61, 813-826. [15] Siepmann, J. I. A method for the direct calculation of chemical potentials for dense chain systems. Mol. Phys. 1990, 70, 1145-1158. [16] de Pablo, J. J.; Laso, M.; Suter, U. W. Simulation of polyethylene above and below the melting point. The Journal of chemical physics. 1992, 96, 2395-2403. [17] Martin, M. G.; Siepmann, J. I. Transferable Potentials for Phase Equilibria. 1. United-Atom Description of n-Alkanes. J. Phys. Chem. B. 1998, 102, 2569-2577. [18] Nath, S. K.; Escobedo, F. A.; de Pablo, J. J. On the simulation of vapor-liquid equilibria for alkanes. The Journal of Chemical Physics. 1998, 108, 9905-9911 [19] Muller, E. A.; Mej´ıa, A. Comparison of United-Atom potentials for the simulation of vapor-liquid equilibria and interfacial properties of long-chain n-alkanes up to n-C100. J. Phys. Chem. B. 2011, 115, 12822-12834. [20] van Leeuwen, M. E.; Smit, B. Molecular simulation of the vapor-liquid curve of methanol. J. Phys. Chem. 1995, 99, 1831. [21] The Siepmann Group http://chem-siepmann.oit.umn.edu/siepmann/trappe/index.html (last accessed Aug 2017) [22] Wick, C. D.; Martin, M. G.; Siepmann, J. I. Transferable potentials for phase equilibria. 4. Unitedatom description of linear and branched alkenes and alkylbenzenes. J. Phys. Chem. 2000, 104, 80088016. 16


Page 16 of 32

Page 17 of 32 [23] Chen, B.; Potoff, J. J.; Siepmann, The JournalCarlo of Physical Chemistry J. I. Monte calculations for alcohols and their mixtures with alkanes. Transferable potentials for phase equilibria. 5. United-atom description of primary, secondary and tertiary alcohols. J. Phys. Chem. B. 2001, 105, 3093-3104. 1 2 [24] Stubbs, M.; Potoff, J. J.; Siepmann, J. I. Transferable potentials for phase equilibria. 6. United-atom 3 description for ethers, glycols, ketones and aldehydes. J. Phys. Chem. B. 2004, 108, 17596-17605. 4 5 [25] Wick, C. D.; Stubbs, J. M.; Rai, N.; Siepmann, J. I. Transferable potentials for phase equilibria. 7. 6 United-atom description for nitrogen, amines, amides, nitriles, pyridine and pyrimidine. J. Phys. 7 8 Chem. B. 2005, 109, 18974-18982. 9 10 [26] Lubna, N.; Kamath, G.; Potoff, J. J.; Rai, N.; Siepmann, J. I. Transferable potentials for phase equilib11 ria. 8. United-atom description for thiols, sulfides, disulfides, and thiophene. J. Phys. Chem. B. 2005, 12 109, 24100-24107. 13 14 [27] Kamath, G.; Robinson, J.; Potoff, J. J. Application of TraPPE-UA force field for determination of 15 vapor-liquid equilibria of carboxylate esters. Fluid Phase Equilibria. 2006, 240, 46-55. 16 17 [28] Maerzke, K. A.; Schultz, N. E.; Ross, R. B.; Siepmann J. I. TraPPE-UA force field for acrylates and 18 Monte Carlo simulations for their mixtures with alkanes and alcohols. J. Phys. Chem. B. 2009, 113, 19 20 6415-6425. 21 [29] Jorgensen, W. L. Optimized Intermolecular Potential Functions For Liquid Alcohols. J. Phys. Chem. 22 23 1986, 90, 1276-1284. 24 25 [30] Briggs, J. M.; Matsui, T.; Jorgensen, W. L. Monte Carlo simulations of liquid alkyl ethers with the 26 OPLS potential functions. J. Comp. Chem. 1990, 11, 958-971. 27 28 [31] van der Spoel, D.; et al. Gromacs User Manual version 4.5.6, www.gromacs.org (2010). 29 30 [32] Ryckaert, J. P.; Bellemans, A. Molecular dynamics of liquid n-butane near its boiling point. Chem. 31 Phys. Lett. 1975, 30, 123-125. 32 33 [33] Berendsen, H. J. C.; van der Spoel, D.; van Drunen, R. GROMACS: A message passing parallel 34 molecular dynamics implementation. Comput. Phys. Commun. 1995, 91, 43-56. 35 36 [34] Lindahl, E.; Hess, B.; van der Spoel, D. J. GROMACS 3.0: A package for molecular simulation and 37 trajectory analysis. Mol. Model. 2001, 7, 306-317. 38 39 [35] van der Spoel, D.; Lindahl, E.; Hess, B.; Groenhof, G.; Mark, A. E.; Berendsen, H. J. C. GROMACS: 40 Fast, flexible, and free. J. Comput. Chem. 2005, 26, 1701-1718. 41 42 [36] Hess, B.; Kutzner, C.; van der Spoel, D.; Lindahl, E. GROMACS 4: Algorithms for highly efficient, 43 44 load-balanced, and scalable molecular simulation. J. Chem. Theory Comput. 2008, 4, 435-447. 45 [37] Nosé, S. A molecular dynamics method for simulations in the canonical ensemble. Mol. Phys. 1984, 46 47 52, 255-268. 48 49 [38] Hoover, W. G. Canonical dynamics: Equilibrium phase-space distributions. Phys. Rev. A. 1985, 31, 50 1965-1697. 51 52 53 54 17 55 56 57 58 59 ACS Paragon Plus Environment 60

The Journal of Physical [39] Parinello, M.; Rahman, A. Polymorphic transitions in Chemistry single crystals: A new molecular dynamics method. J. of Applied Physics. 1981, 52, 7182-7190. 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

[40] Tuckerman, M. E. Statistical mechanics: Theory and molecular simulations. Ed. Oxford University Press, 2010. [41] Hess, B.; Bekker, H.; Berendsen, H. J. C.; Fraaije, J. G. E. M. LINCS: A linear constraint solver for molecular simulations. J. Comput. Chem. 1997, 18, 1463-1472. [42] Hess, B. P-LINCS: A parallel Linear Constraint Solver for Molecular Simulation. J. Chem. Theory Comput. 2008, 4, 116-122. [43] Darden, T.; York, D.; Pedersen, L. Particle mesh Ewald: An N-log(N) method for Ewald sums in large systems. J. Chem. Phys. 1993, 98, 1089-1092. ´ [44] Orea, P.; Lopez-Lemus, J.; Alejandre J. Oscillatory surface tension due to finite-size effects. J. Chem. Phys. 2005, 123, 114702-114813. [45] Lide, D. R.; Haynes, W. M. editors, CRC Handbook of Chemistry and Physics, 90th ed. CRC Press, Boca Raton, FL, 2009. [46] Yaws, C. L. Chemical Propperties Handbook, McGraw-Hill, 1999. [47] Horn, H. W.; Swope, W. C.; Pitera, J. W.; Madura, J. D.; Dick, T. J.; Hura, G. L.; Head-Gordon, T. Development of an improved four-site water model for biomolecular simulations: TIP4P-Ew. J. Chem. Phys. 2004, 120, 9665-9678. [48] National Institute of Standards and Technology http://webbook.nist.gov/ (last accessed Sep 2017). [49] Neumann, M. Dipole moment flutuation formulas in computer simulations of polar systems. Mol. Phys. 1983, 50, 841-858. ˜ [50] Salas, F. J.; Méndez-Maldonado, G. A.; Nu´ nez-Rojas, E.; Aguilar-Pineda, G. E.; Dom´ınguez, H.; Alejandre, J. Systematic procedure to parametrize force fields for molecular fluids. J. Chem. Theory and Comput. 2015, 11, 683-693. ˜ [51] Pérez de la Luz, A.; Méndez-Maldonado, G. A.; Nu´ nez-Rojas, E.; Bresme, F.; Alejandre, J. A. New Force Field of Formamide and the Effect of the Dielectric Constant on Miscibility. J. Chem. Theory Comput. 2015, 11, 2792-2800. ˜ [52] Aguilar-Pineda, J. A.; Méndez-Maldonado, G. A.; Nu´ nez-Rojas, E.; Alejandre, J. Parametrisation of a force field of acetamide for simulations of the liquid phase. Mol. Phys. 2015, 113, 17-18. [53] Osipov, O. A. On the relationship between the dielectric permeability of a polar liquid and the dipole moment. Zh. Fiz. Khim., 1957, 31, 1542-1546. [54] Schneider, R. L. Physical properties of some organic solvents. Eastman Org. Chem. Bull., 1975, 47, 1-12. [55] Cole, K. S.; Cole, R. H. Dispersion and Absorption in Dielectrics: I - Alternating Current Characteristics. J. Chem. Phys., 1941, 9, 341-352.

18


Page 18 of 32

Page 19 of 32 [56] Lomba, L.; Giner, B.; Lafuente, C.; The JournalS.;ofArtigas, Physical H. Chemistry Martin, Thermophysical properties of three compounds from the acrylate family J. Chem. Eng. Data. 2003, 58, 1193. 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

[57] Sastry, N. V.; Patel, S. R.; Soni, S. S. Excess molar volumes, excess isentropic compressibilities, excess viscosities, relative permittivity and molar polarization deviations for methyl acetate +, ethyl acetate +, butyl acetate +, isoamyl acetate +, methyl propionate +, ethyl propionate +, ethyl butyrate +, methyl methacrylate +, ethyl methacrylate +, and butyl methacrylate + cyclohexane at T = 298.15 and 303.15 K J. Mol. Liq. 2013, 183, 102. [58] Rewar, G. D.; Bhatnagar, D. Microwave absorption and relaxation processes of ternary mixtures of non-rigid polar liquids Indian J. Pure Appl. Phys. 2001, 39, 707. [59] George, J.; Sastry, N. V. Thermophysical properties of binary mixtures of methyl methacrylate+diethers (ethyl, isopropyl, and butyl) at 298.15 and 308.15 K. Int. J. Thermophys. 2003, 24, 1697. ´ [60] Lopez-Lemus, J.; Alejandre, J. Thermodynamic and transport properties of simple fluids using lattice sums: bulk phases and liquid-vapour interface. Mol. Phys. 2002, 100, 2983-2992. [61] Brown, D.; Neyertz, S. A general pressure tensor calculation for molecular dynamics simulations. Mol. Phys. 1995, 84, 577-595. [62] Wohlfarth, Ch. Surface Tension of Pure Liquids and Binary Liquid Mixtures. (Supplement to IV/16). Lechner, M. D. (Ed.), 1997, VII, 439. [63] DOW Corporation https://www.dow.com/en-us/products/Acrylates (last accessed Aug 2017). [64] Kartsev, V. N.; Zabelin, V. A. Isothermal Compressibilities of Liquid Straight-chain Alcohols. Russ. J. Phys. Chem. 1978, 52, 1221-1222. [65] Kulikov, D.; Verevkin, S. P.; Heintz, A. Determination of Vapor Pressures and Vaporization Enthalpies of the Aliphatic Branched C5 and C6 Alcohols. J. Chem. Eng. Data. 2001, 46, 1593-1600. [66] Richard, A. J. Isothermal compressibility of methyl ketones and the density coefficient of polarizability for several compounds by ultracentrifugation. J. Phys. Chem. 1978, 82, 1265-1268. ´ ´ R.; Jochim, A. L.; Martyna, G. L. A Liouville[67] Tuckerman, M. E.; Alejandre, J.; Lopez-Rend on, operator derived measure-preserving integrator for molecular dynamics simulations in the isothermal-isobaric ensemble. J. Phys. A: Math. Gen. 2006, 39, 5629-5651. ´ ´ R.; Martyna, G.L.; Tuckerman, M. E. Measure-preserving [68] Yu, T-Q.; Alejandre, J.; Lopez-Rend on, integrators for molecular dynamics in the isothermal-isobaric ensemble derived from the Liouville operator Chem. Phys. 2010, 370, 294-305.

19


The Journal of Physical Chemistry APPENDIX A. Dihedral angle functions 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

OPLS and Ryckaert-Bellemans dihedral angles We provide the relations for mapping the coefficients Fi and Gi to Cn in terms of φ for the OPLS and RB functions given by equations (2) and (4), respectively. By using trigonometric identities for cos(iφ) the results are, C0 = F0 + F1 + 2F2 + F3 C1 = F1 − 3F3 C2 = −2F2 C3 = 4F3 C4 = 0 C5 = 0 We note that with these values of Cn , ψtrans =0, therefore, every coefficient Cn have to be multiplied by (-1)n to give ψtrans =180. These new values of Cn were used in the topology file to perform a simulation with the gromacs package. The relations between the coefficients Gi from equation (3) with those of the RB function Cn in terms of φ are, C0 = G0 − G2 + G4 − G6 C1 = G1 − 3G3 + 5G5 C2 = 2G2 − 8G4 + 18G6 C3 = 4G3 − 20G5 C4 = 8G4 − 48G6 C5 = 16G5 The Cn coefficients have to be multiplied by (-1)n to give ψtrans =180. These new values of Cn were used in the topology file to perform a simulation with the gromacs package. In addition, to check consistency between OPLS and RB functions, NPT MD simulations are performed at 298.15 K and 1 bar on 1000 methyl tert-butyl ether molecules. The OPLS parameters in kJ/mol taken from Ref.24 for the angle CH3-C-O-CH3 are F0 = 0, F1 =6.03, F2 =-1.36 and F3 =4.64. In order to be used in the simulation with the Gromacs package, these values have to be scaled by 2 to obtain the Fig parameters, see the manual.31 The RB coefficients Cn , calculated with Fi values from the OPLS function in terms of φ and scaled by (−1)n , are C0 = 7.95, C1 =7.89, C2 =2.72, C3 =-18.56, C4 = C5 =0. We used a time step of 0.002 ps with τT = 0.5 ps and τP = 1.0 ps for the the Nosé-Hoover thermostat and ParrinelloRahman barostat, respectively, see Appendix B. The dihedral angle distributions from OPLS and RB functions are the same as shown in Fig. 7. As expected, the trans and ±gauche are the stable conformations. The results of liquid density, dielectric constant, heat of vaporization and compressibility factor are given in Table 1. They are the same, within simulation error.

20


Page 20 of 32

Page 21 of 32 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Figure 7: Dihedral angular distribution for the angle CH3 − C − O − CH3 of methyl tert-butyl ether at 298.15 K and 1 bar using the OPLS (dashed line) and Rycakert-Bellemans (continuous line) functions. The diehedral coefficients used in equations (2) are taken from Ref.24

Property ρ(kg/m3 ) ∆H(kJ/mol) κT (1/GP a)

OP LS 748.6 4.6 27.7 1.91

RB 748.4 4.5 27.8 1.91

Table 1: Results for methyl tert-butyl ether at 298.15 K and 1 bar using the OPLS and Ryckaert-Bellemans functions for the dihedral potential. The estimated errors are: 0.3 kg/m3 for density, 0.3 for dielectric constant, 0.2 kJ/mol for ∆Hvap and 0.03 GPa−1 for κT .

21


APPENDIX B. Simulation details

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Page 22 of 32

The TraPPE UA force field was parametrized using MC methods in different ensembles (GEMC, GCMC and NPT). Therefore it is important to discard methodological aspects with the molecular dynamics results from this work. Testing the correct solution of the equation of motion in molecular dynamics is not easy because the algorithm used in the present calculations does not have a constant of motion40 that can be followed to choose the optimum time step, thermostat and barostat parameters. In this Appendix we analyze the effect that system size, time step, thermostat and barostat parameters have on the MD simulations results. 1. System size NPT MD simulations at 298.15 K and 1 bar are performed on methyl tert-butyl ether with a time step of 0.002 ps, τT = 0.5 ps and τP = 1.0 ps to analyze the effect that number of molecules has on several properties. It is seen that results for 500 and 1000 molecules are the same within simulation error, see Table 2.

Number of molecules 250 500 1000

ρ [kg/m3 ] 748.7 748.8 748.4

Dielectric constant 4.6 4.5 4.5

∆Hvap [kJ/mol] 28.2 27.9 27.8

κ [1/GPa] 1.89 1.94 1.91

Table 2: Results of NPT MD simulations for methyl tert-butyl ether using different number of molecules. The estimated errors are: 0.5 kg/m3 for density, 0.3 for dielectric constant, 0.2 kJ/mol for ∆Hvap and 0.04 GPa−1 for κ.

2. Time step, thermostat and barostat parameters Initial NPT MD calculations for ethyl acetate and vinyl acetate at 298.15 K and 1 bar with parameters taken from Ref.27 gave us a liquid density around 1100 kg/m3 for both molecules compared with the reported value of around 900 kg/m3 at 300 K, see Fig. 2 of Ref.27 The difference was due to a missprint in the sign of coeffcient C1 of the dihedral angle CH3-O-C=O in Table 3 of that reference, instead of 2194 K it should be -2194 K.? Both molecules contain that angle. The MD results obtained with the new parameter were 881.2 kg/m3 for ethyl acetate and 920.6 kg/m3 for vynil acetate. The results obtained with NPT MC simulations? at 300 K and 1 bar were slightly smaller, 873.33 kg/m3 and 904.06 kg/m3 , respectively. It is clear that changes in the coefficients of the dihedral angle lead to a great change in the liquid density. After that finding, we decided to make a careful analysis to determine if also the time step, thermostat and barostat parameters have an effect on several properties of ethyl acetate. The results of density, dieletric constant, compressibility factor and heat of vaporization are given in Table 3 for the Berendsen termostat and barostat. It is found that ρ, and ∆Hvap are almost independent of those parameters. The main change is observed in κT where the property does not converge even at a large value of τP =5. The Berendsen barostat does not sample properly the volumen fluctuations. The results for the same properties are given in Table 4 for the Nosé-Hoover thermostat and Parrinello-Rahman 22


Page 23 of 32 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

∆t (ps)

τT (ps)

0.001

0.5

0.001

1.0

0.002

0.5

0.002

1.0

The Journal of3 )Physical τP (ps) ρ (kg/m Chemistry κT (1/GPa) 0.5 880.7 2.3 0.37 1.0 881.6 2.3 0.27 2.5 881.2 2.3 0.20 5.0 881.1 2.3 0.15 0.5 881.2 2.3 0.31 1.0 881.5 2.3 0.23 2.5 881.6 2.3 0.19 5.0 881.4 2.3 0.12 0.5 881.5 2.4 0.34 1.0 881.3 2.4 0.29 2.5 881.0 2.3 0.18 5.0 880.1 2.3 0.12 0.5 881.7 2.3 0.32 1.0 881.6 2.3 0.27 2.5 881.2 2.3 0.19 5.0 881.2 2.3 0.09

∆Hvap (kJ/mol) 33.3 33.2 33.2 33.2 39.2 33.3 33.4 33.3 33.3 33.4 33.4 33.4 33.4 33.4 33.4 33.4

Table 3: Calculated properties of ethyl acetate for different values of ∆t, τT and τP using the Berendsen barostat. The estimated errors are: 0.5 kg/m3 for density, 0.3 for dielectric constant, 0.2 kJ/mol for ∆Hvap and 0.04 GPa−1 for κ. barostat. All the properties converge to the same value, within simulation error, when τP =1.0 or larger. The results of ρ, and ∆Hvap are the same as those obtained with the Berendsen thermostat and barostat while the value of κT is around 5 times larger. The same effect of simulation parameters is expected to be found in other molecules. The values used in all the simulations performed in this work were ∆t=0.002 ps, τT =0.5 ps and τP =1.0 ps. The Nosé-Hoover thermostat does not give the correct canonical distribution of velocities40 and therefore it might fail to give the correct NPT ensemble when it is coupled with the Parrinello-Rahman barostat. The Nosé-Hoover chain method gives the correct NVT and NPT ensembles67, 68 for systems with and without holonomic constraints. That method with holonomic constraints has not been implemented in the Gromacs package. In the present study the bond distances are kept constant. We have developed a simulation with an own code using the NPT Nosé-Hoover chain method for 500 molecules of ethyl acetate at 300 K and 1 bar with the thermostat and barostat parameters of 0.5 ps and 2.0 ps, respectively to compare the volume fluctuations with those from Nosé-Hoover/Parrinello-Rahman methods. The calculation was performed using a chain of 3 thermostats per molecule and a global barostat. The bond distances were kept constant to the equilibrium value using a rolled SHAKE-RATTLE algorithm.68 The values of density and isothermal compressibility were 877 ± 7 kg m−3 and 1.10 GPa−1 , respectively, in good agreement with those given in Table 4.

23


The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

∆t (ps)

τT (ps)

0.001

0.5

0.001

1.0

0.002

0.5

0.002

1.0

τP (ps) 0.5 0.8 1.0 2.5 5.0 0.5 0.8 1.0 2.5 5.0 0.5 0.8 1.0 2.5 5.0 0.5 0.8 1.0 2.5 5.0

ρ (kg/m3 ) 1101.0 1093.9 881.0 880.9 880.9 1101.6 881.6 881.1 880.9 880.9 1099.9 1090.4 881.6 881.0 880.9 1098.3 1086.4 881.2 881.1 881.2

1.1 1.1 2.3 2.3 2.3 1.1 2.3 2.3 2.3 2.3 1.2 1.3 2.3 2.3 2.3 1.1 1.3 2.3 2.3 2.3

κT (1/GPa) 6.34 1.32 1.13 1.09 1.16 6.51 1.11 1.10 1.22 1.12 2.08 1.23 1.11 1.09 1.18 2.12 1.21 1.12 1.13 1.09

Page 24 of 32

∆Hvap (kJ/mol) 39.4 39.4 33.2 33.2 33.2 39.7 33.5 33.5 33.5 33.5 40.2 40.6 33.5 33.5 33.5 38.5 38.8 33.6 33.6 33.6

Table 4: Calculated properties of ethyl acetate for several values of ∆t, τT and τP using the ParrinelloRahman barostat. The estimated errors are: 0.5 kg/m3 for density, 0.3 for dielectric constant, 0.2 kJ/mol for ∆Hvap and 0.04 GPa−1 for κ.

24


Page 25 of 32 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45



The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46


Page 26 of 32

Page 27 of 32 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44



The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46


Page 28 of 32

Page 29 of 32 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44



The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44


Page 30 of 32

Page 31 of 32


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51


Page 32 of 32

Force Field Benchmark of the TraPPE UA for Polar Liquids: Density

Recommend Documents