SAFT-γ Force Field for the Simulation of Molecular Fluids. 5. Hetero-Group Coarse-Grained Models of Linear Alkanes and the Importance of Intra-molecular Interactions

The SAFT-γ Mie group-contribution equation of state [Papaioannou et al., J. Chem. Phys., 140, 054107 (2014)] is used to develop a transferable coarse-grained (CG) force-field suitable for the molecular simulation of linear alkanes. A hetero-group model is fashioned at the resolution of three carbon atoms per bead in which different Mie (generalized Lennard-Jones) interactions are used to characterize the terminal (CH3CH2-CH2-) and middle (-CH2-CH2-CH2-) beads. The force field is developed by combining the SAFT-γ CG top-down approach [Avendaño et al., J. Phys. Chem. B, 115, 11154 (2011)] using experimental phase-equilibrium data for n-alkanes ranging from nnonane to n-pentadecane to parametrize the inter-molecular (non-bonded) bead-bead interactions, and a bottom-up approach relying on simulations with the higher resolution TraPPE united-atom (UA) model [Martin and Siepmann, J. Phys. Chem. B, 102, 2569 (1998)] to establish the intra-molecular (bonded) interactions. The transferability of the SAFT-γ CG model is assessed from a detailed examination of the properties of linear alkanes ranging from n-hexane (n-C6H14) to n-octadecane (n-C18H38), including an additional evaluation of the reliability of the description for longer chains such as nhexacontane (n-C60H122) and a prototypical linear polyethylene of moderate molecular weight (n-C900H1802). A variety of structural, thermodynamic, and transport properties are examined, including the pair distribution functions, vapour-liquid equilibria, interfacial tension, viscosity, and diffusivity. Particular focus is placed on the impact of incorporating intra-molecular interactions on the accuracy, transferability and representability of the CG model. The novel SAFT-γ CG force field is shown to provide a reliable description of the thermophysical properties of the n-alkanes, in most cases at a level comparable to the that obtained with higher resolution models.

interfacial tension, viscosity, and diffusivity. Particular focus is placed on the impact of incorporating intra-molecular interactions on the accuracy, transferability and representability of the CG model. The novel SAFT-γ CG force field is shown to provide a reliable description of the thermophysical properties of the n-alkanes, in most cases at a level comparable to the that obtained with higher resolution models.

