FULL PAPER

WWW.C-CHEM.ORG

A Polarizable Empirical Force Field for Molecular Dynamics Simulation of Liquid Hydrocarbons Oliwia M. Szklarczyk, Stephan J. Bachmann, and Wilfred F. van Gunsteren* Electronic polarizability is usually treated implicitly in molecular simulations, which may lead to imprecise or even erroneous molecular behavior in spatially electronically inhomogeneous regions of systems such as proteins, membranes, interfaces between compounds, or mixtures of solvents. The majority of available molecular force fields and molecular dynamics simulation software packages does not account explicitly for electronic polarization. Even the simplest charge-on-spring (COS) models have only been developed for few types of molecules. In this work, we report a polarizable COS model for cyclohexane, as this molecule is a widely used solvent, and for linear alkanes, which are also used as solvents, and are the precursors of lipids, amino acid side chains, carbohydrates, or nucleic acid backbones. The model is an extension of a nonpolarizable united-atom model

for alkanes that had been calibrated against experimental values of the density, the heat of vaporization and the Gibbs free energy of hydration for each alkane. The latter quantity was used to calibrate the parameters governing the interaction of the polarizable alkanes with water. Subsequently, the model was tested for other structural, thermodynamic, dielectric, and dynamic properties such as trans/gauche ratios, excess free energy, static dielectric permittivity, and self-diffusion. A good agreement with the experimental data for a large set of properties for each considered system was obtained, resulting in a transferable set of polarizable force-field parameters for CH2, C 2014 Wiley Periodicals, Inc. CH3, and CH4 moieties. V

Introduction

standardly used in biomolecular simulation studies. They also require the adaptation of the corresponding simulation software packages.[13,14] When extending a nonpolarizable model or a force field with atomic polarizability one would like to keep both types of force fields as compatible with each other as possible, such that parts of a system may be treated with explicit and other parts with implicit electronic polarizability. Based on this requirement, a number of polarizable models for solvents such as water,[15–17] methanol,[18] ethylene glycol,[19] chloroform,[20] and carbontetrachloride,[21] have been developed that are compatible with the GROMOS nonpolarizable biomolecular force fields.[22,23] Because of the abundance of polarizable electronic degrees of freedom in the solvent part of biomolecular systems such as protein or DNA in aqueous solution, polarization of these degrees of freedom can hardly be modeled quantum chemically, which leaves the introduction of classical atomic or molecular polarizability as the only road to a more faithful representation of polarization properties of the solvent. For solute degrees of freedom, the use of quantum-chemical models, as in quantum-mechanical (QM)/molecular-mechanical simulations, is a viable alternative to represent electronic

Simulation of the motions in molecular systems has become a standard method of investigation during the past decades, in particular for biomolecular systems.[1,2] These motions are governed by interactions between nuclei and electrons and should as such be treated using quantum-chemical methodology. Due to their computational complexity, quantumchemical models can only be applied to relatively small systems and for short simulation periods. This led to the use of empirical interaction functions for the atoms of a system and simulation based on classical, Newton’s equations of motion. These empirical interaction functions or force fields generally assign fixed values of partial charges to atoms in such a way that the average effect of the eliminated electronic degrees of freedom is accounted for. Such an implicit treatment of electronic polarization may well be justified for spatially homogeneous systems such as neat liquids but may be a poor approximation when considering spatially electronically inhomogeneous systems containing proteins, membranes, interfaces between different compounds, or mixtures of solvents. Protein stability and function is influenced by the polarizable environment.[3–9] The polarizability affects the conformational properties of peptides near membranes[10] and plays a role in ion transport through channels. The need to account explicitly for electronic polarizability in empirical force fields has long been recognized but the development of polarizable force fields that are applicable to many different compounds, in particular biomolecules, is a daunting task. Only during the past decades, general biomolecular force fields that include atomic polarizabilities have become available.[11,12] These are not yet extensively tested and are not yet

DOI: 10.1002/jcc.23551

O. M. Szklarczyk, S. J. Bachmann, W. F. van Gunsteren Laboratory of Physical Chemistry, Swiss Federal Institute of Technology ETH, 8093, Z€ urich, Switzerland E-mail: [email protected] Contract grant sponsor: National Center of Competence in Research (NCCR) in Structural Biology; Contract grant sponsor: Swiss National Science Foundation; Contract grant number: 200020-137827; Contract grant sponsor: European Research Council; Contract grant number: 228076 C 2014 Wiley Periodicals, Inc. V

Journal of Computational Chemistry 2014, 35, 789–801

789

FULL PAPER

WWW.C-CHEM.ORG

