Transferable Anisotropic United-Atom Force Field Based on the Mie

Article pubs.acs.org/JPCB

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, Germany

Downloaded via UNIV OF CALIFORNIA SANTA BARBARA on July 4, 2018 at 00:44:16 (UTC). See https://pubs.acs.org/sharingguidelines for options on how to legitimately share published articles.

S Supporting Information *

ABSTRACT: A new transferable force field parametrization 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 toward 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 parametrization 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, 1-octene, 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 PC-SAFT 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, Tr = 0.55−0.97. For substances that were not members of the objective function, 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 parametrization of effective pairwise force fields is particularly difficult to estimate from quantum mechanical calculations. A critical test for the parametrization of the van der Waals contribution to the force field is the 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, © 2015 American Chemical Society

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 technique4,5 and configurational bias Monte Carlo6,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 co-workers.1,11−19 The TraPPE force field is based on intramolecular and electrostatic models of OPLS-UA (optimized potentials for liquid simulations) proposed by Jorgensen2 and uses a Lennard-Jones approach for the van der Waals interactions. The Lennard-Jones parameters were fitted to liquid density and to critical temperatures of pure Received: February 9, 2015 Revised: July 30, 2015 Published: August 4, 2015 11695

DOI: 10.1021/acs.jpcb.5b01354 J. Phys. Chem. B 2015, 119, 11695−11707

Article

The Journal of Physical Chemistry B 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 OPLSUA 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 field25,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 models29 that show promising results). Toxvaerd30 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 parametrization 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 alkanes31 and olefins35 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 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 Lennard-Jones 12-6 potential.45 They obtained very good results using a Buckingham exponential-6 potential, where the repulsive part of the potential is different from the Lennard-Jones potential. With a rationale similar to Toxvaerd’s,30 they shifted the CH3 group outward 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 ρ(Tcoex).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 optimization of the LennardJones parameters of a united-atom force field to those for n-alkanes (i.e., CH3 and CH2), the problem is 4-dimensional. Most studies have parametrized 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 variety of atomistic force fields, Ungerer and co-workers31,35 used a direct gradient approach for only small molecules, and Hülsmann 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 state50 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 PC-SAFT 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−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 the 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 section briefly explains the optimization algorithm, and the last section presents and discusses the simulation results of the TAMie force field in comparison with experimental data and with other force fields.

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 centers of interaction (center of force) for both CH3 and CH2 groups 11696


Article

The Journal of Physical Chemistry B 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 to simply place 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 ⎡⎛ σ ⎞ n ⎛ σ ⎞ m ⎤ ij ij u(rij) = cn ϵij⎢⎜⎜ ⎟⎟ − ⎜⎜ ⎟⎟ ⎥ ⎢⎣⎝ rij ⎠ ⎝ rij ⎠ ⎥⎦

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

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. Table 1. Parameters of Intramolecular Potential for TAMie Force Field, Taken from Literature Bond Length1,12

(1)

⎛ n ⎞⎛ n ⎞ m / n − m ⎜ ⎟⎜ ⎟ ⎝ n − m ⎠⎝ m ⎠

CH2CH2 CH3CH2 CH(sp2)CHz(sp2) z ∈ [1,2] CH3CH(sp2) CH3CH3 ethane Bond Angle8,12,54,55

(2)

where rij denotes the distance between two interaction sites i and j, and ϵij and σij are the well depth and diameter parameter, respectively. 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 three bonds. Parameters for interactions between unlike types of interaction sites are determined using Lorentz−Berthelot combining rules,52,53 as σij = (σii + σjj)/2

ϵij =

ϵiiϵjj

kΘ/kB/K

114.0 CHxCH2CHy x,y ∈ [2,3] CHzCHCHy z ∈ [1,2], y ∈ [2,3] 119.7 Torsional Potential1,2,12 torsion type

62 500 70 420

c0/kB/K

c1/kB/K

