Transferable Anisotropic United-Atom Force Field Based on the Mie

Subscriber access provided by University of Sussex Library

Article

Transferable Anisotropic United-Atom Force Field Based on the Mie Potential for Phase Equilibrium Calculations: n-Alkanes and n-Olefins Andrea Hemmen, and Joachim Gross J. Phys. Chem. B, Just Accepted Manuscript • DOI: 10.1021/acs.jpcb.5b01354 • Publication Date (Web): 04 Aug 2015 Downloaded from http://pubs.acs.org on August 12, 2015

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 45

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

Transferable Anisotropic United-Atom Force Field based on the Mie Potential for Phase Equilibrium Calculations: n-Alkanes and n-Olefins Andrea Hemmen and Joachim Gross∗ Institute of Thermodynamics and Thermal Process Engineering, University of Stuttgart, Pfaffenwaldring 9, 70569 Stuttgart E-mail: [email protected] Phone: +49 (0)711 685 66103. Fax: +49 (0)711 685 66140

∗

To whom correspondence should be addressed

1

ACS Paragon Plus Environment


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

Abstract A new transferable force field parameterization for n-alkanes and n-olefins is proposed in this work. A united-atom approach is taken, where hydrogen atoms are lumped with neighboring atoms to single interaction sites. A comprehensive study is conducted for alkanes, optimizing van der Waals force field parameters in 6 dimensions. A Mie n-6 potential is considered for the van der Waals interaction, where for n-alkanes we simultaneously optimize the energy parameters ǫCH3 and ǫCH2 as well as the size parameters σCH3 and σCH2 of the CH3 (sp3 ) and CH2 (sp3 ) groups. Further, the repulsive exponent n of the Mie n-6 potential is varied. Moreover, we investigate the bond length towards the terminal CH3 -group as a degree of freedom. According to the AUA (anisotropic united-atom) force field, the bond-length between the terminal CH3 -group and the neighboring interaction site, should be increased by ∆l compared with the carbon-carbon distance in order to better account for the hydrogen atoms. The parameter ∆l is considered as a degree of freedom. The intramolecular force field parameterization is taken from existing force fields. A single objective function for the optimization is defined as squared relative deviations in vapor pressure and in liquid density of propane, n-butane, n-hexane and n-octane. A similar study is also done for olefins, where the objective function includes 1-butene, 1-hexene, 1octene, cis-2-pentene, and trans-2-pentene. Molecular simulations are performed in the grand canonical ensemble with transition matrix sampling where the phase equilibrium properties are obtained with the histogram reweighting technique. The 6-dimensional optimization of strongly correlated parameters is possible, because the analytic PCSAFT equation of state is used to locally correlate simulation results. The procedure is iterative but leads to very efficient convergence. An implementation is proposed, where the converged result is not affected (disturbed) by the analytic equation of state. The resulting Transferable Anisotropic Mie-potential (TAMie) force field shows average relative deviations in vapor pressure of 1.1% and in liquid density of 0.9% for alkanes, and 2% and 1.5% for olefins, respectively, in a wide range of (reduced) temperature, from Tr = 0.55 to 0.97. For substances that were not member of the objective function,

2


Page 2 of 45

Page 3 of 45

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 TAMie force field enables predictions of phase equilibrium properties with good accuracy.

1

Introduction

The utility of molecular simulations from an engineering perspective lays in the possibility of predicting static thermodynamic properties, dynamic properties, and microscopic properties and processes with force fields as the only required input for all species. Predictions are particularly valuable for species that are not sufficiently characterized experimentally. Transferable classical force fields provide an appealing perspective for correlating and predicting such properties. Among all force field contributions, the van der Waals parameterization of effective pairwise force fields is particularly difficult to estimate from quantum mechanical calculations. A critical test for the parameterization of the van der Waals contribution to the force field are phase equilibrium properties. It was shown that force fields that accurately represent phase equilibrium properties of fluids are often also well-suited for predicting microscopic (structural) properties. 1,2 Many transferable force fields for predicting vapor-liquid phase equilibrium properties adopt a united-atom approach, where hydrogen atoms are lumped with neighboring atoms to so-called united-atom groups. A CH3 -group, say, is then represented as a single interaction site. In early studies, Siepmann et al. 3 reproduced and predicted critical properties for long chain alkanes using Gibbs ensemble technique 4,5 and configurational bias Monte Carlo 6,7 and proposed a united-atom parametrization referred to as SKS model. 8–10 One of the most popular transferable force fields is the TraPPE (transferable potentials for phase equilibria) model proposed by Siepmann and coworkers. 1,11–19 The TraPPE force field is based on intramolecular and electrostatic models of OPLS-UA (optimized potentials for liquid simulations) proposed by Jorgensen 2 and uses a Lennard-Jones approach for the van der Waals

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

interactions. The Lennard-Jones parameters were fitted to liquid density and to critical temperatures of pure fluids and TraPPE thus accurately describes vapor-liquid coexistence in a temperature-density projection but systematically overestimates vapor density and vapor pressure. 1,12 Jorgensen et al. 2,20–23 proposed the OPLS-UA parametrization for simple organic compounds up to proteins. 24 In the OPLS-UA force field for hydrocarbons, 2 the Lennard-Jones parameters are optimized to describe liquid density and enthalpy of vaporization, so that OPLS-UA achieves good agreement to experimental data for small components (C2 − C6 ). Nath et al. 25–28 proposed the NERD force field for different organic compounds. The Lennard-Jones parameters for hydrocarbons and olefins in the NERD force field 25,26 were adjusted to orthobaric densities and second virial coefficients of n-alkanes and to orthobaric densities of olefins and simulation results show a good agreement with experimental coexistence properties for short and intermediate-length components. Nath et al. 25 raise the question of whether united-atom force fields adjusted to vapor-liquid data can also describe transport properties with reasonable accuracy (and they point to AUA-type models, 29 that show promising results). Toxvaerd 30 argued that a united-atom site that represents CH3 , for example, should not necessarily be placed at the position of the carbon atom. In order to better account for the hydrogen atoms, he defined spacial offsets for united-atom interaction sites (where hydrogen is part of the united-atom group). A modified and extended parameterization lead Ungerer and Boutin et al. 31–44 to the AUA4 (anisotropic united-atom) force field for a variety of organic compounds. The AUA4 force field for chain alkanes 31 and olefins 35 was optimized on the basis of equilibrium properties (vapor pressure, enthalpy of vaporization and liquid density). Subsequently, the torsional potential was adjusted to reproduce different transport properties without critically sacrificing the equilibrium properties of the model. 38 All of the above mentioned force fields describe intermolecular interactions with the Lennard-Jones 12-6 potential, whereby in each force field (TraPPE, OPLS-UA, NERD, AUA4) the diameter parameter and energy parameter were optimized to reproduce the

4


Page 4 of 45

Page 5 of 45

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


desired thermodynamic properties. The transferability of the united-atom parameters for each force field was analyzed and confirmed. The NERD force field thereby assumes the transferability only for alkanes with four or more carbon atoms, whereas the CH3 -group was individualized for propane and propene. Errington and Panagiotopoulos realized that it is not possible to correlate coexisting density data and vapor pressure data simultaneously with good accuracy using a LennardJones 12-6 potential. 45 They obtained very good results using a Buckingham exponential-6 potential, where the repulsive part of the potential is different to the Lennard-Jones potential. With a similar rational as Toxvaerd, 30 they shifted the CH3 -group outwards in order to account for the hydrogen atoms, thereby improving agreement with experimental data of ethane. Potoff and Bernard-Brunel considered Mie n − 6 potentials to describe intermolecular interactions and showed that increasing repulsive exponents n and simultaneously elevating the value of energy parameter ǫ favorably shifts the vapor pressure curve for a given saturation density ρ(T coex ). 46 Force field parameters were adjusted to liquid density and vapor pressure data and very good agreement was found for the simulation results with experimental data for short, intermediate-length and long n-alkanes. Adjusting transferable force field parameters to vapor-liquid equilibria is somewhat tedious, when executed in a simultaneous optimization problem. The problem is multidimensional and is numerically noisy, because of the statistical uncertainties in simulation results. In practice, most studies have avoided the full problem and decomposed the iteration into sets of fairly uncorrelated subproblems. For optimizing the Lennard-Jones parameters of a united-atom force field to n-alkanes (i.e. CH3 and CH2 ) the problem is 4 dimensional. Most studies have parameterized the two parameters of the CH3 -group independently, e.g. by considering only ethane for determining the two Lennard-Jones parameters of the CH3 group. 1,2,46 More recent studies apply advanced numerical methods to simultaneously optimize force field parameters: Faller et al. 47 applied the Nelder-Mead simplex algorithm to optimize a

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