Introduction
Molecular-dynamics and Monte Carlo computer simulation methods are very dependable techniques for studies of the structural and thermophysical and properties of real molecu-lar systems provided that an appropriate force field is available to describe the interactions between the particles. The structural simplicity of alkanes coupled with their practical significance have lead to large body of simulation studies and the development of numerous reliable force fields [1][2][3][4][5] . Linear alkanes constitute a convenient starting point for the development of force fields of more complex molecules such as macromolecules, polymers, and other compounds containing alkyl tails. As a result the development of models which are transferable for these classes of system over wide ranges of thermodynamic conditions and molecular weight continue to be a topic of significant interest. Molecular models of alkanes traditionally invoke all-atom (AA) 1,2 , united-atom (UA) 3,4,6-9 , or anisotropic united-atom (AUA) 5,10 representations. The AA level of resolution is clearly the more detailed and realistic, while the use of UA and AUA models reduces the number of interaction sites being considered resulting in an advantageous reduction in the overall computational overhead of the simulation 4 . Popular and reliable UA force fields are now available for the representation of fluid-phase equilibria of hydrocarbons including, for example, the SKS 6 , TraPPE 4,7 , and NERD 8 models, where the Lennard-Jones (12-6) potential is employed to treat non-bonded interactions between the UA beads. The Mie (n-6), generalized Lennard-Jones, potential has also been used to represent the the interactions between the UA beads in alkanes and perfluoroalkanes 9 ; a variation of the repulsive exponent (and range) of the potential allows for an improved description of the coexistence propeties.
Force fields at the AA and UA level of resolution are well suited for the simulation of the thermophysical properties of systems comprising a few thousand small molecules, but are not practical in studies involving complex phenomena such as self-assembly and micro phase separation, or macromolecular particles that can contain hundreds of thousands of atoms and span the micrometer scale. Lower-resolution coarse-grained (CG) models, where the interaction sites (beads) are taken to incorporate larger numbers of atoms, offer a significant advantage for systems characterized by large length and time scales. The use of CG force fields is becoming increasingly popular in bridging the gap between the atomistic scale and the mesoscale [11][12][13][14][15][16] .
Within a CG formalism the interaction site comprising a group of atoms or functional groups is represented with an effective potential designed to retain an accurate description of the target properties of interest. CG force fields are computationally more efficient, allowing one to explore phenomena involving longer timescales over larger length scales 11 .
However, this computational gain is achieved at the expense of a loss in chemical resolution in the molecular model, loss of transferability over thermodynamic states and/or for different systems, and representability 17,18 , whereby the properties of the system of interest cannot all be represented simultaneously.
The traditional procedure for the development of CG models is the so-called bottom-up approach whereby, once the level of resolution of the coarse graining (the number of atoms assigned per bead) is defined, information obtained at the more detailed molecular level is mapped at the less-detailed CG level 19 . This type of bottom-up approach usually relies on knowledge of the structural and thermodynamic properties obtained from simulations of high-resolution AA or UA models, which serve as a reference target for the development of the coarser models 20 . Key examples of bottom-up approaches include iterative Boltzmann inversion 21,22 and force-matching methods 19,[23][24][25] . The former is based on matching pair correlation functions at both the atomistic and CG levels, while the latter relies on minimizing the differences between the forces acting on the CG sites and the corresponding forces on the reference atomic system.
The level of resolution of the coarse graining will have an effect on the accuracy of the description of a given property, as will the parameters characterizing the bonded and nonbonded interactions of the CG force field. The influence of different levels of resolution in the CG procedure has been recently assessed for different mapping schemes of n-alkanes (with ≥ C 12 ) 26 . Taking n-dodecane as benchmark, one finds that as the degree of coarse gaining is increased (corresponding to models of progressively poorer resolution), the description of the thermodynamic properties deteriorates. It is known that CG models developed using bottom-up approaches can suffer from a lack of representability which prevents one from capturing properties across different thermodynamic states, in part due to the fact that typically only a single thermodynamic state is used in the mapping. An improvement in the representability of the CG model can be obtained by using a multi-state iterative Boltzmann inversion scheme rather than the single-state method as has been demonstrated for the structural properties of n-dodecane 27 .
In top-down CG methods the parameters characterizing the CG force field are instead determined by direct optimization of the description of experimental observables. This type of approach has been used very successfully in the development of the MARTINI force field 28,29 , where the partitioning free energies between polar and apolar phases of a variety of chemical species are employed in the parametrization. A resolution of four heavy atoms per bead (referred to as a 4:1 mapping) is typically employed for the MARTINI CG force field, and the standard Lennard-Jones (12-6) potential is taken to represent the form of the non-bonded interaction between all of the beads. The bonded intra-molecular interactions of aliphatic hydrocarbon chains are obtained from an analysis of the angular distributions and configurational entropies of pure hydrocarbons in the liquid state using the UA GROMOS force field 30 at the 4:1 level of resolution. A related parameterization for the n-alkanes at the 2:1 mapping level of resolution has recently been developed with the particle-swarm optimization technique 31 , where a good description of the heat of vaporization, the surface tension and the diffusion coefficient is demonstrated. The MARTINI force field is finding widespread use in simulation studies of biomolecular systems as a result of the group-contribution spirit of the approach which allows one to construct arbitrary molecules with off-the-shelf CG interaction parameters.
In this context it is important to mention the popular CG models of n-alkanes presented by Klein and co-workers 32,33 which are also developed with both top-down and bottom-up approaches. In this case the non-bonded parameters are determined using target experimental data for the liquid density and surface tension, while the parametrization of the bonded interactions is carried out using simulations at the atomistic scale. Chain molecules of various lengths were considered in the development of the force field but typically only at a single temperature. Of particular relevance to our current work is their use of the Mie (λ rλ a ) potential to represent the non-bonded interactions between the CG beads. By removing the constraint of a fixed repulsive exponent at the Lennard-Jones value of λ r = 12, one can obtain an improved description of the macroscopic observables of interest. In the case of the n-alkanes a Mie (9-6) potential was employed to parametrized the CG interactions 32 , and subsequently used to represent the aliphatic tail of phospholipid molecules 33 .
A similar approach was later followed by Maerzke and Siepmann 34 to develop a TraPPE-CG model for the n-alkanes in which the intra-molecular parameters were obtained using an iterative Boltzmann inversion procedure, and the inter-molecular interactions between CG beads, modelled with a potential of the Mie (λ r -6) form, were parametrized using vapourliquid coexistence data obtained from simulations of the corresponding TraPPE-UA model. The resulting force field was shown to provide an accurate description of the vapour-liquid equilibria for pure n-alkanes ranging from n-hexane to n-triacontane and for binary mixtures such as n-hexane + n-hexatriacontane.
These top-down approaches deliver CG models with much better transferability, but it is notable that the parameterization of the force field requires a number of computationally expensive molecular simulations to be carried out for intermediate systems characterized by suboptimal unbonded parameters, until the final set of parameters is obtained when the convergence criterion of the optimisation algorithm is met. A different approach based on the use of a free-energy function that explicitly takes into account the separate energetic and entropic contributions of the interactions has been proposed in order to improve the transferability of the models, using the linear n-alkane series as benchmark 35 .
In addressing the challenge of developing CG models that are robust in terms of their transferability and representability, the use of an molecular-based equation of state (EOS) has been shown to present an distinct advantage 36,37 . Providing the EOS yields a reliable and accurate description of the macroscopic properties of the underlying molecular model used in the development of the CG representation, the algebraic formalism can be used to determine the non-bonded parameters of the force-field with minimal computational effort, allowing one to consider large sets of target experimental data simultaneously in the estimation of the model parameters. The statistical associating fluid theory (SAFT) 38-42 framework provides a family of algebraic EOSs with a formal link between a detail molecular model and the thermodynamic properties of fluids of associating chain molecules. The PC-SAFT 43 EOS was used early on 44 to determine the non-bonded force-field parameters of n-alkanes CG beads by tuning the model parameters to reproduce experimental data for the liquid density, enthalpy of vaporisation, and vapour pressure. The lack of a direct correspondence between the PC-SAFT model and the underlying force field required the use of an empirical iterative simulation procedure to refine the parameters of the model in order to provide simulated properties which are in sufficient agreement with the target experimental data. A reliable force field can nevertheless be obtained in this manner providing an accurate description of the coexistence densities, vapour pressure, and caloric properties of a number of n-alkanes over a range of temperatures within just three simulation iterations. It should however be noted that because of the homonuclear nature of the underlying PC-SAFT model, different force-field parameters between the CG beads were obtained for each of the n-alkanes studied.
The introduction of the Mie potential to represent the basis of the inter-molecular interactions between the molecular segments and the extension of the high-temperature perturbation expansion to third-order in the SAFT-VR Mie EOS 45 provides a description of the thermodynamic properties of chain molecules formed from Mie beads which is of comparable accuracy to computer simulations for the underlying model. The use of an independent specification of the exponents describing the repulsive and attractive contribution of the Mie potential, allows for a versatile description of the inter-molecular interactions leading to a reliable description of the thermophysical properties, as demonstrated in a number of studies 9,33,34,36,[45][46][47][48] . The recasting of the SAFT-VR Mie formalism as the SAFT-γ Mie group-contribution approach 49 based on a heteronuclear model of molecules comprising Mie beads of different types (corresponding to different chemical functionality) has proved to be equally successful in terms of its applicability and accuracy in the description of the thermodynamic properties of a wide variety of systems and provides a framework to consider species with distinct chemical moieties explicitly 50 . The development of CG models using a top-down methodology with the SAFT-γ Mie EoS, hereafter referred to as SAFT-γ CG force fields, has been presented in a series of contributions for a variety of molecular fluids including one-site models of carbon dioxide 36,51,52 , homonuclear CG bead models of n-alkanes 47,53 , greenhouse gases and refrigerants 47 , natural gas condensates 54 , water and aqueous mixtures [55][56][57][58] , and heteronuclear models of aromatic compounds 46 , polystyrene polymers 59 , perfluoroalkylalkanes 60 , and amphiphilic systems comprising nonionic 61 , light-switching 62 , and super spreading surfactants 63 . The development and use of the SAFT-γ CG force fields in molecular simulation is reviewed in ref. 37 . A generalized methodology for the development of force field parameters for simple pure fluids based on a corresponding states approach has also been described 64,65 .
The development of SAFT-γ CG force fields for n-decane (n-C 10 H 22 ) and n-eicosane (n-C 20 H 42 ) represented as homonuclear models with beads of equivalent size presented in ref. 47 is of particular relevance to our current work. More specifically, the n-alkanes were modelled as chains of tangentially bonded identical spherical beads (three beads in the case of n-decane, and six beads in the case of n-eicosane), and, importantly, the chains were treated as fully-flexible; i.e., no intra-molecular bond bending or torsional interactions were considered. Similarly, a homonuclear fully-flexible chain model was used to represent solutions of polystyrene in n-alkanes 59 . The first-order thermodynamic perturbation theory (TPT1) of Werthiem 66-71 which is the basis for the representation of the chain contribution in SAFT formalism, does not explicitly take into account the configurational structure and flexibility (or rigidity) of chain molecules. For relatively simple fluids, most of the macroscopic thermodynamic properties are quite insensitive to the precise description of the microscopic molecular structure. This approximate treatment tends to break down, however, for more complex molecules such as amphiphiles or polymers, as is particularly apparent in the case of transport and interfacial properties.
The SAFT-γ CG force fields developed for fully flexible models of n-alkanes 47 provide a good representation of the experimental fluid-phase equilibria (vapour pressure and saturated density) and vapour-liquid interfacial tension over a broad range of thermodynamic conditions, with larger deviations notable for the description of the vapour pressure at lower temperature, especially for the longer chains. An empirical rescaling of the inter-molecular parameters can be applied a posteriori to improve the description of the saturation properties.
It is, however, generally accepted that even in the case of relatively simple molecules such as n-alkane, a high-fidelity representation of force field incorporating both the inter-molecular and intra-molecular interactions (bond-stretching, bending, and torsional contributions) is required to reliably represent the physical properties of the system 3-5 .
In our current work we present the SAFT-γ CG force-field for n-alkanes based on heteronuclear beads at a resolution of three heavy atoms per bead (3:1 mapping), incorporating bonded as well as non-bonded interactions. We follow a a bottom-up approach to determine the characteristic constants of the bonded potential from UA simulations, and a top-down approach using the SAFT-γ Mie EOS to characterise the parameters of the non-bonded interaction using experimental fluid-phase coexistence data (saturated density and vapour pressure). We highlight the importance of the implementation of a hetero-segmented model in the development of a transferable force-field for linear alkanes. A detailed analysis of the impact of the inclusion of the bonded potential is presented. This contribution adds to our group's body of work 36,46,47,55,[58][59][60] on the development of the SAFT-γ coarse-grained force field for the simulation of molecular fluids.
The paper is organized as follows: Details of the CG models are described in the next section, where the methodologies used to determine the parameters for both inter-and intra-molecular interactions are discussed. In Section 4 we present the results of molecular-dynamics simulations for the structural, thermodynamic, and transport properties of a nalkanes of varying chain length, including polyethylene chains of 300 CG beads (corresponding to a molecular weight of 12,602 g mol −1 ), making appropriate comparisons with experimental data and the corresponding description with atomistic (UA) models. Our concluding remarks are summarised in Section 5. Details of the simulation methodology are provided in section 3.