c2/kB/K

c3/kB/K

0

355.03

−68.19

791.32

688.5

86.36

−109.77

−282.24

CHxCH2CH2CHy x, y ∈ [2,3] CHzCHCH2CHy z ∈ [1,2], y ∈ [2,3]

CHxCHCHCHy (cis) x,y ∈ [2,3] CHxCHCHCHy (trans) x,y ∈ [2,3]

d0/kB/K

ϕ/rad

12 400 13 400

π 0

Table 2. Parameters for Van der Waals Interactionsa 12-6 LJ

(4)

TAMie

pseudoatom

ϵ/kB/K

σ/Å

n

ϵ/kB/K

σ/Å

n

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

103.215 45.44

3.746 4.051

12 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

(3)

Energy parameter ϵ, size parameter σ, and repulsive exponent n the proposed of TAMie model and for an optimized 12-6 Lennard-Jones potential. a

The remaining parameters of the TAMie potential are summarized in Table 2. Furthermore, optimized 12-6 LennardJones parameters are proposed, and were optimized to reproduce saturated vapor pressure and liquid densities of vapor−liquid equilibria using intramolecular potentials and constants as listed above.

(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, we use the OPLS-UA model,2 according to

3. MOLECULAR SIMULATION TECHNIQUE This section first reviews the Monte Carlo simulation technique for phase equilibria of pure components and binary mixtures and, second, 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,

utor(ϕ) = c0 + c1[1 + cos(ϕ)] + c 2[1 − cos(2ϕ)] + c3[1 + cos(3ϕ)]

1.54 1.54 + 0.2 1.33 1.54 + 0.2 1.54 + 2·0.2 Θ0/deg

bond type

Prefactor cn in eq 2 can also be defined for groups of cross-wise different powers nii and njj. Potoff and Bernard-Brunel46 propose an arithmetic combination rule nij = 0.5(nii + njj) to be considered in eq 2. For the intramolecular force field, we adopt a parametrization 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 potential54 ubend(θ ) = kθ /2(θ − θ0)2

r0/Å

bond type

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

(7)

(6) 11697


Article

The Journal of Physical Chemistry B μ, 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 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 2-fold: first, the calculation time for each window is low, and all windows can run in parallel. Second, 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 cannot cover wide ranges of temperatures away 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 transitionmatrix Monte Carlo64 (TMMC) scheme in a multicanonical ensemble (N1, N2, V, T) for a binary mixture of components 1 and 2. We divide the N1-N2-surface into ΔN1 × ΔN2 windows of 11 × 11 molecules. Exactly 1/121st part of the MCproduction steps are spent for every N1, N2-combination, sampling the transition probabilities for insertion- and deletiontrials. 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 3.2. Simulation Details. The grand canonical Monte Carlo simulations were conducted with box lengths L of cubic simulation boxes chosen as 39 Å for (C3−C6), 43 Å for (C8), up to 61 Å for (C24). The Mie potential is calculated for all interaction pairs up to a radial cutoff of 14 Å with analytical long-range corrections56,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 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 Westen49 (and by Schacht et al. for mixtures66) is modified and is applied to optimize the intermolecular force field parameters. The PC-SAFT equation of state50,67,68 is used to correlate and subsequently predict results of molecular simulations, allowing for 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 (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 2-center Lennard-Jones fluids with an embedded point dipole69 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. Müller, Jackson, Adjiman, and Galindo together with colleagues used the SAFT-γ model73−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 constrained. 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 classical atomistic force fields (on a united-atom level), whereas SAFT models do not resolve the same molecular detail. It is possible, nonetheless, to support 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 results 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 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 unitedatom 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 PC-SAFT parameters {m̂ , σ̂, ϵ̂} is given in three equations. These equations define the three pure component parameters of the PC-SAFT model, i.e., segment 11698


Article

The Journal of Physical Chemistry B

eqs 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 16-6 potential leads to values of parameter ϕϵ in eq 10 quite different from unity. Nonetheless, the convergence of the algorithm did not deteriorate. Equations 9 and 10 provide some rationale for this robustness. Equation 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 canceled 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 is in the range Tr = 0.55−0.97. Force field parameters of different united-atom 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.

number m̂ , segment diameter parameter σ̂, and segment energy parameter ϵ̂, according to united ‐ atomtypes

m̂ =

∑

Nm i 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 ‐ atomtypes

̂ σ )3 = m̂ (σϕ

3 Nm i i σii

∑ i

(9)

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 ‐ atypes u ‐ atypes

̂ σ )3 (ϵ̂ϕϵ) = m̂ 2(σϕ

∑

∑

i

j

3 NN i jσij ϵij

(10)

The two parameters ϕσ and ϕϵ in eqs 9 and 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 toward the chosen group-contribution parameters m̂ i of the PC-SAFT model. With the three relations, eqs 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 observables, Ωsim j (p), where index j counts the observed data points, j = 1,..., Nexp. Now, the two parameters ϕσ,k and ϕϵ,k of the PC-SAFT model according to eqs 9 and 10 (for every substance k considered in the objective function) are adjusted, so that the PC-SAFT model very closely represents the properties Ωsim j (p) . As a modification to our previous work49 we now capture all deviations between Ωsim j (p) and PC-SAFT model in a vector Ωdev . Next, the PC-SAFT model j minimizes the objective function 1 f (p) = exp N

⎛ Ωsim(p) − Ωexp ⎞2 j j ⎟ ∑ ⎜⎜ exp ⎟ Ω j ⎠ j=1 ⎝

5. RESULTS AND DISCUSSION 5.1. n-Alkanes. The homologous series of n-alkanes is important for transferable force fields, because the pair potentials 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 Å) shown in Figure 1. The 4-dimensional iteration converges quickly, and it can be terminated after 4−5 steps, usually. While the objective function for the optimization is defined as relative mean squared deviations, according to eq 11, in Figure 1 and in the entire section, we report average absolute deviations (AADs), as

N exp

(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 j may appear as a subtlety. But in the Appendix we show that the final 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 PC-SAFT

100 AAD %= exp N

N exp

∑ j=1

exp Ωsim j (p) − Ω j

Ωexp j

The parameters {σCH3, ϵCH3, σCH2, ϵCH2} are strongly correlated, and the swift convergence (Figure 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 11699


Article


Table 3. Summary of Average Absolute Deviations (AADs) in Vapor Pressure, Liquid Densities, and Apparent (SystemSize 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 %

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 Å and the data are for propane, n-butane, n-hexane, and n-octane.

(simultaneously adjusting parameters to a collection of substances), it is known that local minima exist.79 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 nonsmooth. The appearance of local minima for a certain range of parameters (and for merely Mie based force fields) is thus a result of how the objective function is defined. From all our observations, we have the conjecture that our objective function is 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. 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. Figure 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, with vapor pressure corresponding to Figure 2.

n

Δl/Å

psat

ρL

Tc

12 12 12 13 13 13 14 14 14 15 15 15 16 16 16

0 0.1 0.2 0 0.1 0.2 0 0.1 0.2 0 0.1 0.2 0 0.1 0.2

5.68 4.97 3.79 4.46 3.50 2.94 1.87 1.21 1.07 1.12 1.34 1.62 1.79 2.11 3.33

2.38 2.02 1.70 1.79 1.78 1.54 1.47 1.23 0.90 0.99 0.97 0.73 0.85 0.75 0.75

3.40 3.06 2.48 2.38 2.44 2.15 1.90 1.66 1.12 1.23 1.19 0.80 0.72 0.52 0.70

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 from 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 parametrization strongly sacrifices the critical point. It has previously been reported that it is not possible with the Lennard-Jones potential and Δl = 0 Å to simultaneously represent all of the above properties with good accuracy. The value of the objective function of about 4.0 AAD% reflects this problem. Figure 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. Figure 3 shows

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 }. 11700


Article


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 LennardJones (njj = 12) potential. Considering all of the above, we propose a force field with n = 14 and Δl = 0.2 Å for the n-alkane family. Figures 5 −10 present coexisting densities and saturated vapor pressures of selected n-alkanes as obtained from the 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 Å. Parameters {σCH3, ϵCH3, σCH2, ϵCH2} are in the optimized state.

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 parametrization. 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 Figure 2 shows). Figure 4 illustrates how the deviations in

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

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] Å and n = 14. Parameters {σCH3, ϵCH3, σCH2, ϵCH2} are in an optimized state.

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

the vapor pressure of n-hexane change, with varying Δl. We note that although only the results for n-hexane are shown in Figures 3 and 4, the parameter regression was done for 4 alkanes simultaneously in all cases. A summary of errors in vapor pressure, liquid density, and 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 Supporting Information. According to Figure 2 and Table 3, we find equally good regression results for n = 15 and Δl = 0.2 Å and for n = 14 and Δl = 0.2 Å, when limiting consideration to integer values of n. Figure 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 Å was here considered in order to not overdraw the physical

TAMie force field in comparison to experimental data. Experimental data of the smaller n-alkanes are actually correlation results from accurate multiparameter equations of state as provided by the NIST Chemistry webbook.80 Quasiexperimental data for hexadecane, eicosane, and tetracontane are results from multiparameter correlations provided by the DIPPR database,81 and quasiexperimental values for the vapor density are estimated from the Peng−Robinson eqauation of state.82 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 11701


Article


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

Figure 7. Vapor−liquid coexistence curve for longer n-alkanes. Symbols are simulation results for TAMie force field (△), solid lines correspond to experimental data from refs 80 and 81.

calculated as absolute average deviations with respect to experimental data. Let us first consider n-alkanes that were members 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 Figures 5, 6 and 8, 9). Deviations in vapor pressure are 1.08% for the TAMie force field, and 3.79% for the PBB parametrization. We mention that for these numbers we have performed own simulations for the PBB force field because we achieved better results for the PBB model (at low temperatures) than the results published by the developers46 of the force field. We note further that our parametrization for n = 16 and Δl = 0 Å is slightly different from the PBB parametrization due to a different objective function. For n-pentane, n-decane, and 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. Figures 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.,83 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 Figures 7 and 10 indicate. Figures 5 and 8 also show 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, 0.1, and 0.2 Å. The best result was also found for n = 14 and 0.2 Å 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 Figures 5 and 8 since our optimization procedure essentially confirmed a parametrization proposed by Potoff and BernardBrunel.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

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

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