polarization. In biomolecular systems, the next most abundant type of molecules, apart from the solvent, are lipids, which make these molecules the next type for which a polarizable model compatible with the GROMOS biomolecular force fields is to be developed. As aliphatic tails constitute a major part of lipids and membranes, their force-field parameters can be developed using the measured thermodynamic and dielectric properties of liquid alkanes as calibration data. There are many force fields available for n-alkanes and cycloalkanes. Most of them use a fixed-charge model for the treatment of electrostatic interactions, and thus, neglect the effect of electronic polarization and the dielectric properties of alkanes. The ones treating electronic polarizability explicitly use one of the three different methods to introduce electronic polarizability explicitly in molecular simulations.[24] In the fluctuating charge (FQ) or chemical potential equalization, or electronegativity equalization model,[25–27] the polarizability is introduced by changing the size of the atomic partial charges in response to changes in the local electric field. The FQs can be assigned fictitious masses and treated as additional dynamic degrees of freedom, or their values may be obtained at every time step of the simulation by minimization of the electrostatic energy of the system for the given configuration. The FQ approach was applied to obtain a polarizable model for aliphatic hydrocarbons in.[28–31] In the inducible point-polarizable dipole moment (PPD) approach,[32,33] each polarizable atom or site is assigned an induced ideal dipole moment proportional to the electric field at the site and to the atomic or site polarizability. The electric field can be determined through a self-consistent field calculation in which the field equations are solved iteratively until a certain level of convergence is achieved. An alternative method to compute the electric field is the extended Lagrangian method,[34] in which the polarizable degrees of freedom are included as dynamic variables. An example of a polarizable force field for alkanes based on the PPD model can be found in reference.[35] The charge-on-spring model (COS) (or Drude’s oscillator, or shell model)[36,37] accounts for electronic polarizability of atoms by introducing a massless virtual charge on a spring at each polarizable atom or site which carries in addition to its fixed partial charge, a charge of opposite sign to that of the virtual charge. The force constant of the spring is determined by the given atom or site polarizability, and charge and virtual charge do not interact with each other through Coulomb forces, only through the spring forces. The COS model has the very useful feature of avoiding the complex evaluation of dipole–dipole interactions—only pairwise Coulombic interactions are calculated. Thus, it is straightforwardly implemented in most (bio-)molecular force fields and software. The computational cost of the additional charge–charge interactions can be kept low by assigning COS-sites to a subset of atoms, for example, only to atoms that are strongly polarizable, as is done in the work presented here. The induced COS dipoles are calculated from the total electric field at their charge sites, which is determined self-consistently at each time step.[24] Only one example of a polarizable force field using the COS 790

Journal of Computational Chemistry 2014, 35, 789–801

model for alkanes is known to us.[38] Some properties of short, up to six carbon atoms, n-alkanes, cyclopentane, and cyclohexane, obtained using this COS polarizable Charmm model are reported in.[39–41] In the present work, we propose a polarizable COS force field for liquid alkanes compatible with the GROMOS force fields as a step toward the development of such a polarizable model for lipid membranes. We follow the standard strategy as used for the determination of the parameters for the other compounds of the GROMOS force field as described in.[42] These force fields treat aliphatic CHn groups as united atoms, which for alkanes have a permanent partial charge equal to zero. Their van der Waal’s interaction parameters had been determined such as to reproduce measured values of various structural and thermodynamic quantities of alkanes in the condensed phase.[43] Atomic polarizabilities were added to the carbon united atoms and subsequently the C12 van der Waal’s parameters governing the interaction with simple point charge (SPC) water molecules[44] were varied to reproduce measured values of the free enthalpy of hydration of the alkanes. Most of the properties calculated are in very good agreement with the experimental data. The set of parameters for the aliphatic groups or united atoms CH2, CH3, and CH4 can be straightforwardly used in a polarizable model for lipids, amino acids, and other larger molecules.

Methodology The nonpolarizable GROMOS force-field parameters for the CH2, CH3, and CH4 united atoms of n-alkanes can be found in,[43,45] and the nonpolarizable water model used is the SPC model.[44] The atomic polarizabilities ai for the carbon atoms were taken from.[46] Explicit electronic polarization of n-alkanes was modeled by a simple COS on the carbon united atoms, that is, without off-atom sites[24] or damping of the electric field.[17] The virtual charge on the spring was 28e, as in.[15,37] As a polarizable model for liquid water the COS/G2 off-site model[16] was used. The number of alkane molecules per system was kept at 512; thus, the number of atoms and effectively the computational cost varied between the systems. For the calculations of the free enthalpy of hydration, we used one alkane molecule solvated in a cubic box of water molecules. The number of water molecules in the box varied between the systems and was determined such as to prevent the periodic images from interacting with each other. We used 999 water molecules for methane, 1500 for butane, 2500 for hexane and cyclohexane, 3300 for nonane, 4500 for dodecane, and 6000 for pentadecane. All simulations were performed with the aid of the GROMOS simulation package.[13,45,47] For all systems considered, the pressure was set to 1 atm and the temperature to 298.15 K. Only for methane and butane, the reference temperature was set to 111.67 K and 272.66 K, respectively, which are the boiling point temperatures of these compounds. These temperatures are referred to as reference temperatures Tref in this work. Bond lengths in the alkanes and the geometry of the water molecules were kept rigid using the SHAKE method[48] with a geometric tolerance of 1024. WWW.CHEMISTRYVIEWS.COM

FULL PAPER

WWW.C-CHEM.ORG

Table 1. Summary of the experimental values of properties analysed for six n-alkanes (C1 to C15) and cyclohexane (c-C6), the subscript indicating the number of carbon atoms. Values at 1 atm and 298.15 K except for methane (111.67 K) and butane (272.66 K). Static dielectric permittivities of liquids at 293 K and 1 atm. Property

C1

C4

C6

C9

C12

C15

c-C6

rt (%) q (kg m23) DHvap (kJ mol21) DGhydr (kJ mol21) DFexc (kJ mol21)

— 422.6[a] 8.2[a] 8.1[e] 21.8[g] 4.5[f ] 13.3[a] 55.8[a] 22.7[a] 3.5[a] 1.632[b] 1.6[b]

60.0[d] 600.7[a] 22.4[a] 9.0[e] 10.1[g] 11.8[f ] 14.9[a] 134.2[a] 21.9[a] 1.8[a] 1 .770[b] 6.0[i] 5.00[j] 0.20[a]

64.0[d] 654.8[a] 31.6[a] 10. 7[e] 14.4[g] 17.0[f ] 17.9[a] 194.0[a] 17.9[a] 1.4[a] 1.882[b] 4.2[c]