variety of atomistic force fields, Ungerer and coworkers 31,35 used a direct gradient approach for only small molecules and H¨ ulsmann et al. 48 developed a gradient based optimization workflow for the automated development of molecular models. In a previous study, van Westen et al. showed that analytic equations of state can be utilized for optimizing the force field parameters. 49 It was shown that the PC-SAFT equation of state 50 can be used to correlate molecular simulation data using the force field parameters of the simulations. The deviation between calculated and experimental values of some defined thermodynamic observable can then be minimized by varying the force field parameters within the PC-SAFT model. Multidimensional parameter optimizations can be conducted in milliseconds using the analytic PC-SAFT equation of state. Since the PCSAFT model is not exact in representing the molecular simulations, the approach has to be repeated iteratively. The appeal of this method is that the iterative procedure converges in a few steps only (typically 3 to 5 steps), even for multidimensional problems. This study comprehensively analyzes and parametrizes a force field for n-alkanes and n-olefins with anisotropic Mie n − 6 potentials. The position of the CH3 -group, the repulsive exponent n of the Mie potential and the energy and size parameters of the Mie potential are considered as degrees of freedom. The objective function for the optimization problem considers vapor pressure and liquid densities. The proposed force field is referred to as Transferable Anisotropic Mie (TAMie) force field and it is analyzed for static thermodynamic properties. This paper is organized as follows: the molecular model together with intramolecular and intermolecular potential functions and force field parameters are presented first, followed by a description of the simulation method and simulation details. The next chapter briefly explains the optimization algorithm and the last chapter presents and discusses the simulation results of the TAMie force field in comparison with experimental data and with other force fields.

6


Page 6 of 45

Page 7 of 45

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


2

Molecular Model

A united-atom approach is used in which carbon atoms and their adjacent hydrogen atoms are defined as a single interaction site. The center of interaction of the methyl end group CH3 (sp3 ) is moved outward to better account for the effect of the bonded hydrogen atoms. The interaction site is shifted by ∆l in the direction of the carbon-carbon bond axis. The parameter ∆l is here considered as a degree of freedom that is varied to ensure good agreement with the physical properties of real substances. This molecular approach stems from Toxvaerd’s anisotropic united-atom model, 30 where the center of interaction (center of force) for both, CH3 and CH2 groups, are moved outward to better reflect the presence of the bonded hydrogen atoms. Moving the methylene CH2 interaction site from the position of the carbon, however, leads to a more cumbersome configurational bias problem for inserting and deleting molecules. The additional complexity is introduced because the angle of the carbon to carbon interaction site, is only defined, once the following interaction site is placed. Although an efficient implementation of a configurational bias method has been proposed by Smit et al., 8 we choose for simply placing the methylene interaction site on the center of the carbon atom. Van der Waals (pair) interactions are described with the Mie potential, a generalized Lennard-Jones (n-m) potential, as

u(rij ) = cn ǫij

σij rij

n

−

σij rij

m

(1)

with prefactor cn that ensures a minimum value of −ǫij for the attractive well, as cn =

n n−m

m n n−m m

(2)

where rij denotes the distance between two interaction sites i and j, and ǫij and σij are the the well depth and diameter parameter, respectively.

7



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 45

The term with power m represents the attractive part of the van der Waals interactions due to dynamically induced dipole-dipole interactions. The dispersive exponent is here fixed to m = 6. This value has some physical significance, since the leading term of the dispersion energy expansion is of power m = 6. 51 The repulsive exponent n is in our work varied from 12 to 16 to study its influence on vapor liquid equilibria. The intermolecular potential u(rij ) governs the interactions between united-atom groups that belong to two different molecules and between united-atom groups within one molecule that are separated by more than two bonds. Parameters for interactions between unlike types of interaction sites are determined using Lorentz-Berthelot combining rules, 52,53 as

σij = (σii + σjj ) /2

(3)

√ ǫij = ǫii ǫjj

(4)

Prefactor cn in Eq. (2) can also be defined for groups of cross-wise different powers nii and njj . Potoff and Bernard-Brunel 46 propose an arithmetic combination rule nij = 0.5(nii +njj ) to be considered in Eq. (2). For the intramolecular force field, we adopt a parameterization available in literature. 1,2,8,12,54,55 We consider fixed bond lengths between the interaction sites, since bond vibrations are of high frequency and of low amplitude and their effect is unimportant for many fluid properties. 56 Bending angles between united-atom groups are generated according to the harmonic bending potential 54 ubend (θ) = kθ /2 (θ − θ0 )2

(5)

with kθ , θ, θ0 as force constant, bending angle and zero-Kelvin angle, respectively. For the torsional potential between four neighboring interaction sites of n-alkanes and of α-olefins

8


Page 9 of 45

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


we use the OPLS-UA model, 2 according to

utor (φ) = c0 + c1 [1 + cos(φ)] + c2 [1 − cos(2φ)] + c3 [1 + cos(3φ)]

(6)

whereas, the torsion around a double bond in cis-/trans configurations of β-olefins is described through a harmonic potential, and thus becomes a ‘trapped’ degree of freedom, 12 with utor (φ) = d0 /2 (φ − φ0 )2

(7)

where d0 is a force constant and φ0 denotes the zero-Kelvin angle. Both parameters were obtained by adjusting the energies to fully relaxed structures. 12 In order to comprehensively present the proposed force field, we summarize all model constants in Table 1 along with the corresponding literature references. The remaining parameters of the TAMie potential are summarized in Table 2. Furthermore, optimized 12-6 Lennard-Jones parameters are proposed, that where optimized to reproduce saturated vapor pressure and liquid densities of vapor-liquid equilibria using intramolecular potentials and constants as listed above.

3

Molecular Simulation Technique

This chapter firstly reviews the Monte Carlo simulation technique for phase equilibria of pure components and binary mixtures and, secondly, summarizes all relevant simulation details.

3.1

Simulation Method

This study started using the grand canonical Monte Carlo (GCMC) simulation technique as described by Castillo-Sanchez et al. 57 The chemical potential µ, needed as a specification (in addition to the volume V and temperature T ) is estimated from the PC-SAFT equation of state. We subdivide the N -space into windows of width ∆Nk = 5 or ∆Nk = 10, where 9



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 10 of 45

Table 1: Parameters of intramolecular potential for TAMie force field, taken from literature

bond length 1,12 r0 / ˚ A

bond type CH2 − CH2 CH3 − CH2 CH(sp2 ) = CHz (sp2 ) z ∈ [1, 2] CH3 − CH(sp2 ) CH3 − CH3 ethane

1.54 1.54 + 0.2 1.33 1.54 + 0.2 1.54 + 2·0.2 bond angle 8,12,54,55

bond type CHx − CH2 − CHy x, y ∈ [2, 3] CHz = CH − CHy z ∈ [1, 2], y ∈ [2, 3]

Θ0 /deg

kΘ /kB /K

114.0 119.7

62500 70420

torsional potential 1,2,12 torsion type CHx − CH2 − CH2 − CHy x, y ∈ [2, 3] CHz = CH − CH2 − CHy z ∈ [1, 2], y ∈ [2, 3]

c0 /kB /K

c1 /kB /K

c2 /kB /K

c3 /kB /K

0 688.5

355.03 86.36

-68.19 -109.77

791.32 -282.24

d0 /kB /K

φ /rad

12400 13400

π 0

CHx − CH = CH − CHy (cis) x, y ∈ [2, 3] CHx − CH = CH − CHy (trans) x, y ∈ [2, 3]

k is an index counting the windows. 58 All windows run in parallel, sampling the probability distribution Pk (N ) using a transition-matrix scheme. 59,60 The advantage of dividing the problem into small windows is twofold: first, the calculation time for each window is low and all windows can run in parallel. Secondly, suitable simulation parameters, such as maximum displacement or number of configurational bias steps are defined for a narrow range of molecule number N . The phase equilibrium properties for the given temperature are determined a posteriori using a histogram reweighing scheme. 61,62 The simulation method leads to accurate phase equilibrium properties, but since a single temperature is considered in all windows, the histogram reweighing can not cover wide ranges of temperatures away 10


Page 11 of 45

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


Table 2: Parameters for van der Waals interactions: energy parameter ǫ, size parameter σ, and repulsive exponent n the proposed of TAMie model and for an optimized 12-6 Lennard-Jones potential 12-6 LJ

TAMie

pseudoatom

ǫ/kB /K

σ /˚ A

n

ǫ/kB /K

σ /˚ A

n

CH3 (sp3 ) CH2 (sp3 ) CH2 (sp2 ) olefins CH(sp2 ) olefins CH3 (sp3 ) ethane

103.215 45.44

3.746 12 4.051 12

136.318 52.9133 100.681 53.9515 130.780

3.6034 4.0400 3.6005 3.8234 3.6463

14 14 14 14 14

from the specified temperature. For every component we therefore performed simulations at 5 temperatures. A modification of the calculation scheme was proposed to us by Panagiotopoulos, where various temperatures are considered in the different windows, thereby allowing the histogram reweighing scheme to cover the entire fluid phase envelope. We have adopted the algorithm in a later stage of this study and present a more detailed analysis of the method in a separate publication. 63 A single mixture is considered in this work, in order to once test the TAMie force field for mixtures. We use a transition matrix Monte Carlo 64 (TMMC) scheme in a multi-canonical ensemble (N1 , N2 , V, T ) for a binary mixture of component 1 and 2. We divide the N1 N2 -surface into ∆N1 × ∆N2 windows of 11 × 11 molecules. Exactly 1/121st part of the MC-production steps are spent for every N1 , N2 -combination, sampling the transition probabilities for insertion- and deletion-trials. Histograms of the particle number distribution are constructed by combining the transition matrix probabilities for each window. The phase equilibrium is a posteriori determined from histogram reweighting method. A more detailed description of the method will be presented in a subsequent study on mixtures. 63

11



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

3.2

Simulation Details