point is estimated at temperatures too high toward 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−0.97, and deviations are 11702


Article

The Journal of Physical Chemistry B (14-6) force field for butane and adopt a Lennard-Jones parametrization proposed by Vrabec et al.84 for methane (σCH4 = 3.7281 Å, ϵCH4 = [148.55]K, n = 12). We consider the combining rules, eqs 3 and 4, without any correction for cross interactions between 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.85 Figure 11 gives the simulation results for the methane−butane

otherwise strongly correlated parameters of CH2(sp2) and CH(sp2) united-atom groups. Reduced temperatures of Tr = 0.55−0.97 were used for all of these olefins. Experimental data for liquid density and saturated vapor pressure are actually results from multiparameter correlations provided by the DIPPR database.81 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 Å, Δl = 0.1 Å, and Δl = 0.2 Å, respectively. Since the additional outward shift Δl ≠ 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 Å (from Δl = 0.0 Å). Simulation results for the TAMie potential along with results of a force field very recently proposed by Potoff and Kamath87 (here abbreviated as PBBK) on transferable Mie potentials for olefins are presented in Figures 12 and 13. Let us first consider the olefins that were member of the objective function (1-butene, 1-hexene, 1-octene, cis-2-pentene,

Figure 11. Vapor−liquid equilibrium of n-butane−methane mixtures at two temperatures, T = 344.26 and 394.26 K. Comparison of molecular simulation results using the TAMie model (△) to experimental data86 (×).

mixture at two temperatures, T = 344.26 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 nonideal mixture. 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 Å and we investigate the last degree of freedom discretized with l = 1.33 Å + Δl for Δl ∈ [0, 0.1, 0.2] Å. 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

Figure 12. Vapor−liquid coexistence curve for short and intermediate n-olefins. Symbols are simulation results for TAMie force field (blue △), and for the PBBK model symbols (red ◇)46,87 and crosses (×) correspond to critical experimental density.80 Solid lines correspond to quasiexperimental data.81,82,88,89

Figure 13. Relative errors (%) in vapor pressure of short and intermediate olefins toward quasiexperimental data.81,88,89 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,87 Diamonds are for 1-butene, squares for 1-pentene, triangles for 1-hexene, and spheres for 1-octene. 11703


Article


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 Å 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.

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 Figure 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 parametrization. Figure 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.,83 and to PBBK46,87 in the Supporting Information. Figures 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 parametrization. That compromise is a result of the chosen objective function. The reason that the TAMie parametrization 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.

■

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 Figure 14. In step 1, the parameters ϕϵ,k and ϕσ,k are adjusted for every substance k using eqs 9 and 10. The sum of squared deviations between the PC-SAFT model and the simulation results is thereby minimized. The deviations (n) Ωdev j (p ) are stored in memory for step 2. Index n thereby counts the iteration steps. The force field parameters are subsequently optimized by minimizing the objective function as predicted by the PC-SAFT equation of state (step 2). For this (n) step, we treat the deviations Ωdev j (p ) as constant. The result is (n+1) . For this set of force field a new force field vector p parameters, p(n+1), new molecular simulations are conducted. This step, labeled “next iteration” in Figure 14, is again followed by step 1. The procedure is repeated until convergence is achieved.

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 methylene 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,

Figure 14. Illustration of the iteration scheme. Top row conceptually shows some experimental observables Ωexp (●), the result of molecular j SAFT (bold red line) . The lower row conceptually illustrates the objective simulations, Ωsim j (blue ■), and results of the PC-SAFT equation of state Ωj function as a function of the force field parameters p (black line) . Step1: Starting from a current estimate of the force field parameters p(n) (where molecular simulations were performed), the parameters ϕϵ,k and ϕσ,k are adjusted for every substance k using eqs 9 and 10. Step 2: The analytic equation of state is used to minimize the objective function f(p). The result serves as a new guess p(n+1), for which molecular simulations are performed (“next iteration”). 11704


Article


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.

We now analyze the approximation of our algorithm. The sim exp 2 1 N exp ⎛ Ω j (p) − Ω j ⎞ ⎟ , as defined objective function f (p) = N exp ∑ j = 1 ⎜ exp Ωj ⎝ ⎠ in eq 11, can be decomposed into f (p) = f SAFT (p) + f dev (p)

■

(12)

S Supporting Information *

This equation expresses the objective function into a part, as predicted by PC-SAFT and a remaining part that accounts for all possible deviations. Equation 12 corresponds to decomposing each observable Ωj (say, a vapor pressure point) according to SAFT Ωsim (p) + Ωdev j (p) = Ω j j (p)

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.jpcb.5b01354. 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 Å,0.1 Å,0.2 Å]; comparison to other common force fields; simulation results for n-olefins for the optimized force field with n = 14 and Δl = 0.0 Å (PDF)