65.0[d] 713.8[a] 46.4[a] 12.8[f ] 22.5[g] 25.0[f ] 22.4[a] 283.1[a] 12.1[a] 1.1[a] 1.960[b] 1.7[c]

65.0[d] 745.2[a] 61.3[a] 15.0[f ] 30.9[g] 33.3[f ] 24.9[a] 376.8[a] 10.1[a] 1.0[a] 2.006[b] 0.8[h]

65.0[d] 764.9[a] 76.2[a] 17.3[f ] 39.7[g] 42.3[f ] 26.7[a] 469.8[a] (sat.pressu.) 9.2[a] — 2.033[b] —

94.0[l] 773.9[a] 33.0[b] 5.2[e] 16.0[g] 18.5[f ] 24.6[a] 153.0[a] 11.5[a] 1.2[a] 2.016[b] 1.4[h]

0.30[a]

0.65[a]

1.36[a]

2.53[a]

0.88[a]

c (mN m21) Cp (J K21 mol21) jT (1021 atm21) aT (K21) (0) D(1025 cm2 s21)

0.12[a]

g(cP) [50,51]

[52]

[70]

[a] from reference [b] from reference [c] from reference [d] from reference [e] from reference [f ] calculated results from reference[72] [g] calculated results from vapor pressures[50] [h] from reference[73] [i] calculated value from reference[74] [j] calculated value from reference[35] [k] calculated value from reference[75] [l] calculated value from reference[76]

Liquid-phase simulations For the pure liquid-phase simulations, 512 single n-alkane molecules in the all-trans configuration were randomly placed in a cubic box at the experimental density of each system (Table 1). The 512 single cyclohexane molecules were in the chair configuration. The systems were energy minimized using minimum image periodic boundary conditions and initial velocities were assigned from a Maxwell distribution at the reference temperature Tref . To obtain a correct ratio of trans versus gauche conformations in the aliphatic n-alkane chains, we performed a pre-equilibration of the systems by carrying out constant volume molecular dynamics (MD) simulations for 10 ns omitting the torsional-angle potential-energy term of the force field which allowed faster transition from the all-trans starting configuration to a configuration with both trans and gauche molecular conformations. The last configuration from this preequilibration at temperature Tref was used as a starting structure for the pure liquid-phase simulations with the torsionalangle potential-energy term switched back on. The system was further equilibrated for 1 ns at a constant temperature Tref and a pressure of 1 atm before performing the liquid phase MD production runs of 10 ns under these isothermalisobaric conditions. The systems were weakly coupled to a temperature bath[59] with a coupling time of 0.1 ps, separately for the translational and rotational degrees of freedom. The pressure was kept at 1 atm by weak coupling to a pressure bath with a coupling time of 0.5 ps. The isothermal compressibility of each system was set to the experimental value (Table 1).[49,50] The equations of motion were integrated with the leap-frog algorithm and a time step of 2 fs. Triple range cutoff radii of 0.8/1.4 nm were used to treat van der Waal’s and electrostatic interactions. The interactions within the shortrange cut-off were calculated every time step, whereas the interactions in the intermediate range were calculated every five steps, from an updated pairlist. The long-range interactions were approximated by a continuum reaction field with a

[43]

[71]

dielectric permittivity of the continuum equal to the experimental dielectric permittivity of the specific alkane system.[51] The configurations and energies during the production runs were saved every 2 ps. Gas-phase simulations The initial configurations for the gas-phase simulations were taken from the last configuration of the liquid-phase simulations. The molecules were set 50 nm apart from each other and allowed to equilibrate for 100 ps before performing the gas-phase stochastic dynamics (SD) production runs. The gas phase was simulated using a SD algorithm for the integration of the equations of motion.[60] The friction coefficient was set to 9.1 ps21. The total simulation time was 10 ns for each alkane system, from which the average gas-phase energies of the molecules could be calculated. Free energy calculations The Helmholtz excess free energy of liquid n-alkanes was calculated using thermodynamic integration (TI) by switching off all nonbonded interactions. The system was changed from a fully interacting pure liquid of 512 alkane molecules in a cubic box (k 5 0) to 512 molecules without nonbonded interactions (k 5 1), at constant temperature and volume. This was followed by an analogous procedure in the gas phase for a single alkane molecule undergoing SD. The excess free energy DFexc was calculated as  DFexc 5

1

Nmol

ð1 dk 0

D @HðkÞ E @k

 ð1 D @HðkÞ E : 2 dk liq @k gas 0

(1)

The water boxes for the calculations of the free enthalpy of hydration consisted of a number of water molecules large enough to prevent the periodic images of the solute from interacting with each other. The systems were energy Journal of Computational Chemistry 2014, 35, 789–801

791

FULL PAPER

WWW.C-CHEM.ORG

minimized, then thermalized at 60 K followed by five 100-ps equilibration periods at constant temperature and volume, after each of which the reference temperature was raised by 60 K till 298.15 K. Subsequently, the systems were equilibrated for 200 ps at constant temperature and pressure using a compressibility of 4.575 3 1024 kJ21 mol nm3 and a value of 87.8 for the permittivity of the reaction field, the dielectric permittivity of the COS/G2 water model,[16] in the COS simulations and a value of 66.6, the permittivity of the SPC water model,[61] in the nonpolarizable model simulations. An energy minimized alkane molecule was then solvated in the preequilibrated water box. The solvated systems were minimized, and then, thermalized and equilibrated analogously to the water boxes, with the positions of the atoms of the alkane molecule harmonically restrained with a force constant of 2.5 3 104 kJ mol21 nm22. The restraining force was gradually reduced to zero while thermalizing. The Gibbs free energy or free enthalpy of hydration DGhydr was calculated using TI. The system was changed from one alkane molecule in a box of water molecules (k50) to one alkane molecule without nonbonded interactions in a box of water (k51) under isothermal-isobaric conditions. The dependence of the Hamiltonian upon the coupling parameter used to switch off or on particular interactions is given in.[13,62] To account for the intramolecular nonbonded interactions, the free enthalpy change of switching off the nonbonded interactions for a single molecule in the gas phase was calculated. The hydration free enthalpy DGhydr was calculated as: ð 1 DGhydr 5