The grand canonical Monte Carlo simulations were conducted with box lengths L of cubic simulation boxes chosen as 39 ˚ A for (C3 -C6 ), 43 ˚ A for (C8 ), up to 61 ˚ A for (C24 ). The Mie A with analytipotential is calculated for all interaction pairs up to a radial cutoff of 14 ˚ cal long range corrections 56,65 where a radial distribution function of g(r) = 1 is assumed. The chemical potential µk for a given temperature Tk is calculated with PC-SAFT. 57 The N -coordinate, that is the coordinate along which the number of molecules N fluctuates, is subdivided into windows of with ∆Nk = 10 for a window k at moderate densities and ∆Nk = 5 for higher densities. Trial moves that would lead to N outside of the window are trivially rejected. For windows of width ∆Nk = 10, we used 8 M (million) MC steps for equilibration and 38 M production steps, whereas for windows of ∆Nk = 5 we defined 12 M MC steps for equilibration and 27 M steps for the production. Five different Monte Carlo moves are performed including particle translational displacement 27.5 %, rotation 27.5 %, particle insertion 15 %, particle deletion 15 %, and particle regrow 15 % with their appropriate probabilities, respectively. The bias potential ω(N ) determined from transition matrix sampling in the course of the simulation, is updated every 500k MC steps and histogram data Hk (N, E) is collected for every step.

4

Optimization Algorithm

The optimization algorithm described by van Westen 49 (and by Schacht et al. for mixtures 66 ) is modified and is applied to optimize the intermolecular force field parameters. The PC-SAFT equation of state 50,67,68 is used to correlate and subsequently predict results of molecular simulations, allowing for a a swift convergence of multidimensional force field parameters. The modification of the algorithm introduced here is essential because it ensures the analytic equation of state does not affect the objective function in the minimum. Once convergence is achieved, the equation of state does not act on the objective function

12


Page 12 of 45

Page 13 of 45

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


(or the gradients of the objective function) as we show in the Appendix. Analytic equations of state have successfully been applied for force field development. Stoll et al. developed accurate analytic expressions for coexisting densities and pressure as a function of temperature for a very comprehensive set of molecular simulations of 2center Lennard-Jones fluids with an embedded point dipole 69 or point quadrupole. 70 Using these expressions, Vrabec et al. 71 and Stoll et al. 72 subsequently identified reliable force field parameters for a large set of real substances. M¨ uller, Jackson, Adjiman and Galindo together with colleagues used the SAFT-γ model 73–75 for predicting force field parameters with very good results. 76–78 The underlying molecular model of SAFT-γ (and for PC-SAFT) is coarse-grained compared with the force field considered in molecular simulations. The molecules are represented as chains of spherical segments where bond-angles, for example, are not constraint. The resulting force field is thus coarse-grained, but is obtained without the need to perform iterative molecular simulations. Our approach is different, because we target at classical atomistic force fields (on a unitedatom level), whereas SAFT models do not resolve the same molecular detail. It is possible, nonetheless, to support in the parameter optimization for such force fields. The approach is inevitably iterative, but leads to very efficient convergence of multidimensional van der Waals parameters. The idea is to first conduct molecular simulations for an initial guess of force field parameters. We then approximately relate parameters of the force field to parameters of the PC-SAFT equation of state. The result of molecular simulations are subsequently correlated with the equation of state. The so-obtained analytic equation of state is used to minimize an objective function (i.e. deviations to experimental data of a set of substances). The resulting force field parameters serve as a next guess for conducting molecular simulations. The force field iteration of highly correlated van der Waals parameters converges in a few steps. A graphical representation of the algorithm is given in the Appendix. The force field parameters used in molecular simulations can be approximately related to PC-SAFT parameters. 49 Details of how the interaction sites are arranged according to the

13



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 14 of 45

force field (defined in terms of bond-lengths and bond angles) are thereby coarse grained using a parameter mi . Parameter mi , for each united-atom group, i describes the contribution of a united-atom group (force field) to the effective chain length m ˆ in PC-SAFT. This parameter mi is taken from a group contribution approach. 79 The relation between the van der Waals force field parameters {σij , ǫij } and the PCSAFT parameters {m, ˆ σ ˆ , ǫˆ} is given in three equations. These equations define the three pure component parameters of the PC-SAFT model, i.e. segment number m, ˆ segment diameter parameter σ ˆ and segment energy parameter ǫˆ, according to united-atom types

m ˆ =

X

Ni m i

(8)

i

where Ni is the number of united-atom groups of type i. With the assumption that a molecule holds the same ’volume’ in both frameworks (force field approach and PC-SAFT) the segment diameter parameter σ ˆ of PC-SAFT is related to the force field segment diameter parameters σii of all united-atom groups i, as united-atom types

m ˆ (ˆ σ φσ ) 3 =

X

Ni mi σii3

(9)

i

The relation between the energy parameter ǫˆ of PC-SAFT and the van der Waals force field parameter ǫ can be derived from a first order perturbation theory approach, which leads to u.-a. u.-a.

m ˆ 2 (ˆ σ φσ )3 (ˆǫφǫ ) =

types types X X i

Ni Nj σij3 ǫij

(10)

j

The two parameters φσ and φǫ in Eq. (9) and Eq. (10), respectively, are introduced to correct the approximations inherent in these relations. Both parameters, φσ and φǫ , are close to unity for Lennard-Jones potentials. 49 It was shown, that the optimization is not sensitive towards the chosen group-contribution parameters m ˆ i of the PC-SAFT model.

14


Page 15 of 45

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


With the three relations, Eq. (8-10), we have expressed the PC-SAFT model in terms of the van-der-Waals parameters of the force field. For a compact notation, we introduce p as the vector containing the force field parameters σii and ǫii of all adjustable groups i. The procedure for optimizing the force field parameters is iterative and requires starting values for all van-der-Waals force field parameters. Molecular simulations are conducted with the starting values of p and the result is given as a set of exp observables, Ωsim . Now, the j (p), where index j counts the observed data points, j = 1, ...N

two parameters φσ,k and φǫ,k of the PC-SAFT model according to Eq. (9) and (10) (for every substance k considered in the objective function) are adjusted, so that the PC-SAFT model 49 very closely represents the properties Ωsim we j (p). As a modification to our previous work dev now capture all deviations between Ωsim j (p) and PC-SAFT model in a vector Ωj . Next,

the PC-SAFT model minimizes the objective function exp

N exp 2 1 X Ωsim j (p) − Ωj f (p) = exp N Ωexp j j=1

(11)

using the force field parameters p as degrees of freedom. The calculated values of Ωsim j (p) are thereby corrected by the vector Ωdev to account for inaccuracies of the PC-SAFT model. j The resulting vector p is the input for the next iteration, i.e. for conducting new molecular simulations. The here proposed modification of correcting the calculated results of PC-SAFT by small offsets Ωdev may appear as a subtlety. But in the appen dix we show that the final j result of our algorithm, through this modification, is unaffected by the analytic equation of state. While approximations of the equation of state do not affect the objective function in the minimum, as the appendix clarifies, we note that model deficiencies can very well weaken the convergence behavior. Because we are going to apply the PC-SAFT equations Eq. (8-10) to force fields of Mie n-6 potentials, for which the PC-SAFT model is not explicitly suited, it is worthwhile to describe some observations. Applying the algorithm to fluids with, say, a Mie

15



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

16-6 potential leads to values of parameter φǫ in Eq. (10) quite different to unity. Nonetheless the convergence of the algorithm did not deteriorate. Eq. (9) and (10) provide some rational for this robustness. Eq. (10), for example, is derived by applying a perturbation approach to a fused sphere model representing the force field and a tangent sphere model underlying the PC-SAFT model. Both perturbation approaches lead to correlation integrals which were cancelled out of Eq. (10) approximately. 49 The absolute value of one of these correlation integrals changes considerably with the repulsive Mie n-parameter, but the dependence on density and temperature changes only moderately. The change of absolute value is captured by the parameter φǫ . The manner of how the van der Waals force field parameters {σij , ǫij } are correlated, however, is weakly dependent on the repulsive Mie n-parameter. In this work, we define the objective function Eq. (11) to include vapor pressure data and liquid density data of several species as observables Ωj (all with equal weight). The reduced temperature ranges from Tr = 0.55 to 0.97. Force field parameters of different unitedatom groups within one molecule are highly correlated. To obtain optimized transferable force field parameters we therefore consider experimental vapor-liquid equilibrium properties (liquid density and saturation pressure) of different pure components simultaneously. For the homologous series of n-alkanes, for example, the species considered in the objective function are propane, n-butane, n-hexane and n-octane.

16


Page 16 of 45

Page 17 of 45

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


5

Results and Discussion

5.1

n-alkanes

The homologous series of n-alkanes is important for transferable force fields, because the pair potential of the CH2 and the CH3 groups are not superpositioned by significant static polar contributions. In this study, we investigate the force field parameters of n-alkanes in 6 dimensions, i.e. the four Mie-potential parameters {σCH3 , ǫCH3 , σCH2 , ǫCH2 }, as well as the repulsive exponent of the Mie potential, {n}, and the outward shift of the CH3 -group, {∆l}. We thereby simultaneously consider the properties of 4 n-alkanes, namely propane, n-butane, n-hexane and n-octane. For optimizing the Mie parameters, σi and ǫi , we use the algorithm described in section 4, that allows simultaneous optimization in a simple way. The iteration progress for the 4-dimensional optimization of the parameters {σCH3 , ǫCH3 , σCH2 , ǫCH2 } is for a representative case (n = 16, ∆l = 0.2 ˚ A) shown in Fig. 1. The 4 dimensional iteration converges quickly and it can be terminated after 4 to 5 steps, usually. While the objective function for the optimization is defined as relative mean squared deviations, according to Eq. (11), in Fig. 1 and in the entire chapter, we report average absolute deviations (AAD), PN exp Ωsim (p)−Ωexp j j as AAD% = N100 . exp j=1 Ωexp j