(13)

Equations 12 and 13 define the quantities f dev(p) and Ωdev j (p); the equations are exact by definition. For the parameter vector (n) of the n-th iteration, the quantities Ωdev j (p ) are determined in step 1 (Figure 14) of our method. As the only assumption of our algorithm, we approximate the quantities Ωdev j (p) ≈ n Ωdev j (p ) as constant for step 2 of our method; i.e., we locally approximate SAFT (n) Ωsim (p) + Ωdev j (p) ≈ Ω j j (p )

■

*Phone: +49 (0)711 685 66579; +49 (0)711 685 66103. Fax: +49 (0)711 685 66140. E-mail: [email protected].

(14)

Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS 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.

■

REFERENCES

(1) Martin, M. G.; Siepmann, J. I. Transferable Potentials for Phase Equilibria. 1. United-Atom Description of n-Alkanes. J. Phys. Chem. B 1998, 102, 2569−2577. (2) Jorgensen, W. L.; Madura, J. D.; Swenson, C. J. Optimized Intermolecular Potential Functions for Liquid Hydrocarbons. J. Am. Chem. 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. (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.: Condens. 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, 2-methylhexane, 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

(15)

Here, 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 n-alkanes in the objective function, so that

p ̂ ∈ IR8 (also with 8 values of ϕσ,k and ϕϵ,k). The 8 × 8 matrix [∇ϕp̂]T is of rank 8. In fact, it would become a diagonal matrix when the factor ϕσ,k in eq 10 was absorbed into ϕϵ,k. Consequently, from eq 15 we have ∇p̂ f dev = 0. A second application of the product rule of differentiation leads to ∇p̂ f dev = [∇p̂p]T ·∇p f dev

AUTHOR INFORMATION

Corresponding Author

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). Step 1 of Figure 14 minimizes squared deviations Ωdev = j Ωsim − ΩSAFT . The parameters of the SAFT equation of state j j entering ΩSAFT (m̂ k, σ̂k, ϵ̂k) are functions σ̂k(p, ϕσ,k) and ϵ̂k(p, j ϕϵ,k), whereas m̂ k is 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 ∇ϕ f dev = [∇ϕ p ̂]T ·∇p̂ f dev