dk 0

D @HðkÞ E @k

ð1 D @HðkÞ E  : 2 dk gas @k liq 0

(2)

The Gibbs free energy or free enthalpy of solvation DGsolv of a methanol molecule in a box of alkane molecules was also calculated using TI. The system was changed from one methanol molecule in a box of alkane molecules (k50) to one methanol molecule without nonbonded interactions in a box of alkane molecules (k51) under isothermal-isobaric conditions. In the COS simulations, a polarizable COS/M methanol molecule[18] was solvated in a box of polarizable alkanes, and in the nonpolarizable simulations, a nonpolarizable methanol molecule was solvated in a box of nonpolarizable alkanes. Boxes of 999 methane molecules and 512 molecules for other alkane systems were used to solvate methanol. The systems were constructed, thermalized, and equilibrated analogously to the systems in the TI calculations of DGhydr. The isothermal compressibility and the dielectric permittivity of the medium outside the cut-off radius were set to the experimental values of these properties for each alkane system. For methane and butane, the simulations were performed at the boiling point temperatures of these alkanes. For all the TI simulations, 21 k-points were used, 0.05 apart in the range [0, 1]. We used the two-step method for the TI calculations of the polarizable systems, in which the electrostatic interactions (charges and polarizabilities) and the Lennard–Jones interactions were switched off separately. The reason for using a two-step procedure is that when Lennard– 792

Journal of Computational Chemistry 2014, 35, 789–801

Jones interactions are very weak and the repulsion does not keep the atoms well apart, the r 21 Coulomb interaction, if present, may cause huge forces. The procedure was the following. In the first step (first TI simulation from k 5 0 to 1), no soft-core interaction was used (aLJ 50; aC 50), and only the electrostatic interactions were switched off. In the second step (second TI simulation from k 5 0 to 1), the Lennard–Jones interactions were switched off with a softness parameter aLJ 50:5.[63,64] After integrating over h@H=@ki, the two integrals were added. In some cases, a simpler one-step procedure could be used, in which the Lennard–Jones and Coulomb interactions were switched off in parallel, using the softness parameters aLJ 50:5 and aC 50:5 nm2 . We used the one-step method for the TI calculations of the nonpolarizable systems, and the two-step TI method for the calculation of the free energies and enthalpies of the polarizable models, because of the larger electrostatic forces in the COS model. TI calculations of the DGhydr for the COS systems were performed with the GROMOS96 Fortran code.[60] In the liquidphase simulations, at each k-point 100 ps equilibration and up to 400 ps sampling were performed. TI calculations of the D Ghydr for the nonpolarizable systems were performed with the GROMOS11 code.[13,45,47] In the liquid-phase simulations, 100 ps of equilibration and 2 ns of sampling were performed at each k-point. The systems were weakly coupled to a temperature bath[59] with a coupling time of 0.1 ps, separately for the translational and rotational degrees of freedom. Additionally, Langevin coupling was used with a friction coefficient of 9.1 ps21. In case of the DGsolv of methanol in alkanes, all simulations were performed with the GROMOS11 code, with 100- ps equilibration and 1-ns production time for each k-value. The systems were weakly coupled to a temperature bath[59] with a coupling time of 0.1 ps, separately for the solute and solvent. All the gas-phase TI calculations were performed for a single molecule coupled to a Langevin thermostat with temperature set to 298.15 K and the friction coefficient set to 9.1 ps21. The setup was as for the liquid phase except that at each k-point 100-ps equilibration and 100-ns sampling were performed, and the dielectric permittivity of the reaction field was set to 1.

Analysis Trans/gauche ratios hri. The spatial structure of aliphatic chain molecules is determined by the energy difference between the gauche and trans configurations reflected by the ratio of trans to gauche configurations. The geometric definition of the percentage of trans configurations of torsional angle / was calculated as[65]:

/5½2180o ; 1180o 

X

100  n/ X t ; hri5 X n/g n/t 1

for

/t 5j/j > 120o ;

(3)

/g 5j/j < 120o where n/t and n/g are the numbers of trans and gauche configurations, respectively. WWW.CHEMISTRYVIEWS.COM

FULL PAPER

WWW.C-CHEM.ORG

Table 2. Simulation set-up for the calculation of the heat capacity, the thermal expansion coefficient (NPT simulations), and the isothermal compressibility (NVT simulations). Results for the alkane liquids are shown, as well as the empirical correction to the heat capacity, dCp. Results from the NPT simulations for ethane and propane were used to derive the heat capacity correction. The final averaged values are shown in Table 7. NPT

Alkane system C1

C2

C3

C4

C6

C9

C12

C15

c-C6

NVT

Csim P 21 21

(K)

(kg m23)

90.9 111.5 70.9

454.8 424.0 513.3

44.7 63.0

2.0

163.9 183.7 143.8

573.5 549.0 597.1

58.0 57.8

15.4

210.6 230.9 190.7

604.1 582.1 625.3