Inter-molecular interactions
We consider heteronuclear (hetero-group) segmented chain molecules formed from tangential spherical beads of different types. The interaction between two beads k and l is represented with the Mie potential, where r kl is the centre-centre inter-bead distance, ε kl the well depth of the potential, σ kl the diameter of the bead, λ r kl the exponent controlling the hardness/softness of the repulsive contribution, and λ a kl the exponent controlling the range of the attractive contribution of the interaction. The constant C kl is defined as to ensure that the energy at the minimum of the potential is −ε kl irrespective of the values of the repulsive and attractive exponents.
As depicted in Figure 1 in the case of n-nonane, a 3:1 level of CG resolution is implemented, in which three consecutive backbone carbon atoms (with the corresponding bonded hydrogen atoms) are modelled as one CG bead. In our heteronuclear model two types of CG beads are used corresponding to either terminal (T) beads, representing ( The non-bonded inter-molecular parameters of the force field are estimated using the SAFT-γ Mie group contribution EOS 49 , carrying out a minimization of an objective function F obj based on the unweighed square residuals of experimental data (exp) and the corresponding SAFT-γ Mie EoS description (cal) for the chosen properties. As is common practice in engineering applications we use the vapour (saturation) pressure and saturated-liquid density data over a wide range of temperatures as target properties. The objective function is given as where N Pv and N ρ l are the number of experimental data points for the vapour pressure P v and the saturated-liquid density ρ l of compound i, respectively. The vector of SAFT parameters at the corresponding temperature T of data point q is represented by α = [σ kl , ε kl , λ r kl , λ a kl ].
The attractive range is fixed as λ a kl = 6 in all cases to reflect the London-type dispersion interaction characteristic of these non-polar systems.
The inter-molecular parameters of our heteronuclear model are based on the analysis of data for several n-alkanes, so that a unique set of parameters will, in principle, be transferable to other n-alkanes. While data for an arbitrary number of n-alkanes could be used, here we consider n-C 9 H 20 , n-C 12 H 26 and n-C 15 H 32 as the target compounds in the parameter estimation. In this way one can ensure that the shortest chain (n-hexane n-C 6 H 14 , which comprised only two terminal CG beads) and longer n-alkanes, for which the experimental data is less reliable, do not bias the resulting model.
The estimated SAFT-γ CG inter-molecular parameters are presented in Table 1 Table 2. The unlike parameters for the diameter σ TM and range λ r TM of the model are determined using combining rules 49 : and The deviation of the description obtained with the SAFT-γ Mie EoS (cal) using the parameters presented in Table 1 in comparison to experimental (exp) data is presented in Table 2, in terms of the percentage average absolute deviation %AAD, which is determined Table 1: SAFT-γ CG Mie inter-molecular force-field parameters for beads of type T (terminal, CH 3 -CH 2 -CH 2 -groups) and type M (middle, -CH 2 -CH 2 -CH 2groups) for n-alkanes: σ kl is the bead diameter, ε kl the well depth of the potential, and λ r kl the repulsive exponent of the potential. The attractive exponent λ a kl = 6 is fixed to the London value in all cases. The unlike parameters σ TM and λ r TM are determined using the combining rules defined in Equations 4 and 5, while ε TM is estimated by minimization of Eq. 3 using the SAFT-γ Mie EoS 49 . k B corresponds to the Boltzmann constant. where N X j is the number of data points of a given property X j . Table 2: Percentage average absolute deviations (%AAD) for the vapour pressure P v (T ) and saturated-liquid density ρ l (T ) obtained with the SAFT-γ Mie equation of state 49 and the CG model of Table 1 for selected n-alkanes. The asterisk indicates that experimental data for n-C 6 H 14 and n-C 18 H 38 were not included in the parameter estimation; the calculations in this case are predictive. The transferability of our SAFT-γ Mie CG model is further assessed for the prediction of properties of compounds which were not included in characterising the force-field parameters.
The fluid-phase equilibria is determined with the SAFT-γ Mie EoS for n-C 6 H 14 and n-C 18 H 38 using the optimized parameters given in Table 1 and the corresponding %AADs are also provided in Table 2. For these two compounds the deviations are found to be slightly larger, especially in the case of the predicted vapour pressures. To a certain extent, the tangent-bead model developed here limits our ability to provide a very accurate description of the vapour pressure for a large range of linear alkanes 74 , although it is important to note also that the vapour pressures of long n-alkanes are very low and that a more appropriate measure of the performance of the would be in absolute and not relative terms; we use such a measure for comparison of vapour pressure data later in the current work.

Intra-molecular interactions
The intra-molecular interactions are described with harmonic potentials for both bond stretching and angle bending: where k bond and k angle are the stretching (bond) and bending (angle) constants, respectively, r 0 and θ 0 are the equilibrium bond length and bond angle, respectively, and k B is the Boltzmann constant. The dihedral contribution to the intra-molecular potential is not been included in the force field, as magnitude of the torsional term is found to be lower than k B T .
Non-bonded interactions between pairs of beads separated by more than three bonds are considered.
The intra-molecular potential parameters are obtained by carrying out molecular dynamics (MD) simulations using the TraPPE-UA force field 4 for n-C 6 H 14 , n-C 9 H 20 , and n-C 12 H 26 at a number of thermodynamic conditions (cf. Table 3). We note here that simulations using the NERD-UA force field 75 would lead to very similar results as the intramolecular parameters are identical to those employed in the TraPPE-UA model. A preliminary evaluation of the intra-molecular parameters of the n-alkanes for gas and liquid states at two states (400 and 600 K) lead to values of k bond /k B = 7373 K/Å 2 (r 0 = 3.50Å) and k angle /k B = 2124 K/rad 2 (θ 0 = 159.9 o ) when obtained at a temperature of 400 K, and k bond /k B = 7588 K/Å 2 (r 0 = 3.50Å) and k angle /k B = 2668 K/rad 2 (θ 0 = 157.6 o ) when averaged for states at 400 K and 600 K. These preliminary values of the intra-molecular interactions for the n-alkanes have already been used in some studies with the SAFT-γ Mie CG force-fields of chain molecules 58,60 .
A detailed examination of the intra-molecular parameters is provided over a range of temperatures for gas and liquid states densities (chosen to avoid any possibility of metastability) to improve the transferability of the resulting force filed. Subcritical gas and liquid states Table 3: Thermodynamic density (ρ) and temperature (T ) states and intramolecular parameters (bond distance r 0 , bond stretching constant k bond , angle θ 0 , and bond bending constant k angle ) of the SAFT-γ CG force field for n-C 6 H 14 , n-C 9 H 20 , and n-C 12 H 26 . The TraPPE-UA force field is used as the reference model in the atomistic MD simulations and equation (8) is used to determine the intra-molecular parameters. The results of the parametrization of the intra-molecular parameters at each state are reported together with the overall unweighed average values of each parameter.  Normalized frequency, Normalized frequency, Figure 2: Representative intra-molecular probability distributions for the stretching bonds and bending angles obtained for n-C 12 H 26 at two thermodynamic states (one supercritical, and one subcritical). The description is obtained using the reference TraPPE-UA force field with a 3:1 level of CG mapping. The symbols represent the results obtained from N V T molecular-dynamics simulations, and the continuous curves to the corresponding fits using Eq. 8 with a weighting factor of w(r) = 1 for the bond distributions and w(θ) = sin(θ) for the angular distributions, respectively. and supercritical states for each system are used for the parameterization. The subcritical states presented in Table 3 correspond to temperatures 15% and 25% below the critical temperatures of each n-alkane with corresponding densities that are 15% higher (lower) than the experimental 72 saturated-liquid (vapour) densities to ensure single-phase states are considered. The supercritical states correspond to temperatures 20% above the reported experimental critical temperature for densities corresponding to the critical density.
The intra-molecular potentials are determined by defining effective CG beads, using the  Table 3 using weighted Gaussian distributions P (x) of the form: where x is either the distance x = r between adjacent beads or the the angle x = θ subtended by three consecutive beads, and the mean µ is µ = r 0 or µ = θ 0 , with s 2 the standard deviation. A weighting function w(x) is used in Eq. 8 to take into account the asymmetry of the distributions, with w(r) = 1 for the bond distributions and w(θ) = sin(θ) for the angular distributions, respectively. The spring constants in Eq. 7 are then obtained as k bond ∝ k B T ln P (r) and k angle ∝ k B T ln P (θ)/ sin(θ), respectively 76 .
The bond stretching and bending parameters of the SAFT-γ Mie CG force-field obtained from the UA trajectories are collected in Table 3  The overall averages of the bonding constants over the thermodynamic states given in Table 3 are used to represent the bonded potential of the SAFT-γ Mie CG force field for n-alkanes. As the intra-molecular potential is obtained over a broad range of thermodynamic conditions, one would expect to have a better transferability (and representability) of the properties of n-alkanes. It is important to note that the value of r 0 obtained in the parametrization of the bond-length distributions is not used in the representation of the CG model.
As we are using a tangent chain model within the SAFT-γ Mie description the corresponding bead diameters (cf. Table 1) are used instead.

SAFT-γ Mie CG force field for n-alkanes
The Bond stretching, bending, and torsional potentials are however usually included in simulations of molecular systems 23,29,33,34 and are known to be important in the determination of the properties of dense systems, particularly the structure. The impact of the inclusion of these intra-molecular interactions in our proposed force field is, therefore, of particular interest. With this in mind three SAFT-γ CG models (A, B, and C) are considered (cf. Table 4) to study a number of properties of n-alkane fluids with effective bond and bending constants for each of these models. Table 4: Summary of the effective bond and bending parameters used in our SAFT-γ CG force field models for n-alkanes. Models A, B, and C are described in terms of effective bond k bond,eff and bending k angle,eff constants taken as multiples of the unweighed average values reported in Table 3. The equilibrium bond length r 0 for all the models are taken from the diameters of the CG beads reported in Table 1.  Table 3.

Theoretical Methods
We

Thermodynamic properties
The vapour-liquid equilibrium properties of n-alkanes are determined using a direct coexistence method using MD-N V T simulations. The molecules are arranged in an orthorhombic box with a constant volume V = L x ×L y ×L z , such that L x = L y << L z . After equilibration the density profile along the z direction (normal to the interface) is used to determine the saturated densities of the vapour and liquid phases 84 . The vapour pressures are calculated in two ways: in the first method the saturated gas density from the direct coexistence simulations is used to obtain the vapour pressure via the ideal gas law; in the second method, a single phase MD-N V T simulation using the saturated gas density is carried out and the pressure is obtained through the virial route 85,86 .
The vapour-liquid interfacial tension γ is obtained from the diagonal components of the pressure tensor 87,88 as where P zz , P yy and P xx are the average values of the components of the main diagonal of the pressure tensor. The pre-factor 1/2 denotes the presence of two interfaces.
In the calculations using the TraPPE-UA force field, 7200 beads are used, while in the CG simulations 13800 beads are used, except for the simulations of n-C 60 H 122 and polyethylene where 11500 and 9000 CG beads are used, respectively, in simulations spanning up to 25 ns to allow relaxation of these long molecules. The simulations of these large molecular systems were carried out using the HOOMD-blue CUDA-enabled MD package to accelerate the calculations.

Transport properties
The diffusion coefficient D for n-alkanes ranging from n- where r(t) is the position of the particle at time t and r(0) the reference position.
The viscosity is obtained using non-equilibrium MD (NEMD) simulations in the N V T Larger shear rates were also applied to ensure that the system was not in the shear thickening regime. These shear rates gave a linear regime for the equilibrium viscosity to be calculated by extrapolating to a zero shear rate. These simulations were carried out using GROMACS,

Results and Discussion
In this Section the impact of including intra-molecular interactions on a number of struc-

Structural properties of n-alkanes
The centre-of-mass pair distribution functions g(r) of n-hexane, n-nonane, and n-dodecane obtained from N V T molecular dynamics simulations using SAFT-γ Mie CG force fields A, B, and C for selected thermodynamic states used during the parametrization of the models, are presented in Figure 3. The results are compared with the results obtained from simulations using the TraPPE-UA force field 4 .
We start by analysing the results for n-hexane. In our 3:1 CG model, n-hexane is treated as a dimer with two tangentially bonded terminal (T) beads, thus only models A and C are compared with the UA model. In terms of the contributions to the intra-molecular potential contributions, only bond stretching is considered; the difference between the results for models A and C for n-hexane is, therefore, only the stiffness of the bond potential. The pair distribution functions obtained for n-hexane are shown in Figure 3 Table 3.
the simulation results for models A and C are essentially indistinguishable, suggesting that a bond constant k bond /k B = 6666 K/Å 2 is large enough to keep the beads in our CG model tangential essentially corresponding to a rigid bond. For the lowest temperature (T = 381 K), both CG models exhibit two maxima (at 6.5Å and 11.0Å) close to those characterizing the TraPPE-UA force field. The peak at ∼6.5Å in the CG model is shaper compared to the UA model, due to the development of an additional peak at ∼4.9Å that corresponds to the contact distance between two CG beads. This additional peak disappears as the temperature is increased; for T = 432 K the peak at ∼6.5Å in both CG models is almost as broad as in the case of the UA model. In the of the supercritical state at T = 609 K, the pair correlation functions for both CG models and the UA are essentially indistinguishable.
The pair correlation functions for n-nonane and n-dodecane are shown in Figures 3(b,c).
At subcritical conditions (T = 446 K and T = 505 K), the curves obtained for n-nonane with the TraPPE-UA force field exhibit two shallow peaks located at ∼5.8Åand ∼8.8Å. Both of the SAFT-γ CG models B and C, which take into account the bending contribution to the intra-molecular potential, exhibit the same features as the UA force field although sharper peaks are seen for the CG models. Model A, however, exhibits a single peak (located at

Vapour-liquid equilibria of n-dodecane
The SAFT-γ Mie CG model for n-dodecane comprises four beads: two terminal (T) beads and two middle (M) beads, providing a suitable case of study to assess the influence the intra-molecular interactions in the description of the vapour-liquid equilibria. The saturated densities of the liquid (ρ l ) and vapour (ρ v ) phases, the vapour pressures P v , and the vapourliquid interfacial tensions γ are assessed in Figures 4(a-c). Molecular dynamics simulations are carried using the three SAFT-γ Mie CG models (A, B, and C) presented earlier, and the results are compared to corresponding data obtained from simulations with the TraPPE-UA force field, the SAFT-γ Mie EoS, and with the experimental data 72 . Overall, the three CG models provide a very good description of the experimental data, the TraPPE-UA data, and SAFT-γ Mie EoS. The calculated %AADs with respect to the experimental data for each of the models are presented in Table 5.
In the case of the saturated liquid density ρ l , small differences are observed for the three CG models, the results of the UA model, and the experimental data. As expected these differences are more evident near the critical point. It is encouraging to note that the order of magnitude of the error is comparable for all of the models, including the more detailed UA model, with the advantage of an increase in the speed of the calculations in the case of the CG models due to the reduction in the number of beads used to represent the molecule. Very good agreement is apparent Figure 4 between the simulated results with the CG models and the predictions using the SAFT-γ Mie EoS, confirming the reliability of the theory employed to derive the molecular parameters.
In the case of the vapour pressure, the models considered provide a good description over the entire range of temperatures (from the triple point to the critical point). From a detailed examination of the calculated errors, we note that this is a very sensitive property to calculate via molecular simulation, while the %AAD obtained is very small for the calculations carried out with the SAFT-γ Mie EoS (cf.  Table 5 are given in absolute rather than relative terms for a more appropriate comparison. In the case of the SAFT-γ Mie CG models, the smallest deviations are found for Models B and C, with a slightly larger error apparent for model A (in which a fully flexible CG chain model is used). All of the models considered provide a good description of the vapour-liquid interfacial tension tension (cf. Figure 4 (c)) over the entire temperature range, with small %AAD found in each case. It is encouraging to find that the proposed SAFT-γ Mie CG force field provide good description of the surface tension, with %AAD comparable to that obtained with the TraPPE-UA model; this is an encouraging in terms of the representability of the CG models.
Overall, it is difficult to distinguish between the adequacy of the three SAFT-γ CG models in the reproduction of the fluid-phase coexistence and interfacial tension of n-dodecane.
However, the comparisons provided in Figure 3 clearly highlight the advantage of using a model with angle constraints (CG models B and C) in terms capturing the structural properties of n-alkane fluids. In the calculation of bulk and interfacial tension properties little difference is found between model B and C; a more accurate saturated liquid density and pressure are obtained with model C, but a more accurate interfacial tension is obtained with model B (cf. Table 5). The only difference between these two force fields is in the bond constant, which in the case of model B corresponds to a effectively rigid bond. As a more physical representation of a molecule corresponds to a harmonic bond rather than to a rigid one, model C is adopted at this point to study the transferability of SAFT-γ Mie CG force  Table 5: Percentage absolute average deviation %AAD obtained with the SAFTγ Mie EoS, with molecular dynamics simulations using the SAFT-γ Mie CG models A, B, and C, and with simulations using the TraPPE-UA for the saturated liquid density ρ l , and vapour-liquid interfacial tension γ, and absolute deviations for the vapour pressure P v for n-dodecane and n-hexacontane. The experimental data are from references 72,91 . The deviations in pressure are calculated as simple averages from Dev = 1 | as the values for this property are very small for the systems of interest.

Transferability of the SAFT-γ Mie CG force field for n-alkanes
Following from the previous section, we proceed to test the transferability of the SAFT-γ Mie CG model C to describe the fluid-phase equilibrium properties of a range of n-alkanes vapour pressures as seen in Figure 5(a,b). The results obtained for γ using the SAFT-γ Mie CG model C are shown in Figure 5(c). Good agreement is seen in the case of the longer n-alkanes, even at low temperatures, though we note that the predictions for n-C 6 H 14 exhibit slightly larger deviations for the interfacial tension at all temperatures. This larger deviation is justified given that only data for the longer n-alkanes, which contain terminal and middle beads, are used in the determination of the model parameters. In other work 76 the terminal-terminal non-bonded bead-bead interactions are determined using target data for n-hexane, with the other parameters (those involving the middle-middle and the terminalmiddle interactions) determined using longer n-alkanes.
Overall we find the agreement between the molecular dynamics simulations carried out with our SAFT-γ MIE CG force field (model C) to be very satisfactory when compared with experiment, from which we conclude that that the proposed SAFT-γ CG force field can be successfully applied to a range of n-alkanes, over a wide range of temperatures, pressure,s and densities. The transferability of our force field is assessed for longer n-alkanes in the following section.

Predictions with SAFT-γ Mie CG force field for long-chain
n-alkanes

n-Hexacontane
We carry out molecular dynamics simulations for n-C 60 H 122 using our SAFT-γ Mie CG model C. The results for the fluid-phase coexistence properties and the vapour-liquid interfacial tension as a function of temperature are presented in Figure 6, with appropriate comparisons to correlated pseudoexperimental data from the DECHEMA database 91 .Also included in the comparison are simulation results using the TraPPE-UA force field 92 , and the theoretical predictions obtained using SAFT-γ Mie EoS 49 with the same non-bonded parameters. As can be seen from the figure the agreement between the simualtion data obtained with the SAFTγ Mie CG model C and the experimental data is satisfactory for the properties considered.
This is very gratifying considering that the model was developed for much shorter linear n-alkanes, and it further confirms the transferability of the proposed force field. A summary of the %AAD and deviations using the different models are reported in Table 5. Care should be taken, however, in interpreting these results, as n-alkanes are expected to be thermally unstable at temperatures above 650 K and the reported DECHEMA data 91 are extrapolated correlations.
The coexistence envelope obtained from the CG simulations and the theoretical calculations are very close, and both are seen to slightly underpredict (more so at the higher temperatures) the experimental saturated liquid density. The data reported with the TraPPE-UA force field is also seen to underpredict the experimental saturated liquid density as observed in Figure 6(a), although to a lesser extent. One can also see from Figure 6 (b) that the SAFT-γ Mie EoS calculations underpredict the vapour pressure over the temperature range considered. The enthalpy of vaporization (which is related to the slope of the curve) is however quite well represented; the proposed CG force field leads to a slightly larger slope of the Clausius-Clapeyron (ln P vs 1/T ) curve, overpredicting the pressure at high temperature and underpredicting it at low temperatures. The TraPPE-UA force field provides the best agreement at high temperature although it markedly overpredicts the vapour pressure at low temperatures.
As for the shorter n-alkanes considered ealier, we note that the vapour pressure is a notably difficult property to calculate accurately from molecular simulations in the case of very non-volatile substances; because the vapour pressure is very low, very few molecules are present in the vapour phase, which leads to poor statistics. We conclude that the overall agreement for all of the models considered is fair. The interfacial tension (Figure 6(c)) obtained with the SAFT-γ Mie CG force field is in satisfactory agreement with experiment, although larger deviations are found at lower temperatures than with the TraPPE-AU model, which is very accurate for this property.

Polyethylene
The benefits and importance of developing accurate CG models for polymers has been highlighted in a recent paper 93

SAFT-γ Mie CG force field
Our comparison of SAFT-γ Mie CG force fields (models A, B, and C) presented in section B (cf. Table 4) did not highlight significant differences in the resulting description of the fluidphase behaviour and interfacial tensions of the n-alkanes considered, although it is apparent, and perhaps unsurprising, that the microscopic structure of the fluid is very sensitive to the intra-molecular interactions. Here, we further consider transport properties in order to assess the importance of the bonded interactions in the representability of the SAFT-γ Mie CG force field. The viscosity and diffusivity of n-C 12 H 26 and n-C 24 H 50 (chosen as representative molecules that are small enough for ease of computation, but long enough to study the impact of the bonded potential) with the expectation that these transport properties will be subject to the details of the intra-molecular parameters.

Viscosity
The three SAFT-γ Mie CG models (cf. Table 4) are used to calculate the viscosity of ndodecane at T = 400 K and P = 1 MPa. The results of the molecular-dynamics simulations are presented in Table 6. It is immediately apparent that the fully flexible model A leads to an underprediction (of 17.5%) of the experimental data, whereas models B and C, in which a bonded potential is considered, lead to values considerably closer to the experimental values 72 (7.5% error for the case of model B and 2.5% with model C). This is encouraging considering the coarse-grain nature of the force field, and the fact that it was not parametrized for transport properties. It is, however, difficult to unequivocally assert which of models B and C is better in describing the viscosity of n-dodecane, keeping in mind that the experimental values are reported to carry an uncertainty of up to 5% 99 , and that the two models differ only in the stiffness of the bond potential.
We have also carried out calculations using model C with an additional torsional (dihedral) contribution, using a similar procedure to used in Section 2.2. The value obtained for the torsional constant of the n-alkanes is found to be less than k B T . The incorporation of the torsional contribution is therefore found to have minimal impact the description of the viscosity of n-dodecane (data not shown). One should however acknowledge that the small effect observed in the case of n-dodecane, may be more important in longer n-alkane chains and polymers and should not be ruled out in further investigations. At this point, given the slightly better performance of SAFT-γ Mie CG models C for the prediction of viscosity, we continue to use model C to compute the viscosity over a range of temperatures for both n-dodecane and n-tetracosane. The viscosities of n-C 12 H 26 and n-C 24 H 50 are plotted over a range of temperatures in Figure 8 . The corresponding experimental and simulation data presented are collected in Table 7. Overall we find that the proposed coarse-grain force field provides an accurate prediction of this transport property, confirming the representability of the SAFT-γ Mie CG model C. As can be seen from the figure, the viscosity of n-dodecane is predicted accurately over the entire temperature range considered. Although in the case of n-tetracosane the viscosities obtained with the CG model underpredict the experimental values at lower temperatures (which could be an indication of the onset of solidification), with a slight overprediction of the data at higher temperatures, overall the representation can be considered adequate. The average absolute deviation from the experimental values is less than 5% for ndodecane and 15% for n-tetracosane. The activation energy for viscous flow can be calculated from the expression η = η 0 exp(E/RT ) using the temperature-dependence of the data. Excellent agreement with experiment is found for n-C 12 H 26 , with values of (E/R) sim = 1670 K and (E/R) exp = 1721 K. In the case of n-C 24 H 50 a larger deviation is noted, mostly due to the under-prediction of the viscosity at the lowest temperature considered: an activation energy of (E/R) sim = 1640 K is obtained from the CG simulation data, while the experimental value is (E/R) exp = 2029 K.
The transport properties of n-alkanes using CG models based on fully-flexible chains of Lennard-Jones segments have been reported in a recent study by Galliero 48 . Larger errors were obtained for the viscosity for lower temperatures with a flexible model, and Galliero notes that an additional parameter to take into account the 'rigidity' of the molecule is needed to improve the accuracy of the description. It is unfortunately difficult to make a direct comparison with our model as different chain lengths are considered; we note, however, that the use of the Mie potential (instead of a fixed-range 12-6 Lennard-Jones form) and the inclusion of intra-molecular interactions in our SAFT-γ Mie CG force-field provides a force-field that is transferable over thermodynamic states and that performs well in terms of representing a range of thermophysical properties, including the viscosity.

Diffusivity
The diffusion coefficient is an important transport property that can be easily calculated from molecular-dynamics simulations; the theoretical prediction of the diffusivity of dense fluids is still lacking 48 . The diffusivity serves as a further test of the representability of the SAFT-γ Mie CG force field. The diffusion coefficient for n-C 12 H 26 is calculated first by carrying out N V T molecular-dynamics simulation at one state point (T = 303 K and 759.7 kg m −3 ) with the SAFT-γ Mie CG models A, B, and C (cf. Table 4), in order to assess the impact of the different contributions in the bonded potential, by comparison to experimental data 103 . The results are collected in Table 6 and shown in Figure 9. It is encouraging to find that our SAFT-γ force fields provide a reasonable description of diffusion coefficient even though transport properties are not considered in the parametrization of the intra-and inter-molecular interactions.
In more detail, we note the larger deviation from the experiment for model A. As in the case of viscosity, it is difficult to differentiate between the performance of models B and C. In this case, this is due in part to the large uncertainty in experimental values 100,101 .
Overall the findings allow us to conclude once again that there is indeed an advantage of incorporating the intra-molecular potential in the SAFT-γ CG force field in terms of the computation of transport properties. Moreover, given that model C incorporates a more realistic representation of the bond constant as well as providing a slightly better prediction of the transport properties studied (especially noticeable for the viscosity), we proceed to use model C to predict the diffusivity of a number of n-alkanes as a further challenge to the transferability and representability of the CG force field.
Diffusion coefficients for n-C 8 H 18 to n-C 16 H 34 ) are determined from N V T moleculardynamics simulations at 303 K for the densities reported in Table 8 using the SAFT-γ Mie CG force field Model C. The simulated values are compared to experimental data 103 and with the CG simulations reported by Nielsen et al. 33 in Figure 9. The experimental data presented here only spans a limited range of n-alkane chain lengths, but it suffices to reveal a trend with molecular weights at the selected temperature. As can be seen in the figure, we find good agreement between the simulation data obtained with the SAFT-γ Mie CG model C and the experimental data available. We note the larger (negative) slope in the Table 8: Diffusion coefficient D calculated via N V T molecular-dynamics simulations with the SAFT-γ Mie CG force field, model C, for a number of n-alkanes at T = 303 K. The density ρ used for each system is also given, corresponding to the resulting equilibrium density from N P T simualtions at a temperature of T = 303 K and pressure of P = 1 MPa. lead to an overprediction of the diffusion coefficient corresponding to much faster dynamics.
They also a 3:1 CG level of mapping, but use Lennard-Jones interactions to describe the bead-bead interactions instead of the harder Mie form incorporated in our SAFT-γ Mie CG force field. The choice of a more repulsive bead-bead interaction (λ r kl ∼ 16) allows us to obtain a better prediction of the transport properties. It is important to emphasize that this value of the repulsive exponent is estimated from vapour-pressure and saturated-liquid density equilibrium data alone in the parametrization of the inter-molecular potential; the transport or interfacial properties are not used in the development of the model.

Conclusions
A SAFT-γ Mie CG heteronuclear model for linear alkanes has been developed incorporating intra-molecular (bonded) and and inter-molecular (non-bonded) interactions with the aim of addressing the transferability, representability and robustness of the force field. We have angle contributions, which are obtained following a bottom-up approach using Trappe-UA models; we find a very small effect for the dihedral contribution, which is hence neglected.
The resulting force-field provides a very accurate description of the coexistence properties of the n-alkanes used in the parameterization (n-nonane, n-dodecane, and n-pentadecane) and is transferable in the prediction of the thermodynamic properties of n-alkanes not included in the parametrization, including long linear alkanes and polymers; in particular n-C 60 H 122 and a linear polyethylene of 900 carbons. The force field performs very well in terms in terms of representability. We report simulated values of the transport properties (viscosity and diffusivity) which are shown to be good agreement with experiment for a range of n-alkanes. The level of agreement we find in terms of the diffusivity is particularly encouraging.
A detailed assessment of the impact of the bonded potential in the SAFT-γ Mie CG force field is made by provided comparisons between three models: a pearl-necklace chain (fully flexible with rigid bonds), a semiflexible chain with rigid bonds, and a semiflexible chain with harmonic bonds. Interestingly, little difference between these three cases is found for the fluid-phase behaviour of n-dodecane. As expected, however, the inclusion of the intra-molecular bonded potential is critical in terms of reproducing the structural properties of the system. Furthermore, the fully flexible model is found to present larger deviations in the description of the transport properties simulated; a good prediction of the viscosity and diffusivity is achieved by incorporating the intra-molecular interactions.
In conclusion, a transferable CG force-field is presented for n-alkane fluids, which includes