ASSOCIATED CONTENT

(16)

For the alkanes [∇p̂p] is an 8 × 4 matrix 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 T

11705


Article

The Journal of Physical Chemistry B 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. (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 Equilib. 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) Bhatnagar, 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. Chem. 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. Chem. Soc. 1985, 107, 1489−1496. (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. Comput. Chem. 1990, 11 (8), 958−971. (24) Jorgensen, W. L.; Tirado-Rives, J. The OPLS Force Field for Proteins. Energy Mini- mization for Crystals of Cyclic Peptides and Crambin. J. Am. Chem. Soc. 1988, 110 (6), 1657−1666. (25) Nath, S.; Escobedo, 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 alpha-Olefins. 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 AlcoholAlkane 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. (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. Simul. 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. 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, 118 (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) 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. J. Phys. Chem. 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. Comput. Chem. 1999, 20, 1009. (48) Hülsmann, M.; Köddermann, T.; Vrabec, J.; Reith, D. ROW: A Gradient-based Optimization Workflow for Automated Development of Molecular Models. Comput. Phys. Commun. 2010, 181, 499−513. (49) van 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; Oxford University Press: Oxford, U.K., 1984. 11706


Article

The Journal of Physical Chemistry B (52) Lorentz, H. A. Ü ber die Anwendung des Satzes vom Virial in der kinetischen Theorie der Gase. Annalen der Physik. 1881, 248, 127−136. (53) Berthelot, D. C. R. Sur ur le mélange des gaz. Hebd. Seanc. Acad. Sci. (Paris) 1898, 126, 1703−1855. (54) van der 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; Clarendon Press: Oxford, U.K., 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.; Müller, M. Calculation of Free Energy through Successive Umbrella Sampling. J. Chem. Phys. 2004, 120, 10925− 10930. (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, 4533−4541. (61) Ferrenberg, A. M.; Swendsen, R. H. New Monte Carlo Technique for Studying Phase Transitions. Phys. Rev. Lett. 1988, 61, 2635. (62) Ferrenberg, A. M.; Swendsen, R. H. Optimized Monte Carlo Data Analysis. Phys. 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) Shen, V. K.; Errington, J. R. Determination of fluid-phase behavior using transition-matrix Monte Carlo: Binary Lennard-Jones mixtures. J. Chem. Phys. 2005, 122, 064508. (65) Smit, B. Phase diagrams of Lennard Jones 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. (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 Equilib. 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 LennardJones plus Pointquadrupole fluid. Fluid Phase Equilib. 2001, 179, 339− 362. (71) Vrabec, J.; Stoll, J.; Hasse, H. A set of Molecular Models for Symmetric Quadrupolar Fluids. J. Phys. Chem. B 2001, 105, 12126− 12133. (72) Stoll, J.; Vrabec, J.; Hasse, H. A set of Molecular Models for Carbon Monoxide and Halogenated Hydrocarbons. J. Chem. Phys. 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-γ). J. Chem. Phys. 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 Seg- ments. Fluid Phase Equilib. 2008, 274, 85−104. (75) Papaioannou, V.; Lafitte, T.; Avendaño, C.; Adjiman, C. S.; Jackson, G.; Müller, E. A.; Galindo, A. Group Contribution Methodology based on the Statistical Associating Fluid Theory for Heteronuclear Molecules formed from Mie Segments. J. Chem. Phys. 2014, 140, 054107. (76) Avendano, C.; Lafitte, T.; Galindo, A.; Adjiman, C. S.; Jackson, G.; Müller, E. A. SAFT-γ Force Field for the Simulation of Molecular Fluids. 1. A Single-Site Coarse Grained Model of Carbon Dioxide. J. Phys. Chem. B 2011, 115, 11154−11169. (77) Avendaño, C.; Lafitte, T.; Adjiman, C. S.; Galindo, A.; Müller, E. A.; Jackson, G. SAFT-γ Force Field for the Simulation of Molecular Fluids: 2. Coarse-Grained Models of Greenhouse Gases, Refrigerants, and long Alkanes. J. Phys. Chem. B 2013, 117, 2717−2733. (78) Lafitte, T.; Avendaño, C.; Papaioannou, V.; Galindo, A.; Adjiman, C. S.; Jackson, G.; Müller, E. A. SAFT-γ Force Field for the Simulation of Molecular Fluids: 3. Coarse-Grained Models of Benzene and Hetero-Group Models of n-Decylbenzene. Mol. Phys. 2012, 110, 1189−1203. (79) Sauer, E.; Stavrou, M.; Gross, J. Comparison between a Homoand 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. (80) NIST Chemistry WebBook; 2013; http://webbook.nist.gov/ chemistry. (81) 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. (82) Peng, D. Y.; Robinson, D. B. A New Two-Constant Equation of State. Ind. Eng. Chem. Fundam. 1976, 15, 59−64. (83) 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. J. Am. Chem. Soc. 1999, 121, 3407−3413. (84) Vrabec, J.; Stoll, J.; Hasse, H. A set of Molecular Models for Symmetric Quadrupolar Fluids. J. Phys. Chem. B 2001, 105, 12126− 12133. (85) 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. (86) Sage, B. H.; Hicks, B. L.; Lacey, W. N. Phase Equilibria in Hydrocarbons. Ind. Eng. Chem. 1940, 32 (8), 1085−1092. (87) Potoff, J. J.; Kamath, G. Mie Potentials for Phase Equilibria: Application to Alkenes. J. Chem. Eng. Data 2014, 59, 3144−3150. (88) 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. (89) Smith, B. D.; Srivastava, R. Thermodynamic Data for Pure Compounds: Part A Hydrocarbons and Ketones; Elsevier: New York, 1986

11707


Transferable Anisotropic United-Atom Force Field Based on the Mie

Recommend Documents