72.9 73.5

25.8

248.7 270.6 228.0

628.6 605.9 649.1

88.8 89.3

317.5 297.9 277.5

639.5 656.4 673.2

144.0 144.0 142.0

337.6 317.6 297.7

684.7 699.4 713.6

337.68 317.6 297.8

aT (1021 atm21)

(atm)

q (kg m23)

jT (K21)

3.3 4.3

0.8 214.1 268.5

424.0 442.6 446.6

20.6 20.0

36.2

1.7 1.6

1.3 165.3 387.8

605.9 620.7 640.7

14.9 14.9

57.0

1.3 1.2

1.2 191.0 461.8

656.4 674.8 694.8

12.7 14.8

212.9 214.3

88.2

1.0 1.0

1.5 279.4 630.2

713.6 733.8 753.8

9.0 10.2

717.4 731.0 744.0

286.5 285.5

119.4

0.9 0.9

2.3 297.7 735.3

744.0 765.2 785.2

7.6 8.6

337.5 317.7 297.7

737.8 750.4 763.1

361.1 359.8

150.0

0.8 0.8

0.0 373.9 809.8

763.1 784.9 804.9

6.8 7.7

317.6 297.8 277.4

752.3 768.5 785.0

120.5

62.4

1.1 1.0

1.3 346.5 686.1

768.5 793.9 813.9

8.6 9.6

Enthalpy of vaporization DHvap .

(JK

mol

dCP (JK21 mol21)

)

The free enthalpy of vaporiza-

tion DHvap was calculated as: DHvap ðTÞ5

1 pot 1 pot U ðTÞ2 U ðTÞ1RT; Nmol gas Nmol liq

(4)

where R is the gas constant, T is the temperature, and Upot gas and Upot liq denote the total potential energies in the gas phase and in the liquid phase, respectively. Self-polarization energy Upol . The self-polarization energy Upol

accounts for the cost of the dipole induction and is calculated as follows[66]: Upol 5

NCOS 2 1X ðlind i Þ : 2 i51 ai

lind i

(5) qvi

are the charges multiThe induced dipole moments plied by the displacements Driv of the virtual COS particles. ai denotes the electronic polarizability assigned to the COS-site, that is, carbon united atom, and NCOS is the number of COSsites in the system. Surface tension c. To calculate the surface tension c, we performed five additional NVT simulations with varying random number seeds. A special box was constructed by increasing

the length of the box obtained from the NPT liquid-phase simulations by 45 nm in the z-direction. Each simulation was 1 ns long (preceded by a 100-ps equilibration). The surface tension c was obtained from the diagonal elements of the pressure tensor as[67]: c5

  Lz 1 pzz 2 ðpxx 1pyy Þ ; 2 2

(6)

where Lz is the box length in the z-direction, and pab is the pressure tensor element. The resulting values of the surface tension were then averaged over the five runs. Heat capacity Cp. The isobaric molar heat capacity was calculated using a finite-difference approximation,[68,69]

Cp 5Cpsim 1dCp ;

(7)

where Cpsim is calculated from the simulations as Cpsim 5

Utot;2 2Utot;1 : T2 2T1

(8)

Utot is the total energy per mole, T is the temperature, and d Cp is an empirical correction to the heat capacity calculated per functional group. To evaluate the first term in eq. (7), for Journal of Computational Chemistry 2014, 35, 789–801

793

FULL PAPER

WWW.C-CHEM.ORG

each system 5 ns constant pressure simulations (preceded by a 1-ns equilibration) at two additional temperatures were performed with DT5620 K (Table 2). The temperatures were chosen such that the considered systems would remain in the liquid phase. Thus, every system was simulated at three different temperatures, from which two values of the heat capacity were obtained. The final value of the first term in eq. (7) was the average of the two intermediate values. The second term in eq. (7), dCp , is a first-order empirical correction for the neglect of hydrogen atoms in the united-atom model, for the rigid treatment of bonds within a molecule, as well as for the neglect of QM features of the liquid in the classical approximation. To evaluate the heat capacity correction term, liquid-phase simulations for ethane and propane were performed at three different temperatures, as described above. The correction to the heat capacity is defined per CHn group and calculated from the difference between the experimental and simulated values,

0 q 1   ln q2 1 @V 1 A :  2@ aT 5 T2 2T1 V @T p

The densities and temperatures were calculated from the same constant pressure simulations as used for the calculation of heat capacity. Thus, the final value of the thermal expansion coefficient was averaged over the two intermediate values, as for the heat capacity and the isothermal compressibility.

The average induced molecular dipole moment lind;mol was calculated as follows: Average induced molecular dipole moment lind;mol .

lind;mol 5

Nmol 1 X

Nmol n51 NX COS;n

dCp ðCH4 Þ5Cpexp 2Cpsim

(9)

from methane simulations,

(13)

p

lind;n 5

jlind;n j;

(14)

0

qvi Dri ;

(15)

i51 0

1 dCp ðCH3 Þ5 ðCpexp 2Cpsim Þ 2

(10)

where qvi denotes the charge of the COS-site, Dri is the displacement of the virtual charge qvi from the COS-site, and the sum in eq. (15) is over the number of COS-sites in molecule n.

from ethane simulations, and Static dielectric permittivity 僆(0). The static dielectric permit-

dCp ðCH2 Þ5ðCpexp 2Cpsim Þ22dCp ðCH3 Þ

(11)

from propane simulations.

tivities E(0) were calculated from the linear response to an applied electric field.[71] External electric fields of different strengths were applied to the system and the permittivity was calculated from the linear regression of the box polarization Pz calculated as