The parameters {σCH3 , ǫCH3 , σCH2 , ǫCH2 } are strongly correlated and the swift conver-

gence (Fig. 1) is due to the proposed optimization algorithm. In the optimization, we found no indication of local minima. In preparatory studies we chose various starting values for the force field parameters and confirmed convergence to the same minimum. For parameter identification studies using analytic equations of state (simultaneously adjusting parameters to a collection of substances) it is known that local minima exist. 80 Local minima emerge for cases where, for example, a calculated critical temperature falls below an experimental vapor pressure point, prohibiting this experimental vapor pressure point to be considered in the objective function. The objective function is then non-smooth. The appearance of local minima for a certain range of parameters (and for merely Mie-based force fields) is thus

17



a result of how the objective function is defined. From all our observations, we have the conjecture that our objective function if free of local minima in the range of parameters we investigated. We attribute this robustness to the fact that temperatures above 0.95 times the experimental critical temperature are omitted in our objective function.

20

15

AAD / %

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

5

0 0

1

2

3

4

number of iterations

Figure 1: Example for the typical iteration progress for the 4-dimensional optimization of the parameters {σCH3 , ǫCH3 , σCH2 , ǫCH2 }, for a given repulsive exponent n and outward shift ∆l of the CH3 -group. Here, n = 16 and ∆l = 0.2˚ A and the data are for propane, n-butane, n-hexane and n-octane. We so far discussed the parameters {σCH3 , ǫCH3 , σCH2 , ǫCH2 }. The remaining 2 dimensions, namely the repulsive exponent of the Mie potential, n, and the outward shift of the CH3 -group, ∆l, were investigated by parameter variation. Fig. 2 shows the resulting objective function f as a function of n and ∆l. Note, that each point of the surface is optimized in the 4 dimensions {σCH3 , ǫCH3 , σCH2 , ǫCH2 }, so that a minimum value of f in this figure corresponds to a minimum in the 6 dimensions we here investigate. Table 3 shows the absolute averaged deviation for liquid densities, vapor pressure corresponding to Fig. 2. It is interesting to see, that the Lennard-Jones potential leads to a comparably large value of the objective function of about 4.0% with our definition of the objective function. The resulting Lennard-Jones parameters are significantly different to those of the TraPPE force field (see Table 2). Since our objective function does not stress the T -ρ-projection, we get much better results for the vapor pressure curve but give in on the coexisting densities and in particular our parameterization strongly sacrifices the critical point. It has previously been reported, that it is not possible with the Lennard-Jones potential and ∆l = 0 ˚ A to 18


Page 18 of 45

Page 19 of 45

simultaneously represent all of the above properties with good accuracy. The value of the objective function of about 4.0 AAD% reflects this problem.

4.5 4 3.5

AAD / %

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


3 2.5 2 0.2

1.5 1 12

0.1 13

14

15

16

0 ∆l/˚ A

n

Figure 2: Objective function f as a function of n and ∆l with each point of the surface optimized in the 4 dimensions {σCH3 , ǫCH3 , σCH2 , ǫCH2 }. Fig. 2 suggests that the problem of simultaneously representing the entire p-ρ-T surface of the n-alkane family can be alleviated by increasing the repulsive exponent n, or by increasing ∆l. It is interesting to isolate how the repulsive exponent n acts on the vapor pressure curve. Fig. 3 shows the relative deviations of the calculated vapor pressure from the experimental value for n-hexane. The results for n-hexane are quite representative for all of the n-alkanes of our objective function. The slope of the vapor pressure curve changes strongly when varying the exponent from 12 to 16. Also the liquid and vapor densities and the critical point are strongly improved compared to the Lennard-Jones parameterization. The critical point properties reported in this study are not extrapolated to infinite sized systems. For n-hexane, as an example, Table 3 provides both values , the apparent critical point, which depends on the size of the simulation unit-cell, and a critical point corresponding to an infinite sized system.

Increasing the outward shift of the CH3 group, ∆l, also has a pronounced effect on the description of the n-alkanes (as Fig. 2 shows). Fig. 4 illustrates, how the deviations in the vapor pressure of n-hexane change, with varying ∆l. We note, that although only the results

19



Rel. dev. /% vapor pressure

10 n = 12 n = 13 n = 14 n = 15 n = 16

5

0

-5

-10

300

400

350

450

500

T/K

Figure 3: Relative deviations of calculated vapor pressure from the experimental value for n-hexane for repulsive exponent n ∈ [12, 14, 15, 16] and ∆l = 0˚ A. Parameters {σCH3 , ǫCH3 , σCH2 , ǫCH2 } are in optimized state. for n-hexane are shown in Fig. 3 and Fig. 4, the parameter regression was done for 4 alkanes simultaneously in all cases. A summary of errors in vapor pressure, liquid densities and in the apparent critical temperature is given in Table 3. All the simulated data for varying n and ∆l, and optimized in {σCH3 , ǫCH3 , σCH2 , ǫCH2 }, is provided as a supporting information. 10

Rel. dev. /% vapor pressure

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 20 of 45

∆l=0 ∆ l = 0.1 ∆ l = 0.2

5

0

-5

-10

300

400

350

450

500

T/K

Figure 4: Relative deviations of calculated vapor pressure from the experimental value for n-hexane for outward shift of CH3 with ∆l ∈ [0, 0.1, 0.2]˚ A and n = 14. Parameters {σCH3 , ǫCH3 , σCH2 , ǫCH2 } are in optimized state. According to Fig. 2 and Table 3 we find equally good regression results for n = 15 and ∆l = 0.2 ˚ A and for n = 14 and ∆l = 0.2 ˚ A, when limiting consideration to integer values of n. Fig. 2 indicates that a higher value of ∆l may lead to a somewhat lower error for n = 14. A maximum value of ∆l = 0.2 ˚ A was here considered in order to not overdraw the physical 20


Page 21 of 45

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


significance of the ∆l-parameter, as a parameter that implicitly accounts for the hydrogen interactions in a CH3 -group. We have a slight preference for even numbers of the repulsive exponent n, because that ensures integer valued cross exponents nij when considering mixed force field exponents with, say, a Lennard-Jones (njj = 12) potential. Considering all of the above, we propose a force field with n = 14 and ∆l = 0.2 ˚ A for the n-alkane family. Fig. 5 - 10 present coexisting densities and saturated vapor pressures of selected n-alkanes as obtained from the TAMie force field in comparison to experimental data. Experimental data of the smaller n-alkanes are actually correlation results from accurate multi-parameter equations of state as provided by the NIST Chemistry webbook. 81 Quasi-experimental data for hexadecane, eicosane, and tetracontane are results from multi-parameter correlations provided by DIPPR data base 82 and quasi-experimental values for the vapor density are estimated from the Peng-Robinson eqauation of state. 83 Another recent and very accurate force field proposed by Potoff and Bernard-Brunel (here referred to as PBB) 46 is also included in these Figures. The diagrams confirm a good overall agreement between TAMie force field and experimental data. The critical point is estimated at too high temperatures towards tetracosane. While the force field of Potoff and Bernard-Brunel gives a better description of the critical point, we see a weaker agreement to the vapor pressure data when compared to the TAMie force field. This is a consequence of the selected objective function. The simulation results for all n-alkanes are obtained for reduced temperatures Tr = 0.55 to 0.97 and deviations are calculated as absolute average deviations with respect to experimental data. Let’s first consider n-alkanes that were member of the objective function (propane, n-butane, n-hexane, n-octane). The average deviation in liquid density for the TAMie model is 0.89 % and for PBB it is 0.55 %. For both force fields the deviations of liquid densities to experimental values are small and nearly constant for a wide temperature range (see Fig. 5, 6 and 8, 9). Deviations in vapor pressure are 1.08 % for the TAMie force field, and 3.79 % for the PBB parameterization. We mention that for these numbers we have performed own simulations for the PBB force field because we achieved better results for

21



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 PBB model (at low temperatures) than the results published by the developers 46 of the force field. We note further, that our parameterization for n = 16 and ∆l = 0 ˚ A is slightly different to the PBB parameterization due to a different objective function. For n-pentane, n-decane, n-dodecane as substances that were not part of the objective function, the TAMie force field leads to errors of 0.95 % and 1.98 % for liquid density and vapor pressure, respectively. Fig. 5 - 10 give a graphical representation of these results. Deviations of the PBB force field are 1.15 % for liquid density and 5.95 % for vapor pressure. While TAMie does not represent the critical point as accurately as PBB, vapor pressure and vapor densities are very well represented. A comparison of the TAMie model to other force fields, i.e. TraPPE, 1 NERD, 25 the model of Spyriouni et al., 84 and PBB 46 is given in the supporting information. We have chosen n-pentane and n-decane for this comparison, because these substances were not part of our objective function. For n-hexadecane, eicosane, and n-tetracosane the results of TAMie do not significantly degenerate, as Fig. 7 and 10 indicate. Fig. 5 and 8 also shows coexisting densities and vapor pressures of ethane. The parameters of the Mie potential for ethane were optimized for n = 12, 14, and 16 and for outward shifts of the CH3 -group of ∆l = 0 ˚ A, 0.1 ˚ A, and 0.2 ˚ A. The best result was also found for n = 14 and 0.2 ˚ A with parameter values listed in Table 2. All results are summarized in the supporting information. The errors of the TAMie model are 0.61 % and 0.80 % for vapor pressure and liquid density over the considered temperature range. Methane is not included in Fig. 5 and 8 since our optimization procedure essentially confirmed a parameterization proposed by Potoff and Bernard-Brunel. 46 In order to assess whether a Mie 14-6 potential is compatible with a Lennard-Jones 12-6 potential, we consider the binary mixture of n-butane with methane. We apply the TAMie (14-6) force field for butane and adopt a Lennard-Jones parameterization proposed by Vrabec et al. 85 for methane (σCH4 = 3.7281 ˚ A, ǫCH4 = 148.55 K, n = 12). We consider the combining rules, Eq. (3) and (4) without any correction for cross interactions between