Isothermal compressibility jT . To obtain the isothermal com-

pressibility jT , two additional constant volume simulations were performed, each with a different density of the system with Dq5620 kg m23 (Table 2). The densities were chosen such that the systems stayed in the liquid phase. The systems were simulated for 5 ns (preceded by a 1 ns NVT equilibration). Assuming a linear dependence between two thermodynamic variables and having simulated the system at two different thermodynamic states, jT was calculated from the following equation[70]: 0 q 1       ln q2 1 @V 1 @q @lnq 1 A ; jT 52 5 5 @ p2 2p1 V @p T q @p T @p T

(12)

Pz 5V 21

Nmol NX COS;n X

ð0Þ5114p

794

Journal of Computational Chemistry 2014, 35, 789–801

(17)

calculated[72] according to the Einstein relation

t!1

sion coefficients aT , an expression equivalent to eq. (12) was used,

: Ezext

Self-diffusion constant D. The self-diffusion constant D was

D5 lim Thermal expansion coefficient aT . To obtain the thermal expan-

(16)

where lind,n,i,z is the induced dipole moment along the z-direction of COS-site i in molecule n, and V is the volume of the computational box, versus the applied electric field strength Ezext . Twelve different field strengths Ezext in the z-direction were used: 0.00, 0.01, 0.02, 0.03, 0.05, 0.1, 0.2, 0.5, 1.0, 2.0, 5.0, 7.0, and 10.0 ½e nm22 . The static dielectric permittivity was then calculated from the equation

T

where V denotes the volume, p the pressure, T the temperature, and q denotes the density of the system. Thus, every system was simulated at three different densities at the same temperature, from which two values of the isothermal compressibility were obtained. The final value for the isothermal compressibility was the average of the two intermediate values.

lind;n;i;z ;

n51 i51

hDr 2 ðtÞi ; 2Nd t

(18)

where Nd denotes the number of dimensions and the meansquare-displacement is calculated as WWW.CHEMISTRYVIEWS.COM

FULL PAPER

WWW.C-CHEM.ORG

hDr 2 ðtÞi5

Nat tav 2t 1 X 1 s 5X ðri ðt1sÞ2ri ðsÞÞ2 ; Nat i51 tav 2t s 5 0

where Nat denotes the number of atoms considered, ri is the position of atom i, and tav is the averaging time period. The selfdiffusion constants were calculated from the slope of the meansquare-displacements as a function of time by least-squares fitting.

Table 3. Parameters of the non-polarisable (45A3) and polarisable (COS) alkane models: charge q [e], electronic polarisability aCH2 [(4p0)1023 nm3], Lennard-Jones repulsion parameter [C12]1/2 [1023 (kJ mol21 nm12)1/2] used for the interactions of united atom carbons with other non-polar atoms, and Lennard-Jones dispersion parameter [C6]1/2 [(kJ mol21 nm6)1/2]. The ð2Þ separate Lennard-Jones repulsion parameter [C12 ]1/2 [1023 (kJ mol21 12 1/2 nm ) ] is used for the interactions of united atom carbons with polar atoms, such as the oxygen atom of water.

Shear viscosity g. The shear viscosity g was calculated accord-

Parameter

45A3

COS

ing to the Einstein relation for viscosity using the time integrals of the off-diagonal elements of the pressure tensor pab [72–74]