22


Page 22 of 45

Page 23 of 45

the two components. For calculating the vapor-liquid phase envelope we apply histogram reweighting. The mixture critical points are determined with the Binder cumulant intersection method where analysis of the finite size scaling is done for three different simulation volumes. 86 Fig. 11 gives the simulation results for the methane-butane mixture at two temperatures, T = 344.26 K and 394.26 K. The comparison to the experimental data confirms a very good agreement for this mixture. Although a more comprehensive investigation is needed, the result is promising because it suggests that mixtures of compounds with unlike repulsive exponents n of the Mie n-6 potential can be described with simple combining rules. A very good description for mixtures of substances with different Mie-n parameters was also found by Potoff and Bernard-Brunel. 46 These authors introduced a correction to the combining rule (Eq. (4)), but this was done for a more non-ideal mixture.

hexane

500

pentane butane

400

T/K

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


propane

300

200 0

ethane exp. data [NIST] Potoff, Bernard-Brunel TAMie

100

200

300

400

500

600

700

3

l / kg/m

Figure 5: Vapor-liquid coexistence curve for short n-alkanes. Symbols are simulation results for TAMie force field (triangles) and Mie potential from Potoff and Bernard-Brunel (circles). 46 Solid lines and crosses correspond to experimental data. 81

23



Table 3: Summary of average absolute deviations (AAD) in vapor pressure, liquid densities and apparent (system-size dependent) critical temperature of n-alkanes (propane, n-butane, n-hexane and n-octane) for repulsive exponent n ∈ [12, 14, 15, 16] and outward shift of the CH3 -group ∆l ∈ [0, 0.1, 0.2]. AAD-% n

∆l /˚ A

psat

ρL

Tc

12 12 12

0 0.1 0.2

5.68 4.97 3.79

2.38 2.02 1.70

3.40 3.06 2.48

13 13 13

0 0.1 0.2

4.46 3.50 2.94

1.79 1.78 1.54

2.38 2.44 2.15

14 14 14

0 0.1 0.2

1.87 1.21 1.07

1.47 1.23 0.90

1.90 1.66 1.12

15 15 15

0 0.1 0.2

1.12 1.34 1.62

0.99 0.97 0.73

1.23 1.19 0.80

16 16 16

0 0.1 0.2

1.79 2.11 3.33

0.85 0.75 0.75

0.72 0.52 0.70

700 dodecane decane

600

octane

T/K

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

500

400

300 0

experimental data Potoff, Bernard-Brunel TAMie

200

400

600 3

l / kg/m

Figure 6: Vapor-liquid coexistence curve for intermediate n-alkanes. Symbols are simulation results for TAMie force field (triangles) and Mie potential from Potoff and Bernard-Brunel (circles). 46 Solid lines and crosses correspond to experimental data. 81

24


Page 24 of 45

Page 25 of 45

900 tetracosane eicosane

800

hexadecane

T/K

700 600 TAMie experimental data

500 400 0

100

200

300

400

600

500

700

3

l / kg/m

Figure 7: Vapor-liquid coexistence curve for longer n-alkanes. Symbols are simulation results for TAMie force field (triangles), solid lines correspond to experimental data from. 81,82

16

Pressure ln(p/Pa)

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


experimental data TAMie Potoff, Bernard-Brunel

14

12

10

ethane

8

hexane

2

3

butane pentane

4

propane

5

6

Temperature / 1000/K

Figure 8: Vapor pressure of short n-alkanes. Symbols are simulation results for TAMie force field (triangles) and Mie potential from Potoff and Bernard-Brunel (circles). 46 Solid lines correspond to experimental data. 81

25



16 TAMie Potoff, Bernard-Brunel experimental data

ln(p/Pa)

14 12 10 8

dodecane decane

6

1.5

2

octane

3

2.5

3.5

-1

(T / 1000 K)

Figure 9: Vapor pressure of intermediate n-alkanes. Symbols are simulation results for TAMie force field (triangles) and Mie potential from Potoff and Bernard-Brunel (circles). 46 Solid lines correspond to experimental data. 81

experimental data TAMie

14 12

ln(p/Pa)

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 26 of 45

10 8 tetracosane eicosane

6 1

hexadecane

2

1.5

2.5 -1

(T / 1000 K)

Figure 10: Vapor pressure of longer n-alkanes. Symbols are simulation results for TAMie force field (triangles) and Mie potential from Potoff and Bernard-Brunel (circles). 46 Solid lines correspond to experimental data. 82

26


Page 27 of 45

experimental data TAMie

100

p/ bar

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


50

0 0

0.2

0.4

xmethane

0.6

0.8

Figure 11: Vapor-liquid equilibrium of n-butane-methane mixtures at two temperatures, T = 344.26 K and 394.26 K. Comparison of molecular simulation results using the TAMie model (triangles) to experimental data 87 (crosses).

27



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

5.2

n-olefins

For n-olefins, we distinguish between the terminal united atom group CH2 (sp2 ) of α-olefins and between the united-atom group CH(sp2 ). A β-olefin, such as 2-butene is composed of two CH3 (sp3 ) groups and two CH(sp2 ) groups. Cis and the trans isomers of β-olefins are only distinguished through the torsional potential, Eq. (7). For the n-olefins we perform a 5 dimensional optimization of the force field parameters, with {σCH2 , ǫCH2 , σCH , ǫCH } as four degrees of freedom that are simultaneously optimized. The last degree of freedom is bond length l between the terminal united atom group CH2 (sp2 ) and the group CH(sp2 ). The carbon-carbon bond length is about 1.33˚ A and we investigate the last degree of freedom discretized with l = 1.33 ˚ A + ∆l for ∆l ∈ [0, 0.1, 0.2]˚ A. The repulsive exponent n of the Mie n-6 potential was set to n = 14. The “n-alkane parameters” CH2 (sp3 ) and CH2 (sp3 ) are used unaltered. For the optimization we again define an objective function that measures relative squared deviations in vapor pressure and liquid densities (with equal weight) of 1-butene, 1-hexene, 1-octene, cis-2-pentene and trans-2-pentene. The objective function includes α- and β-olefins in order to decouple the otherwise strongly correlated parameters of CH2 (sp2 ) and CH(sp2 ) united-atom groups. Reduced temperatures of Tr = 0.55 to 0.97 were used for all of these olefins. Experimental data for liquid density and saturated vapor pressure are actually results from multi-parameter correlations provided by the DIPPR database. 82 The optimization gave final values for the objective function of 1.92 %, 1.98 %, and 1.74 % for the three different bond length outward shift ∆l = 0.0 ˚ A, ∆l = 0.1 ˚ A, and ∆l = 0.2 ˚ A, respectively. Since the additional outward shift ∆l 6= 0 in the bond length does not improve the results significantly, as opposed to the n-alkanes, the bond length between CH2 (sp2 ) and CH(sp2 ) is chosen as l = 1.33 ˚ A (from ∆l = 0.0 ˚ A). Simulation results for the TAMie potential along with results of a force field very recently proposed by Potoff and Kamath 88 (here abbreviated as PBBK) on transferable Mie potentials for olefines are presented in Fig. 12 to 13. 28


Page 28 of 45

Page 29 of 45

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


Let’s first consider the olefins that were member of the objective function (1-butene, 1-hexene, 1-octene, cis-2-pentene and trans-2-pentene). Deviations for liquid density from the TAMie force field are 1.38 % whereas PBBK gives errors of 0.64 %. The critical point is again estimated at somewhat too high temperatures for all substances, as Fig. 12 indicates for 1-hexene, trans-2-pentene, and 1-butene. Deviations in vapor pressure for this group of substances are 3.10 % for the TAMie potential and 5.61 % for the PBBK parameterization. Fig. 13 illustrates good agreement between the TAMie force field and the experimental vapor pressure data, for 1-butene, 1-hexene, and 1-octene. For the case of 1-pentene, as a substance that was not part of our objective function, we give a comparison to other force fields, namely to TraPPE, 1,12 to NERD, 25,26 to the model of Spyriouni et al., 84 and to PBBK 46,88 in the supporting information. Fig. 12 and 13 confirm that our force field performs better for vapor pressure data but weaker for liquid density data and critical temperature, compared to the PBBK parameterization. That compromise is a result of the chosen objective function. The reason, why the TAMie parameterization is not overall better than the PBBK model indicates that for olefins the parameters {σCH2 , ǫCH2 } are largely decoupled from the {σCH , ǫCH } parameters of the sp2 -groups. That means that solving the decoupled 2 × 2 problem (determining {σCH , ǫCH } and {σCH2 , ǫCH2 }) closely resembles the full 4 dimensional problem of simultaneously determining all parameters.

29



600 trans-3-octene 1-hexene

500