qCH2 aCH2 [C12(CH2)]1/2 ð2Þ [C12 (CH2)]1/2 [C6(CH2]1/2

0.000 — 5.828 5.828 0.08642

0.000 1.835[a] 5.828 6.29424 0.08642

qCH3 aCH3 [C12(CH3)]1/2 ð2Þ [C12 (CH3)]1/2 [C6(CH3)]1/2

0.000 — 5.162 5.162 0.09805

0.000 2.241[a] 5.162 5.26524 0.09805

qCH4 aCH4 [C12(CH4)]1/2 ð2Þ [C12 (CH4)]1/2 [C6(CH4)]1/2

0.000 — 5.862 5.862 0.1148

0.000 2.593[b] 5.862 5.9788 0.1148

qCH2 r aCH2 r [C12(CH2r)]1/2 ð2Þ [C12 (CH2r)]1/2 [C6(CH2r)]1/2

0.000 — 5.297 5.297 0.08564

0.000 1.835[a] 5.297 5.56185 0.08564

! av 2t X 1 V d 1 s5t 2 gab 5 ðGab ðt1sÞ2Gab ðsÞÞ ; lim 2 kB T t!1 dt tav 2t s50

(19)

(20)

where V denotes the box volume and ðt Gab ðtÞ5 pab ðt0 Þdt 0 :

(21)

0

For isotropic systems, the shear viscosity can be obtained from the average of the three off-diagonal elements of the viscosity tensor, 1 g5 ðgxy 1gxz 1gyz Þ: 3

(22)

To calculate the shear viscosity, additional NVT simulations of 5 ns (preceded by a 100-ps NPT equilibration) were performed at densities obtained from constant pressure simulations at the same temperature. The pressure tensor values were saved every time step. Least-squares fits were performed to calculate gab from the slope of the mean-square displacement (MSD) of Gab ðtÞ as a function of time over the time interval 10–100 ps.

Results The values of model parameters and properties obtained using the GROMOS 45A3 force field are denoted by 45A3 and correspond to the simulations of nonpolarizable alkanes, while the values denoted by COS correspond to the simulations of polarizable alkanes. Because the polarizable and nonpolarizable models only differ in the atomic polarizability parameter ai and the repulsive Lennard–Jones parameter for aliphatic centers interacting with polar atoms such as the oxygen of water, the values of properties for these models will and do only differ for the properties depending on these parameters, that is, ð0Þ and DGhydr . The values for the parameters of the nonbonded interaction of the model are shown in Table 3. The bonded and nonbonded parameters are those for the aliphatic carbon united atoms of the GROMOS 45A3,[43] 53A6,[22] and 54A7[23] force fields. The polarizabilities of the united carbon atoms were set to the theoretical values for the aliphatic groups CH2, CH3, and CH4.[46,75] For cyclohexane, the value for CH2 in linear aliphatic chains was used. The properties of the n-alkane liquids were only calculated for chain lengths 1, 4, 6, 9, 12, and 15 to keep the computa-

[a] from reference[46]; [b] from reference[68]

tional effort limited. As the various properties show smooth behavior as function of chain length, the remaining n-alkanes are expected to behave similarly. The van der Waal’s parameters C12 and C6 for aliphatic united-atom carbons had been determined to reproduce the density and the heat of vaporization of liquid alkanes.[43] The values in Table 4 for the polarizable model are as expected equal to the values for the nonpolarizable model (data not shown) and in agreement with the experimental values, which are reproduced within 1%. Structural properties The structural properties of n-alkanes, apart from the density of the liquid, were studied in terms of the percentage of trans configurations in case of linear alkanes and chair configurations in case of cyclohexane, hri (Table 4). The occurrence of trans configurations is overestimated by approximately 7% for chains longer than butane, compared to the experimental values inferred from NMR experiments. The occurrence of the chair conformation of cyclohexane in the condensed phase agrees with the infrared spectroscopic data within less than 2%. Free enthalpy of hydration, excess free energy, and surface tension. Gibbs free energies of hydration were used as a cali-

bration property for the interaction between CHn groups and polar atoms, such as the oxygen atom of the COS/G2 water model. To reproduce the free enthalpy of hydration of alkanes, Journal of Computational Chemistry 2014, 35, 789–801

795

FULL PAPER

WWW.C-CHEM.ORG

Table 4. Comparison of the properties used in the parameter calibration[43]: density q, heat of vaporisation DHvap, and percentage hri of trans configurations for linear alkanes and chair configurations for cyclohexane. q [kg m23] System

exp

C1 C4 C6 C9 C12 C15 c-C6

422.6 600.7 654.8 713.8 745.2 764.9 773.9

DHvap [kJ mol21]

hri [%]

45A3/COS

exp

45A3/COS

exp

45A3/COS

424.0 605.9 656.4 713.6 744.0 763.1 768.5

8.2 22.4 31.6 46.4 61.3 76.2 33.0

8.2 22.3 31.7 46.7 61.8 76.7 32.9

— 60.0 64.0 65.0 65.0 65.0 94.0

— 58.5 68.1 69.8 70.4 70.8 95.2

new values for the repulsive Lennard–Jones parameters ð2Þ ðC12 ðCHn ÞÞ1=2 of aliphatic atoms for use with those for the oxygen atom of water molecules in the geometric combination rule C12 ðA; BÞ5ðC12 ðAÞÞ1=2 ðC12 ðBÞÞ1=2

(23)

for obtaining Lennard–Jones parameters for pairs of atom types A and B had to be determined (Table 3). Using the ðC12 ðCHn ÞÞ1=2 parameters for the interaction with COS/G2 water resulted in underestimation of the Gibbs free energies of hydration, for example, DGhydr 521:3 kJ mol21 for pentadecane. To make the interaction with water less favorable, the ð2Þ separate ðC12 ðCHn ÞÞ1=2 parameters for the interaction of aliphatic carbons with polar atoms were changed (Table 3). They were increased compared to the nonpolarizable ones by 2% for the CH4 and CH3 groups, by 8% for the CH2 group, and by 5% for the CH2r group. The resulting DGhydr values are listed in Table 5. All the values for the COS model are in good agreement with the experimental ones with deviations ranging from 0.0 to 1.7 kJ mol21. The uncertainty estimates of the hydration free enthalpies for the COS model are larger than those of the simulations of the nonpolarizable systems as the h@H=@ki values were not as well converged due to the shorter simulation time per k-value. The calculated Helmholtz excess free energies of the n-alkanes agree with the experimental values with deviations ranging from 0.1 to 1.7 kJ mol21. The surface tensions are reproduced within 2.8% (Table 5). The calibrated van der Waal’s parameters for the interaction of COS united carbon atoms with polar atoms were tested by calculating the free enthalpies of solvation of a COS/M metha-

nol molecule in a box of COS alkane molecules. The experimental data available for hexane, nonane, and cyclohexane are reproduced within 0.2 kJ mol21 in the simulations using the polarizable model, whereas for the nonpolarizable models the deviations are 1.2, 0.9, and 0.8 kJ mol21. For other alkanes, no experimental values are known to us. However, experimental values of the free enthalpy of solvation of similar small alcohols such as ethanol or propanol indicate that for aliphatic alkane chains with more than six carbon atoms no variation in the value for DGsolv occurs. This phenomenon is reflected by the COS model results presented in Table 6.

Heat capacity, isothermal compressibility, and thermal expansion coefficient Table 2 lists for each system the two values of heat capacity Cpsim obtained from simulations at different temperatures. The empirical heat capacity corrections are also listed. From eqs. (9)–(11), three group corrections were calculated: dCp ðCH4 Þ5 2:0; dCp ðCH3 Þ57:7 and dCp(CH2) 5 10.4 J K21 mol21. The correction increases with increasing chain length as more CH2 groups are present in longer chains. The contribution of the correction term to the total heat capacity is about 30% except for methane, where it is about 4%. The final Cp values were calculated as the sum of the average over the two intermediate Cpsim values and the empirical correction (Table 7). As Cpexp for methane was used to derive the correction dCp(CH4), the final value for the heat capacity of methane is identical with the experimental value. The heat capacities are overestimated for all alkanes except butane with an average relative error of

Table 5. Comparison of the thermodynamic properties. Gibbs free energy of hydration DGhydr, Helmholtz excess free energy DFexc, and surface tension c. DGhydr [kJ mol21] System

exp

C1 C4 C6 C9 C12 C15 c-C6

8.1 9.0 10.7 12.8 15.0 17.3 5.2

45A31SPC 8.3 7.7 9.2 10.5 12.6 14.9 4.5

(0.3) (0.5) (0.7) (1.0) (1.4) (2.0) (0.7)

DFexc[kJ mol21] COS1COS/G2 8.7 9.0 10.6 11.4 16.7 16.5 6.1

(1.1) (1.9) (1.9) (2.8) (3.5) (4.1) (2.1)

c[mN m21]

exp

45A3/COS

exp

45A3/COS

0.4[a] 10.9[a] 15.7[a] 23.8[a] 32.1[a] 41.0[a] 17.3[a]

0.0 9.2 15.1 23.9 32.6 41.7 17.2

13.3 14.9 17.9 22.4 24.9 — 24.6

13.3 14.9 17.9 22.8 25.0 26.4 23.9

(0.0) (0.1) (0.3) (0.5) (0.6) (0.8) (0.0)

[a] average of the experimental data from Table 1.

796

Journal of Computational Chemistry 2014, 35, 789–801

WWW.CHEMISTRYVIEWS.COM

FULL PAPER

WWW.C-CHEM.ORG

Table 6. Comparison of the Gibbs free energy of solvation DGsolv of a single methanol molecule in the alkane liquids. The error estimates are indicated in the parentheses. DGsolv [kJ mol21] System

exp

C1 C4 C6 C9 C12 C15 c-C6

— — 26.2 25.4 — — 26.2[a]

COS/M[18]1COS

45A3145A3 29.4 26.2 25.0 24.5 24.7 24.4 25.4

(0.2) (0.4) (0.4) (0.4) (0.4) (0.6) (0.3)

210.6 27.1 26.2 25.5 25.5 25.5 26.4

(0.2) (0.4) (0.4) (0.5) (0.5) (0.6) (0.4)

[a] average of the experimental data from references[77,78]

about 6%. Only for cyclohexane, the heat capacity was strongly overestimated, by 20%. The isothermal compressibilities jT of the simulated alkanes are underestimated. The average relative error is 21%, with the largest deviation for butane: 32%. The trend of jT decreasing with increasing chain length is reproduced. The thermal expansion coefficients aT agree better with the experimental data, on average within 7.4%, with an overestimation for methane and underestimation for longer chains. The decrease with increasing chain length is reproduced.

Dielectric properties For the nonpolarizable 45A3 alkane force field, E(0) is equal to 1, because the united carbon atoms bear partial charges equal to zero. For the polarizable n-alkanes, the dependence of the polarization on the applied electric field is shown in Figure 1. The linear regression analysis of these slopes for low field strengths, that is, less than 1.0 e nm22, yields values of the static dielectric permittivities in very good agreement with the experimental ones (Table 8). The correlation coefficients of the regression analysis deviated less than 1026 from 1. If the polarizable sites with polarizability a are isotropically and homogeneously distributed with a density of qa sites per unit volume for a nonpolar liquid, that is, consisting of molecules without permanent dipole moment, the static dielectric permittivity ð0Þ can be obtained from the relation[78] ð0Þ5

3=ð4pÞ12aqa 3=ð4pÞ2aqa

(24)

derived by Mossotti in 1850. The ð0Þ values, obtained using in eq. (24) qa values derived from the simulated densities given in Table 4, and a values that are linear combinations of the polarizabilities a for the CHn united atoms given in Table 3 in proportion to the occurrence of CH2, CH3, and CH4 united atoms in the different alkanes, are also given Table 8. They indicate

Table 7. Secondary thermodynamic properties: heat capacity Cp, isothermal compressibility jT, and thermal expansion coefficient aT. Cp [J K21 mol21] System

exp

C1 C4 C6 C9 C12 C15 c-C6

55.8 134.2 194.0 283.1 376.8 469.8 153.0

jT [1025 atm21]

aT [K21]

45A3/COS

exp

45A3/COS

exp

45A3/COS

55.8 125.2 200.0 301.8 405.4 511.0 183.0

22.7 21.9 17.9 12.1 10.1 9.2 11.5

20.3 14.9 13.8 9.6 8.1 7.2 9.1

3.5 1.8 1.4 1.1 1.0 — 1.2

3.8 1.6 1.3 1.0 0.9 0.8 1.0

Figure 1. Electronic polarization Pz as a function of the applied electric field Ezext for the simulations of polarizable alkanes at 298.15 K and 1 atm. BP: simulation at boiling point temperature. Ci: alkane with i C-atoms. c-C6: cyclohexane.

Journal of Computational Chemistry 2014, 35, 789–801

797

FULL PAPER

WWW.C-CHEM.ORG

Table 8. Static dielectric permittivity (0) from experiment, calculated by linear regression of data points for EZext

A polarizable empirical force field for molecular dynamics simulation of liquid hydrocarbons.

Electronic polarizability is usually treated implicitly in molecular simulations, which may lead to imprecise or even erroneous molecular behavior in ...
1KB Sizes 0 Downloads 8 Views