T/K

trans-2-pentene 1-butene

400

300

200 0

exp.data TAMie Potoff, Bernard-Brunel

100

200

300

400

500

600

700

3

ρ / kg/m

Figure 12: Vapor-liquid coexistence curve for short and intermediate n-olefins. Symbols are simulation results for TAMie force field (blue triangles), and for the PBBK model (red diamonds), 46,88 and crosses correspond to critical experimental density. 81 Solid lines correspond to quasi-experimental data. 82,83,89,90

6

sat

rel. error (%) p (T)

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 30 of 45

4

TAMie Potoff, Bernard-Brunel

2 0 -2 -4 -6 0.5

0.6

0.7

0.8

0.9

1

T / Tc

Figure 13: Relative errors (%) in vapor pressure of short and intermediate olefins towards quasi-experimental data. 82,89,90 Symbols are simulation results for TAMie model (blue symbols connected by solid line) and the PBBK force field (red symbols connected by dashed line). 46,88 Diamonds are for 1-butene, squares for 1-pentene, triangles for 1-hexene, and spheres for 1-octene.

30


Page 31 of 45

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


6

Conclusion

This study proposes a new transferable anisotropic force field based on the Mie-potential for n-alkanes and n-olefins. In a united-atom approach, the homologous series of n-alkanes is defined by 2 + 2 parameters of the Mie-potential representing a methyl-group CH3 (sp3 ) and a methylen-group CH2 (sp3 ). Analogously, 2 + 2 additional parameters were defined for linear olefins. These 4 parameters (for each chemical family) were optimized simultaneously, applying an analytic equation of state (PC-SAFT) for improved convergence. The high efficiency of the iteration procedure allowed us to define and optimize 2 further degrees of freedom in the force field, namely the repulsive exponent n of the Mie n-6 potential as well as the outward shift ∆l of the methyl-group from the position of the carbon atom. We propose a Mie 14-6 model where the methyl-group is shifted outward from its carbon position by 0.2˚ A as a suitable force field for n-alkanes. For n-olefins, the exponents of 14-6 also turned out to be suitable, but a shift in bond length did not significantly improve the results. The resulting force field reproduces and predicts phase equilibrium properties of n-alkanes (and a binary mixture) and n-olefins in excellent agreement with experimental data.

31



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

Acknowledgement The authors thank the German Research Foundation (DFG) for financial support of the project within the Cluster of Excellence in Simulation Technology (EXC 310/2) at the University of Stuttgart. We thank Athanassios Panagiotopoulos for very helpful comments on the manuscript.

Supporting Information Available Simulation results for n-alkanes with optimized force field parameters and with varying repulsive exponent of the Mie potential n ∈ [12, 13, 14, 15, 16] and varying outward shift of the CH3 -group ∆l ∈ 0.0 ˚ A, 0.1 ˚ A, 0.2 ˚ A are provided as supporting information. A

comparison to other common force fields is given. We also include simulation results for

n-olefins for the optimized force field with n = 14 and ∆l = 0.0˚ A. This material is available free of charge via the Internet at http://pubs.acs.org/.

7

Appendix: converged results are unaffected by the analytic equation of state

We show that errors caused by the equation of state do not alter the minimum of the objective function. Although an analytic equation of state is used to accelerate the iteration of the force-field parameters, the absolute contribution of the analytic model to the objective function (as well as the gradient-contribution) disappears at the minimum of the objective function. The objective function f (p) is a function of the selected force-field parameters p. For n-alkanes, for example, we selected {σCH3 , ǫCH3 , σCH2 , ǫCH2 } as elements of p. The iteration scheme is visualized in Fig. 14. In step 1, the parameters φǫ,k and φσ,k are adjusted for every substance k using Eq. (9) and (10). The sum of squared deviations between the PC-SAFT

32


Page 32 of 45

Page 33 of 45


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


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 34 of 45

remaining part that accounts for all possible deviations. Eq. (12) corresponds to decomposing each observable Ωj (say, a vapor pressure point) according to SAFT Ωsim (p) + Ωdev j (p) = Ωj j (p)

(13)

Eq. (12) and (13) define the quantities f dev (p) and Ωdev j (p); the equations are exact by (n) ) are definition. For the parameter vector of the n-th iteration, the quantities Ωdev j (p

determined in step 1 (Fig. 14) of our method. As the only assumption of our algorithm, we dev (n) approximate the quantities Ωdev ) as constant for step 2 of our method, i.e. j (p) ≈ Ωj (p

we locally approximate SAFT (n) Ωsim (p) + Ωdev ) j (p) ≈ Ωj j (p

(14)

We will proceed to show that this assumption does not alter the optimum; towards the minimum of the objective function the assumption vanishes. Consider the optimal parameter vector p = p∗ minimizing f (p). The step 1 of Fig. 14 minimizes squared deviations Ωdev = Ωsim − ΩSAFT . The parameters of the SAFT equation j j j (m ˆ k, σ ˆk , ǫˆk ) are functions σ ˆk (p, φσ,k ) and ǫˆk (p, φǫ,k ), whereas m ˆ k is of state entering ΩSAFT j constant. Since in step 1, the values of φσ,k and φǫ,k of each substance k contained in the objective function are determined by minimizing the the squared deviations Ωdev j , we obtain zero gradients of f dev (p∗) with respect to φσ,k and φǫ,k , i.e. we have ∇φ f dev = 0 where ∇φ denotes the derivative to all φ-parameters at p∗. The product rule of differentiation requires T ˆ · ∇pˆf dev ∇φ f dev = ∇φ p

(15)

ˆ Where, for a compact notation, we introduced ∇pˆ as a derivative operator with respect to p as a vector of the relevant PC-SAFT parameters σ ˆk , ǫˆk of all components. For the n-alkanes we considered the 4 force field parameters of the Mie potential p ∈ IR4 and included 4 nalkanes in the objective function, so that p ˆ ∈ IR8 (with also 8 values of φσ,k and φǫ,k ). The 34


Page 35 of 45

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 ˆ is of rank 8. In fact, it would become a diagonal matrix when the factor 8 × 8 matrix ∇φ p

φσ,k in Eq. (10) were absorbed into φǫ,k . Consequently, from Eq. (15) we have ∇pˆf dev = 0. A second application of the product rule of differentiation, leads to T ∇pˆf dev = ∇pˆp · ∇p f dev

(16)

T For the alkanes ∇pˆp is a 8 × 4 matrix is of rank 4. The condition, that a given vector of

force field parameters p leads to different PC-SAFT parameters for the different n-alkane-

members of the objective function is fulfilled. We thus get vanishing gradients, ∇p f dev = 0 with respect to the force field parameters p at the optimum p∗. Using Eq. (12), this condition ensures that the full objective function is indeed at its minimum for ∇p f SAFT = 0, with ∇p f dev = 0 at p∗. In conclusion: once the algorithm converges, the final result is unaffected by the analytic equation of state.

References (1) Martin, M. G.; Siepmann, J. I. Transferable Potentials for Phase Equilibria. 1. UnitedAtom Description of n-Alkanes. J. Phys. Chem. B 1998, 102, 2569–2577. (2) Jorgensen, W. L.; Madura, J. D.; Swenson, C. J. Optimized Intermolecular Potential Functions for Liquid Hydrocarbons. J. Am. Soc. 1984, 106, 6638–6646. (3) Siepmann, J. I.; Karaborni, S.; Smit, B. Simulating the Critical Behaviour of Complex Fluids. Nature 1993, 365, 330–332. (4) Panagiotopoulos, A. Z. Direct Determination of Phase Coexistence Properties of Fluids by Monte Carlo Simulation in a new Ensemble. Mol. Phys. 1987, 61, 813–826. (5) Panagiotopoulos, A. Z.; Quirke, N.; Stapleton, M.; Tildesley, D. J. Phase Equilibria by Simulation in the Gibbs Ensemble. Mol. Phys. 1988, 63, 527–545.

35



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

(6) Siepmann, J. I.; Frenkel, D. Configurational Bias Monte Carlo: a new Sampling Scheme for Flexible Chains. Mol. Phys. 1992, 75, 59–70. (7) Frenkel, D.; Mooij, G. C. A. M.; Smit, B. Novel Scheme to Study Structural and Thermal Properties of Continuously Deformable Molecules. J. Phys. Cond. Matter 1992, 4, 3053–3076. (8) Smit, B.; Karaborni, S.; Siepmann, J. I. Computer Simulations of Vapor-Liquid Phase Equilibria on n-Alkanes. J. Chem. Phys. 1995, 102(5), 2126–2140. (9) Siepmann, J. I.; Martin, M. G.; Mundy, C. J.; Klein, M. L. Intermolecular Potentials for Branched Alkanes and the Vapour-Liquid Phase Equilibria of n-heptane, 2methylhexane, and 3-ethylpentane. Mol. Phys. 1997, 90(5), 687–693. (10) Martin, M. G.; Siepmann, J. I. Predicting Multicomponent Phase Equilibria and Free Energies of Transfer for Alkanes by Molecular Simulation. J. Am. Chem. Soc. 1997, 119, 8921–8924. (11) Martin, M. G.; Siepmann, J. I. Novel Configurational-Bias Monte Carlo Method for Branched Molecules. Tranferable Potentials for Phase Equilibria. 2. United-Atom Description of Branched Alkanes. J. Phys. Chem. B 1999, 103, 4508–4517. (12) Wick, C. D.; Martin, M. G.; Siepmann, J. I. Transferable Potentials for Phase Equilibria. 4. United Atom Description of Linear and Branched Alkenes and Alkylbenzenes. J. Phys. Chem. B 2000, 104, 8008–8016. (13) Chen, B.; Potoff, J. J.; Siepmann, J. I. Monte Carlo 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.

36


Page 36 of 45

Page 37 of 45

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


(14) Stubbs, J. M.; Potoff, J. J.; Siepmann, J. I. Transferable Potential for Phase Equilibria. 6. United-Atom Description for Ethers, Glycols, Ketones, and Aldehydes. J. Phys. Chem. B 2004, 108, 17596–17605. (15) Wick, C. D.; Stubbs, J. M.; Rai, N.; Siepmann, J. I. Transferable Potential for Phase Equilibria. 7. Primary, Secondary, and Tertiary Amines, Nitroalkanes and Nitrobenzenes, Nitriles, Amides, Pyridine, and Pyrimidine. J. Phys. Chem. B 2005, 109, 18974– 18982. (16) Kamath, N. L. G.; Potoff, J. J.; Rai, N.; Siepmann, J. I. Transferable Potential for Phase Equilibria. 8. United-Atom Description for Thiols, Sulfides, Disulfides, and Thiophene. J. Phys. Chem. B 2005, 109, 24100–24107. (17) Ketko, M. H.; Rafferty, J.; Siepmann, J. I.; j. J. Potoff, Development of the TraPPE-UA Force Field for Ethyleneoxide. Fluid Phase Equilibria 2008, 274, 44–49. (18) Maerzke, K. A.; Schultz, N. E.; Ross, R. B.; Siepmann, J. I. TraPPE-UA Fore Field for Acrylates and Monte Carlo Simulations for Their Mixtures with Alkanes and Alcohols. J. Phys. Chem. B 2009, 113, 6415–6425. (19) Bhatnaga, N.; Kamath, G.; Potoff, J. J. Biomolecular Simulations with the Transferable Potentials for Phase Equilibria: Extension to Phospholipids. J. Phys. Chem. B 2013, 117, 9910–9921. (20) Jorgensen, W. L.; Swenson, C. J. Optimized Intermolecular Potential Functions for Amides and Peptides. Structure and Properties of Liquid Amides. J. Am. Soc. 1985, 107, 569–578. (21) Jorgensen, W. L.; Swenson, C. J. Optimized Intermolecular Potential Functions for Amides and Peptides. Hydration of Amides. J. Am. Soc. 1985, 107, 1489–1496.

37



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

(22) Jorgensen, W. L. Optimized Intermolecular Potential Functions for Liquid Alcohols. J. Phys. Chem. 1986, 90, 1276–1284. (23) Briggs, J. M.; Matsui, T.; Jorgensen, W. L. Monte Carlo Simulations of Liquid Alkyl Ethers with the OPLS Potential Functions. J. Comp. Chem 1990, 11(8), 958–971. (24) Jorgensen, W. L.; Tirado-Rives, J. The OPLS Force Field for Proteins. Energy Minimization for Crystals of Cyclic Peptides and Crambin. J. Am. Chem. Soc. 1988, 110(6), 1657–1666. (25) Nath, S.; Escobeda, F. A.; de Pablo, J. J. On the Simulation of Vapour-Liquid Equilibria for Alkanes. J. Chem. Phys. 1998, 108, 9905–9911. (26) Nath, S.; Banaszak, B. J.; de Pablo, J. J. A new United Atom Force Field for alphaOlefins. J. Chem. Phys. 2001, 114, 3612–3616. (27) Nath, S. K.; Khare, R. New Forcefield Parameters for Branched Hydrocarbons. J. Chem. Phys 2001, 115(23), 10837–10844. (28) Khare, R.; Sum, A. K.; Nath, S. K.; de Pablo, J. J. Simulation of Vapor-Liquid Phase Equilibria of Primary Alcohols and Alcohol-Alkane Mixtures. J. Phys. Chem. B 2004, 108, 10071–10076. (29) Toxvaerd, S. Equation of State of Alkanes II. J. Chem. Phys. 1997, 107(13), 5197 –5204. (30) Toxvaerd, S. Molecular Dynamics Calculation of the Equation of State of Alkanes. J. Chem. Phys. 1990, 93(6), 4290–4295. (31) Ungerer, P.; Beauvais, C.; Delhommelle, J.; Boutin, A.; Rousseau, B. Optimization of the Anisotropic United Atoms Intermolecular Potential for n-Alkanes. J. Chem. Phys. 2000, 112, 5499–5510.

38


Page 38 of 45

Page 39 of 45

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


(32) Bourasseau, E.; Ungerer, P.; Boutin, A.; Fuchs, A. H. Monte Carlo Simulation of Branched Alkanes and Long Chain n-Alkanes with Anisotropic United Atoms Intermolecular Potential. Mol. Sim. 2002, 28(4), 317–336. (33) Bourasseau, E.; Ungerer, P.; Boutin, A.; Mackie, A. D. Prediction of Equilibrium Properties of Cyclic Alkanes by Monte Carlo Simulation - New Anistotopic United Atoms Intermolecular Potential - New Transfer Bias Method. J. Phys. Chem. B 2002, 106, 5483–5491. (34) Kranias, S.; Pattou, D.; Levy, B.; Boutin, A. An Optimized Potential for Phase Equilibria Calculation for Ketone and Aldehyde Molecular Fluids. Phys. Chem. Phys. 2003, 5, 4175–4179. (35) Bourasseau, E.; Haboudou, M.; Boutin, A.; Fuchs, A. H.; Ungerer, P. New Optimization Method for Intermolecular Potentials: Optimization of a new Anisotropic United Atom Potential for Olefins: Prediction of Equilibrium Properties. J. Chem. Phys. 2003, 18(7), 3020–3034. (36) Contreras-Camacho, R. O.; Ungerer, P.; Ahunbay, M. G.; Lachet, V.; Perez-Pellitero, J.; Mackie, A. D. Optimized Intermolecular Potential for Aromatic Hydrocarbons Based on Anisotropic United Atoms. 2. Alkylbenzenes and Styrene. J. Phys. Chem. B 2004, 108, 14115–14123. (37) Contreras-Camacho, R. O.; Ungerer, P.; Boutin, A.; Mackie, A. D. Optimized Intermolecular Potential for Aromatic Hydrocarbons Based on Anisotropic United Atoms. 1. Bezene. J. Phys. CHem. B 2004, 108, 14109–14114. (38) Nieto-Draghi, C.; Ungerer, P. Optimization of the Anisotropic United Atoms Intermolecular Potential for n-Alkanes: Improvement of Transport Properties. J. Chem. Phys. 2006, 125, 044517.

39



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

(39) Nieto-Draghi, C.; Bonnaud, P.; Ungerer, P. Anisotropic United Atom Model Including the Electrostatic Interactions of Methylbenzenes. II. Transport Properties. J. Phys. Chem. C 2007, 111, 15942–15951. (40) Nieto-Draghi, C.; Bonnaud, P.; Ungerer, P. Anisotropic United Atom Model Including the Electrostatic Interactions of Methylbenzenes. I. Thermodynamic and Structural Properties. J. Phys. Chem. C 2007, 111, 15686–15699. (41) Biscay, F.; Ghoufi, A.; Goujon, F.; LAchet, V.; Malfreyt, P. Surface Tension of Linear and Branched Alkanes from Monte Carlo Simulations Using the Anisotropic Uited Atom Model. J. Phys. Chem. B 2008, 112, 13885–13897. (42) Biscay, F.; Ghoufi, A.; LAchet, V.; Malfreyt, P. Calculation of the Surface Tension of Cyclic and Aromatic Hydrocarbons from Monte Carlo Simulations using an Anisotropic United Atom Model (AUA). Phys. Chem. Chem. Phys. 2009, 11, 6132–6147. (43) Ferrando, N.; Lachet, V.; Teuler, J.; Boutin, A. Transferable Force Field for Alcohols and Polyalcohols. J. Phys. Chem. B 2009, 113, 5985–5995. (44) Ferrando, N.; Lachet, V.; Perez-Pelliter, J.; d. Mackie, A.; Malfreyt, P.; Boutin, A. A Transferable Force Field to Predict Phase Equilibria and Surface Tension of Ethers and Glycol Ethers. J. Phys. Chem. B 2011, 115, 10654–10644. (45) Errington, J. R.; Panagiotopoulos, A. Z. A new Intermolecular Potential Model for the n-Alkane Homologous Series. The Journal of Physical Chemistry B 1999, 103, 6314–6322. (46) Potoff, J. J.; Bernard-Brunel, D. A. Mie Potentials for Phase Equiilbria Calculations: Application to Alkanes and Perfluoralkanes. J. Phys. Chem. B 2009, 113, 14725–14731. (47) Faller, R.; Schmitz, H.; Biermann, O.; Mueller-Plathe, F. Automatic Parametrization of Force Fields for Liquids by Simplex Algorithm. J. Comp. Chem. 1999, 20, 1009. 40


Page 40 of 45

Page 41 of 45

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


(48) H¨ ulsmann, M.; K¨oddermann, T.; Vrabec, J.; Reith, D. GROW: A Gradient-based Optimization Workflow for Automated Development of Molecular Models. Comp. Phys. Comm. 2010, 181, 499–513. (49) v. Westen, T.; Vlugt, T. J. H.; Gross, J. Determining Force Field Parameters Using a Physically Based Equation of State. J. Phys. Chem. B 2011, 115, 7872–7880. (50) Gross, J.; Sadowski, G. Perturbed-Chain SAFT: An Equation of State Based on a Perturbation Theory for Chain Molecules. Ind. Eng. Chem. Res. 2001, 40, 1244–1260. (51) Gray, C. G.; Gubbins, K. E. Theory of molecular fluids; 1984. ¨ (52) Lorentz, H. A. Uber die Anwendung des Satzes vom Virial in der kinetischen Gastheorie. Ann. Phys. 1881, 12 . (53) Berthelot, D. C. R. Sur ur le m´elange des gaz. Hebd. Seanc. Acad. Sci. (Paris) 1898, 126, 1703–1855. (54) d. Ploeg, P. V.; Berendsen, H. J. C. Molecular dynamics simulation of a bilayer membrane. J. Chem. Phys. 1982, 76(6), 3271–3276. (55) Cornell, W. D.; et al., A Second Generation Forc Field for the Simulation of Proteins, Nucleic Acids, and Organic Molecules. J. Am. Chem. Soc. 1995, 117, 5179–5197. (56) Allen, M. P.; Tildesley, D. J. Computer simulation of liquids; 1987. (57) Sanchez, J. M. C.; Danner, T.; Gross, J. Grand Canonical Monte Carlo Simulations of Vapor-Liquid Equilibria Using a Bias Potential from an Analytic Equation of State. J. Chem. Phys. 2013, 138, 234106. (58) Virnau, P.; M¨ uller, M. Calculation of Free Energy through Successive Umbrella Sampling. J. Chem. Phys. 2004, 120, 10925–10930.

41



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

(59) Errington, J. R. Direct Calculation of Liquid-Vapor Phase Equilibria from Transition Matrix Monte Carlo Simulation. J. Chem. Phys. 2003, 118, 9915–9925. (60) Paluch, A. S.; Shen, V. K.; Errington, J. R. Comparing the Use of Gibbs Ensemble and Grand-Canonical Transition-Matrix Monte Carlo Methods to Determine Phase Equilibria. Ind. Eng. Chem. Res. 2008, 47, 45233–4541. (61) Ferrenberg, A. M.; Swendsen, R. H. New Monte Carlo Technique for Studying Phase Transitions. Phy. Rev. Lett. 1988, 61, 2635. (62) Ferrenberg, A. M.; Swendsen, R. H. Optimized Monte Carlo Data Analysis. Phy. Rev. Lett. 1989, 63, 1195–1198. (63) Hemmen, A.; ; Panagiotopoulos, A. Z.; Gross, J. Grand Canonical Monte Carlo Simulations Guided by an Analytic Equation of State - Transferable Anisotropic Mie Potentials for Ethers. J. Phys. Chem. B 2015, 119, 7087–7099. (64) Paluch, A. S.; Shen, V. K.; Errington, J. R. Determination of Fluid-Phase Behavior Using Transition-Matrix Monte Carlo: Binary Lennard-Jones Mixtures. Ind. Eng. Chem. Res. 2008, 47, 45233–4541. (65) Smit, B. Phase diagrams of LennardJones fluids. J. Chem. Phys. 1992, 96, 8639–8640. (66) Schacht, C. S.; Vlugt, T. H. H.; Gross, J. Using an Analytic Equation of State to Obtain Quantitative Solubilities of CO2 by Molecular Simulations. J. Phys. Chem. Lett. 2011, 2, 393–396. (67) Gross, J.; Sadowski, G. Modeling Polymer Systems Using the Perturbed-Chain Statistical Associating Fluid Theory Equation of State. Ind. Eng. Chem. Res. 2002, 41, 1084–1093. (68) Gross, J.; Sadowski, G. Application of the Perturbed-Chain SAFT Equation of State to Associating Systems. Ind. Eng. Chem. Res. 2002, 41, 5510–5515. 42


Page 42 of 45

Page 43 of 45

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


(69) Stoll, J.; Vrabec, J.; Hasse, H. Comprehensive Study of the Vapour–Liquid Equilibria of the Pure Two-Centre Lennard–Jones plus Pointdipole Fluid. Fluid phase equilibria 2003, 209, 29–53. (70) Stoll, J.; Vrabec, J.; Hasse, H.; Fischer, J. Comprehensive Study of the Vapour–Liquid Equilibria of the Pure Two-Centre Lennard–Jones plus Pointquadrupole fluid. Fluid Phase Equilibria 2001, 179, 339–362. (71) Vrabec, J.; Stoll, J.; Hasse, H. A set of Molecular Models for Symmetric Quadrupolar Fluids. The Journal of Physical Chemistry B 2001, 105, 12126–12133. (72) Stoll, J.; Vrabec, J.; Hasse, H. A set of Molecular Models for Carbon Monoxide and Halogenated Hydrocarbons. The Journal of chemical physics 2003, 119, 11396–11407. (73) Lymperiadis, A.; Adjiman, C. S.; Galindo, A.; Jackson, G. A Group Contribution Method for Associating Chain Molecules Based on the Statistical Associating Fluid Theory (SAFT-γ). The Journal of chemical physics 2007, 127, 234903. (74) Lymperiadis, A.; Adjiman, C. S.; Jackson, G.; Galindo, A. A Generalisation of the SAFT-group Contribution Method for Groups Comprising Multiple Spherical Segments. Fluid Phase Equilibria 2008, 274, 85–104. (75) Papaioannou, V.; Lafitte, T.; Avenda˜ no, C.; Adjiman, C. S.; Jackson, G.; M¨ uller, E. A.; Galindo, A. Group Contribution Methodology based on the Statistical Associating Fluid Theory for Heteronuclear Molecules formed from Mie Segments. The Journal of chemical physics 2014, 140, 054107. (76) Avendano, C.; Lafitte, T.; Galindo, A.; Adjiman, C. S.; Jackson, G.; Muller, E. A. SAFT-γ Force Field for the Simulation of Molecular Fluids. 1. A Single-Site Coarse Grained Model of Carbon Dioxide. The Journal of Physical Chemistry B 2011, 115, 11154–11169.

43



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

(77) Avendano, C.; Lafitte, T.; Adjiman, C. S.; Galindo, A.; Muller, E. A.; Jackson, G. SAFT-γ Force Field for the Simulation of Molecular Fluids: 2. Coarse-Grained Models of Greenhouse Gases, Refrigerants, and long Alkanes. The Journal of Physical Chemistry B 2013, 117, 2717–2733. (78) Lafitte, T.; Avenda˜ no, C.; Papaioannou, V.; Galindo, A.; Adjiman, C. S.; Jackson, G.; M¨ uller, E. A. SAFT-γ Force Field for the Simulation of Molecular Fluids: 3. CoarseGrained Models of Benzene and Hetero-Group Models of n-Decylbenzene. Molecular Physics 2012, 110, 1189–1203. (79) Sauer, E.; Stavrou, M.; Gross, J. Comparison between a Homo- and a Heterosegmented Group Contribution Approach based on the Perturbed-Chain Polar Statistical Associating Fluid Theory Equation of State. submitted (80) Sauer, E.; Stavrou, M.; Gross, J. Comparison between a Homo-and a Heterosegmented Group Contribution Approach Based on the Perturbed-Chain Polar Statistical Associating Fluid Theory Equation of State. Ind. Eng. Chem. Res. 2014, 53, 14854–14864. (81) NIST Chemistry WebBook. 2013; http://webbook.nist.gov/chemistry. (82) Rowley, R. L.; Wilding, W. V.; Oscarson, J. L.; Yang, Y.; Zundel, N. A.; Daubert, T. E.; Danner, R. P. DIPPR Data Compilation of Pure Chemical Properties. Design Institute for Physical Properties. 2009. (83) Peng, D. Y.; Robinson, D. B. A New Two-Constant Equation of State. Ind. Eng. Chem. 1976, 15, 59–64. (84) Spyriouni, T.; Economou, I. G.; Theodorou, D. N. Molecular Simulation of α-Olefins Using a New United-Atom Potential Model: Vapor-Liquid Equilibria of Pure Compounds and Mixtures. Journal of the American Chemical Society 1999, 121, 3407–3413.

44


Page 44 of 45

Page 45 of 45

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


(85) Vrabec, J.; Stoll, J.; Hasse, H. A set of Molecular Models for Symmetric Quadrupolar Fluids. J. Phys. Chem. B 2001, 105, 12126–12133. (86) Perez-Pellitero, J.; Ungerer, P.; Orkoulas, G.; Mackie, A. D. Critical Point Estimation of the Lennard-Jones Pure Fluid and Binary Mixtures. J. Chem. Phys. 2006, 125, 054515. (87) Sage, B. H.; Hicks, B. L.; Lacey, W. N. Phase Equilibria in Hydrocarbons. Ind. Eng. Chem. Ind. Ed. 1940, 32(8), 1085–1092. (88) Potoff, J. J.; Kamath, G. Mie Potentials for Phase Equilibria: Application to Alkenes. J. Chem. Eng. Data 2014, 59, 3144–3150. (89) Wagner, W. New Vapour Pressure Measurements for Argon and Nitrogen and a new Method for Establishing Rational Vapour Pressure Equations. Cryogenics 1973, 13(8), 470. (90) Smith, B. D.; Srivastava, R. Thermodynamic data for pure compounds: Part A Hydrocarbons and Ketones; 1986.

45


Transferable Anisotropic United-Atom Force Field Based on the Mie

Recommend Documents