Molecular Dynamics Simulation of Pure n-Alkanes and Their Mixtures

Article Cite This: J. Phys. Chem. B XXXX, XXX, XXX−XXX

pubs.acs.org/JPCB

Molecular Dynamics Simulation of Pure n‑Alkanes and Their Mixtures at Elevated Temperatures Using Atomistic and CoarseGrained Force Fields Konstantinos D. Papavasileiou,†,‡ Loukas D. Peristeras,† Andreas Bick,‡ and Ioannis G. Economou*,†,§

Downloaded via NOTTINGHAM TRENT UNIV on July 17, 2019 at 19:54:34 (UTC). See https://pubs.acs.org/sharingguidelines for options on how to legitimately share published articles.

†

Institute of Nanoscience and Nanotechnology, Molecular Thermodynamics and Modelling of Materials Laboratory, National Center for Scientific Research “Demokritos”, Aghia Paraskevi, Attikis, GR-15310 Athens, Greece ‡ Scienomics SARL, 16 rue de l’Arcade, 75008, Paris, France § Chemical Engineering Program, Texas A&M University at Qatar, Education City, P.O. Box 23874, Doha, Qatar S Supporting Information *

ABSTRACT: The properties of higher n-alkanes and their mixtures is a topic of significant interest for the oil and chemical industry. However, the experimental data at high temperatures are scarce. The present study focuses on simulating n-dodecane, n-octacosane, their binary mixture at a n-dodecane mole fraction of 0.3, and a model mixture of the commercially available hydrocarbon wax SX-70 to evaluate the performance of several force fields on the reproduction of properties such as liquid densities, surface tension, and viscosities. Molecular dynamics simulations over a broad temperature range from 323.15 to 573.15 K were employed in examining a broad set of atomistic molecular models assessed for the reproduction of experimental data. The well-established united atom TraPPE (TraPPE-UA) was compared against the all atom optimized potentials for liquid simulations (OPLS) reparametrization for long n-alkanes, L-OPLS, as well as Lipid14 and MARTINI force fields. All models qualitatively reproduce the temperature dependence of the aforementioned properties, but TraPPE-UA was found to reproduce liquid densities most accurately and consistently over the entire temperature range. TraPPE-UA and MARTINI were very successful in reproducing surface tensions, and L-OPLS was found to be the most accurate in reproducing the measured viscosities as compared to the other models. Our simulations show that these widely used force fields originating from the world of biomolecular simulations are suitable candidates in the study of n-alkane properties, both in the pure and mixture states.

1. INTRODUCTION Normal alkanes (n-alkanes) are compounds of great importance for the oil and chemical industries as they are abundant in multicomponent fuel, oil, and wax fluid mixtures,1 which are heavily utilized in the manufacturing of many consumer goods and cosmetics.2,3 Moreover, they have served as model compounds in gaining insight into the fundamental driving forces of biologically related mechanisms of protein folding and binding,4 given the aliphatic content of several amino acid side chains. n-Alkanes comprise the alkyl chains of several fatty acids and are vital in controlling lipid self-assembly and micelle formation in water.5,6 Finally, they are important components of “leaf waxes” and therefore find uses in environmentally related processes as biomarkers.7,8 The knowledge of vapor−liquid equilibria (VLE) properties including composition and density of each coexisting phase along with the related interfacial and rheological properties such as surface tension and viscosity, are central in controlling the individual processes involved for processing and production. Specifically, the vapor−liquid (VL) surface tension, γVL, of pure and mixed liquids is an important property for almost every separation process involving mass transfer between the vapor and liquid phases. Controlling and © XXXX American Chemical Society

predicting transport properties such as the viscosity of nalkanes and their mixtures is critical in rational product design and process optimization particularly in the petrochemical industry. These properties can be determined through both experimental and molecular simulation techniques. Experimental efforts have been dedicated on obtaining accurate data on these properties for decades. Experiments are often time consuming, expensive, and overall become a very tedious task, particularly as the n-alkane chain length increases and especially at conditions relevant to some industrial processes (i.e., at elevated temperatures and pressures), because of rapid thermal degradation above 650 K.9,10 One of the earliest reported works on surface tension is that of Jasper et al., who examined n-pentane to n-eicosane at temperatures between 273 and 353 K,11,12 while Doolittle and Peterson studied the liquid densities and viscosities of ten n-alkanes from n-pentane to n-tetrahexacontane at atmospheric pressure and temperatures in the range of 263−573 K.13 Several others have reported on these properties for pure nReceived: March 26, 2019 Revised: June 6, 2019 Published: June 28, 2019 A

DOI: 10.1021/acs.jpcb.9b02840 J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B alkanes,14−16 while many studies focused on their binary and ternary mixtures.17−23 For example, Queimada et al.24 performed measurements of liquid density, viscosity, and surface tension for several mixtures of heavy (n-eicosane, ndocosane, n-tetracosane) with smaller n-alkanes (n-heptane, ndecane, or n-hexadecane) from 293.15 up to 343.15 K at atmospheric pressure. Recently, Koller et al.25 determined the viscosity and surface tension of pure n-dodecane, n-octacosane, their binary mixture, and the SX-70 wax for temperatures between 323 and 573 K from surface light scattering experiments close to saturation conditions. On the other hand, molecular simulation has been established in recent years as a suitable tool to obtain such properties for n-alkanes and has greatly advanced in accurately reproducing the structural, dynamic, and thermodynamic properties that are macroscopically measured in actual physical systems, often acting complementary to experiments. Intermolecular interaction potential models or force fields lie at the core of molecular simulation. The availability of experimental data is not only a source for their parametrization but also provides the necessary means to evaluate and validate their performance. In this respect, numerous molecular dynamics (MD) and Monte Carlo (MC) simulation studies have been reported on liquid−vapor interfacial properties and surface tension of n-alkanes using a great variety of force fields.26−30 Harris31 was one of the first to present MD studies examining n-decane and n-eicosane using the united atom (UA) optimized potentials for liquid simulations (OPLS) force field (OPLS-UA)32 to calculate the liquid−vapor interfacial structure and surface tension. Nicolas and Smit33 computed liquid densities and surface tensions of n-hexane, n-decane, and n-hexadecane at different temperatures, comparing the OPLSUA32 and the Smit−Karaborni−Siepmann (SKS)34 force fields. They found that liquid densities were better described by SKS at high temperatures while surface tension was underestimated by about 15% at low temperatures.33 Ibergay et al.35 performed Gibbs ensemble MC simulations with the transferable potential for the phase equilibria (TraPPE-UA)36 force field in the canonical ensemble (NVT) to study methane, n-pentane, and n-decane interfacial properties of the VL interface over the temperature range 120−510 K. Their results highlight the chemical equilibrium of the liquid−vapor interface by calculating a local chemical potential including the appropriate long-range correction profiles by applying an extended approach of the “test-area” method developed by Gloor et al.37 Ismail et al.38 studied the surface tension of normal and branched alkanes using the fully atomistic (AA) OPLS variant (OPLS-AA)39 and the exponential-6 force field of Smith et al.,40,41 suggesting that the handling of long-range interactions is an important factor in determining surface tensions for alkanes. Mendoza et al.42 performed MD simulations in the NVT ensemble over the range of 200− 525 K for ethane to n-hexadecane at the VL interface using the NERD43 and TraPPE-UA36 force fields, revealing the good agreement of the TraPPE-UA model with experimental data, while NERD slightly overestimated γVL.42 Müller and Mejiá 10 performed MD simulations for n-decane, n-eicosane, nhexacontane, and n-decacontane using the SKS,34 NERD,43 and TraPPE-UA36 force fields. Their results showed that the TraPPE-UA calculated phase equilibria properties of n-decane and n-eicosane exhibited the lowest deviations compared to the experiment, but for n-hexacontane and n-decacontane the reliability of the all UA potentials decreases, underpredicting

saturated liquid density and vapor pressure, with TraPPE-UA potential exhibiting the lowest liquid density deviations.10 Chilukoti et al.44 used the NERD model to calculate density profiles, interfacial structure properties, self-diffusion coefficients, as well as surface tension values for n-butane, n-hexane, and n-decane from 295 to 450 K, showingin agreement to Mendoza et al.42that this model follows the temperature dependence of γ but yields slightly overestimated values. Galicia-Andrés and Medeiros45 studied n-hexane, n-heptane, noctane, n-nonane, and n-decane in the range of 300−420 K via the square gradient theory and MD simulations with the TraPPE-UA36 force field, showing that surface tension values were reliably predicted using a very large 2.4 nm cut-off without including tail corrections. Leonard et al.46 performed MD simulations of pure n-heptane, n-pentadecane, and nhexadecane with the CHARMM3647 additive and CHARMM Drude polarizable48 force fields at various biologically relevant temperatures, using the Lennard-Jones (LJ) particle-mesh Ewald (LJ-PME) method49−52 to calculate long-range LJ interactions. Their results showed that the implementation of LJ-PME improves agreement with the experiment for predictions made with these force fields of properties such as density, surface tension, and viscosity.46 Recently, Morrow and Harrison53 performed MD VL equilibrium simulations of several compounds typically found in diesel and jet fuel, among which n-dodecane, n-hexadecane, n-octadecane using the TraPPE-UA,36 OPLS-AA,39,54 and CHARMM55 force fields with the LJ-PME method. Their results showed that compared to the experiment, TraPPE-UA calculated critical temperatures (Tc) are the most accurate and surface tension values are overpredicted at all temperatures, while OPLS-AA produces surface tension lower than that of the experiment, which deviates by a larger extent as the temperature increases. Viscosity of n-alkanes is typically obtained using nonequilibrium MD (NEMD)56−66 or equilibrium MD (EMD) using the Einstein or Green−Kubo (GK) relations,54,67−70 covering a broad range of linear and branched n-alkanes with chain lengths from n-C10 up to n-C40. For more on the individual methods the reader is directed to the relevant literature71−73 as well as to the work of Hess.74 Briefly, NEMD simulations offer a computationally efficient means of calculating viscosity by employing arduous physics. Most of the previous efforts to study the viscosity of nalkanes make use of various force fields for the description of the system, which ultimately determine the accuracy of predictions. For example, Mundy et al.75 studied n-decane by performing both EMD-GK and NEMD simulations using a flexible variant of the TraPPE-UA force field at T = 480 K and ρ = 0.6136 g cm−3, pointing out the need for long simulation times in order to obtain good statistics with the GK method, as it converges slowly for long n-alkanes. Allen and Rowley76 performed NEMD NVT viscosity simulations of branched and linear n-alkanes (n-butane, n-decane, and n-dodecane) using the united-atom Ryckaert and Bellemans model,77 OPLS-UA32 and OPLS-AA39 force fields in the range of 150−523 K. Their results shown for n-alkanes in particular, that UA models produce similar results and predict viscosity values below the experimental measurementswith the exception of nbutaneespecially at higher densities. They also noted that UA and AA OPLS models produce quite different viscosity results despite yielding equivalent thermodynamic properties, which were attributed to the lack of detailed hydrogen atom representation, as hydrogen atom intermolecular potentials B


Article

The Journal of Physical Chemistry B have a strong effect on viscosity.76 Dysthe et al.78 performed EMD simulations for n-butane, n-decane, n-hexadecane, and 2methylbutane at different state points for seven different UA models using the GK formalism. Their results showed that the viscosities were increasingly underestimated with increasing packing fraction (i.e., the ratio of the density of the fluid to the “close packed” density of the fluid) and increasing chain length, highlighting the importance of a force field’s torsional potential coupling to viscosity.78 Gordon investigated a modified n-6 LJ potential in its ability to reproduce the transport properties of n-alkanes over a range of temperatures, pressures, and chain lengths, showing that below n-C10, TraPPE-UA calculated liquid densities are slightly overestimated, while GK calculated viscosity values in the NVT ensemble were systematically below experimental values.79 Ungerer et al.9 calculated viscosity from EMD by employing the Einstein relation using the TraPPE-UA36 and the anisotropic UA80 models, concluding that these force fields underestimate viscosity and modifying the torsion potential improves the results without compromising phase equilibrium properties.9 Payal et al.81 used EMD simulations to calculate the shear viscosity of n-decane and n-hexadecane at ambient as well as at high temperature−high pressure conditions using TraPPE-UA36 and the AA model of Tobias, Tu, and Klein (TTK-AA). 82 Their results showed that both models reproduce experimental densities, while viscosity values obtained from the TraPPE-UA and the TTK-AA force fields were found to be approximately 20% lower and 10% higher compared to their experimental counterparts, respectively. The authors concluded that the AA representation is crucial under extreme conditions in which hydrogen−hydrogen interactions are expected to be important.81 Makrodimitri et al.83 performed long MD simulations with the TraPPE-UA36 force field using the GK formulation to study the viscosity of pure nalkanes from n-C8 up to n-C96 as well as of the binary and sixcomponent model n-alkane mixtures relevant to the gas-toliquids process.83 Recently, Ewen et al.84 performed MD simulations on n-hexadecane, evaluating the performance of several AA and UA force fields. Their results showed that density prediction was generally accurate compared to the experiment for all force fields considered, while L-OPLS54 reproduced the most accurate viscosity values (within 15%) at the conditions simulated.84 From this literature overview it becomes apparent that the accuracy of the molecular models describing the intra and intermolecular interactions that are employed in simulations are of outmost importance for attaining accurate liquid density, surface tension, and viscosity values. To this extent, the main goal of this work was to benchmark the performance of wellestablished atomistic and coarse-grained (CG) molecular models developed primarily for biomolecular simulations in systems of interest to the chemical/petrochemical industry such as n-alkanes and their mixtures. Specifically, the LOPLS,54 Lipid14,85 and MARTINI86 models were benchmarked against the data reported by Koller et al.25 and the TraPPE-UA36 force field on the reproduction of the aforementioned properties for n-dodecane (n-C12), n-octacosane (n-C28), their binary mixture at a n-C12 mole fraction of 0.3 (referred to simply as n-C12/n-C28), as well as a representative model mixture of the commercially available SX-70 wax, for temperatures in the range of 323.15−573.15 K, mostly at ambient pressure.

2. SIMULATION METHODS 2.1. Force Fields. In this work, four different force fields covering AA, UA, and CG approaches were employed, namely, the fully atomistic L-OPLS and Lipid14, TraPPE-UA, and the MARTINI models, respectively. The L-OPLS force field is a reparametrization of the original OPLS-AA in order to accurately reproduce liquid properties of long hydrocarbons.54,87 Lipid1485 is the AMBER force field for lipids and features optimized alkyl chain dihedral angle terms and LJ parameters, allowing for tensionless simulations of such systems and offering enhanced compatibility and transferability to a wide range of physical systems, ranging from proteins and lipids to organic compounds.85,88,89 Partial charges for both molecules were generated according to the original Lipid14 derivation procedure,85 the details of which can be found elsewhere.90 The TraPPE-UA force field36 was employed for the UA representation of the n-alkanes examined, as it has been shown to be accurate for the reproduction of liquid densities in the pure state and in mixtures over a wide range of temperatures and pressures. Coarse Grained Molecular Dynamics (CGMD) simulations were performed using the MARTINI force field (version 2.0).86 The MARTINI force field was also originally developed for lipid simulations91 and has proved quite popular in simulations of not only lipids,92 but also biopolymers,93 proteins,94 and nanoparticles.95 It is based on a four-to-one mapping, that is, four heavy atoms are represented by a bead or single interaction center. n-Alkanes were parameterized using the standard C1 bead mapping,86 with realistic masses for each bead. The default combination rules for the size parameter σ and well depth ε of the nonbonded LJ interactions between sites of different types were used (i.e., arithmetic and geometric means for σ and ε parameters for the TraPPE-UA and Lipid14 force fields, respectively, geometric means for both σ and ε parameters for the L-OPLS and MARTINI force fields). Structures, partial charges, and force field parameters of all the n-alkane molecules examined are collected and presented in Tables S1−S6 of the Supporting Information document. 2.2. Simulation Details. All simulations were prepared and executed by means of the Scienomics MAPS 4.2 software package.96 Initial conformations of all systems were constructed using the amorphous builder module, according to which the constituent molecules are inserted into the cubic boxes with a modified configuration bias scheme.97 The number of molecules used was 280 for n-C12, 120 for n-C28, 44 n-C12, and 101 n-C28 for the n-C12/n-C28 binary mixture, respectively. For the model representing the SX-70 wax, the nalkanes used and their respective molecular numbers are collected in Table 1. All MD simulations for the systems in question were performed with the GROMACS 2016.5 software98−103 using the double precision compilation and imposing periodic boundary conditions in all directions. For the TraPPE-UA and Lipid14 simulations, the cut-off distance for LJ interactions was set to 1.4 nm, while the LOPLS LJ interactions were treated using a switch function104 between 1.1−1.3 nm, in accordance to the original work of Siu et al.54 Beyond the cut-off, long-range dispersion corrections to both the internal energy and to the pressure were considered, except in the case of the surface tension simulations. The PME method105,106 was applied for the calculation of long-range C


Article

The Journal of Physical Chemistry B

liquid density converged to a mean value, along with the energy of the system. Before performing the NVT surface tension and viscosity production runs, all systems were subjected to simulations for the evaluation of their liquid density over a period of 4 ns in the NPT ensemble. As far as n-C12 is concerned, the pressure was set to 1 atm in the range 323.15−473.15 K, while it was gradually increased in the range 498.15−573.15 K in order to maintain the liquid phase. A constant pressure of 1 atm was also used for pure n-C28 and all mixtures examined. The reported properties were calculated using coordinate sets representing the average liquid density value predicted by each force field at each state (Figure 1). Calculated liquid densities along with structural data and n-alkane chain relaxation times were obtained from a total of 20 000 frames.

Table 1. Composition of the Model SX-70 Wax n-alkane

number of molecules

n-C20 n-C24 n-C28 n-C32 n-C36 n-C40 n-C44 n-C48 n-C52 n-C56 n-C60 n-C80

1 19 128 260 265 188 93 32 9 2 1 1

electrostatic interactions at the same cut-off distance for LOPLS and Lipid14 simulations. Neighbor searching was performed using the “Verlet” scheme, with a list created every 5 steps using a length of 1.5 nm. The LINCS algorithm107 was used to constraint all bonds in TraPPE-UA and the bonds involving hydrogen atoms in the L-OPLS and Lipid14 simulations, respectively. The integration time step used in the TraPPE-UA simulations was 2 fs, while for the LOPLS and Lipid14 simulations it was set to 1 fs. For the simulations using the MARTINI force field, a time step of 20 fs was used, and a neighbor searching was performed using the “group-scheme” with a list length of 1.4 nm created every 5 steps.108 LJ interactions and forces were smoothly shifted to zero using the GROMACS shifting function between 0.9 and 1.2 nm. Simulation parameters are summarized in Table 2. Temperature and pressure were maintained constant using a Nosé−Hoover thermostat109 and Parrinello−Rahman barostat,110 respectively, with coupling constants equal to 0.2 and 2 ps−1 for the L-OPLS, Lipid14, and TraPPE-UA simulations, while for MARTINI these were set equal to 2 and 12 ps−1, respectively. Pressure coupling was isotropic in all directions. Prior to production runs, all systems were subjected to the steepest descent energy minimization for 20 000 steps for the elimination of any close contacts between atoms or beads. Short equilibration simulations of 100 ps and 1 ns in the NVT and isobaric−isothermal (NPT) ensembles, respectively, were performed in order for pure n-C12, n-C28, and the n-C12/n-C28 binary mixture systems to gradually reach their target temperatures and equilibrate their density. The model SX-70 mixture was equilibrated for 40 ns. During this period, the

Figure 1. Equilibrated conformations of n-C12 representing the average liquid density of each force field at 323.15 K. The figure was prepared by means of Scienomics MAPS 4.2.96

2.3. Surface Tension. The VL surface tension, γVL, NVT simulations were performed in orthorhombic boxes of dimensions Lx = Ly = 1/2Lz, with periodic boundary conditions imposed in all directions (Figure S1a). The simulation cells had an equilibrated size sufficiently large enough to avoid finite-size effects and allow the use of a 2.3 nm cut-off. For example, in the case of n-C28 (the longest chain tested), for the entire range of temperatures examined, the box size ranged from 4.7 to 5.0 nm, which is larger than the rootmean-square end-to-end distance (2.05 ± 0.05 nm at the 398.15 K) of the chains. This is also the case for the rest of the systems examined. For n-C12 and temperatures above 473.15 K, systems of greater size were employed by quadrupling the number of n-C12 molecules (i.e., 1120) and tripling the simulation cell along the z axis, in order to ensure the formation of a stable interface (i.e., Lz = 3Lx). VL surface tension γVL was then evaluated from a single 40 ns run per state

Table 2. Summary of the Simulation Settings for the Force Fields Used parameter time step (fs) RCoulomb (nm) Coulomb method RLJ (nm) nonbonded cut-off scheme neighbor list length (nm) neighbor list search (steps) constraints dispersion correction thermostat τT (ps) barostat τp (ps)

TraPPE-UA

MARTINI

2

20

1.4 Verlet 1.5 5 all-bonds EnerPres Nosé−Hoover 0.2 Parrinello−Rahman 2

0.9−1.2 group 1.4 5

Nosé−Hoover 2 Parrinello−Rahman 12 D

L-OPLS

Lipid14

1 1.3 PME 1.1−1.3 Verlet 1.5 5 h-bonds EnerPres Nosé−Hoover 0.2 Parrinello−Rahman 2

1 1.4 PME 1.4 Verlet 1.5 5 h-bonds EnerPres Nosé−Hoover 0.2 Parrinello−Rahman 2 DOI: 10.1021/acs.jpcb.9b02840 J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B point by using the diagonal stress tensor elements (Pxx, Pyy, and Pzz) through the following expression31 γVL =

Lz 2

jijP − Pxx + Pyy yzz jj zz zz 2 k {

respectively. It should be noted here that the maximum relaxation time calculated from the end-to-end vector autocorrelation function is smaller than 0.5 ns. The integration time was based on the same criterion suggested by Zhang et al.,118 that is, the time up to which the standard deviation of the running integrals σ is less than 40% of the corresponding viscosity value. On the other hand, nonequilibrium methods calculate shear viscosities either by imposing a linear velocity profile along the shear direction in the simulation box (p-SLLOD119,120 and SLLOD121,122 algorithms), applying a momentum flux to obtain the velocity gradient123 or by applying a force profile to the simulation box and then measuring the viscosities from the fluctuations of induced momentums (periodic perturbation method, PP).124,125 Hess’s review74 demonstrated that the PP method is efficiently applicable and accurate. In this respect, Zhao et al.126,127 have successfully applied this method in the prediction of several molecular liquids, including n-alkanes, with promising results. In the PP method, an external force is applied to the system, the magnitude of which can be chosen such that the effects are much easier to measure than the internal fluctuations. Important factors affecting the method’s performance, besides the quality of the force field in use, are the wavelength and the magnitudes of the applied periodic forces, as presented in detail in the review by Hess.74 Briefly, a periodic external force a(z), a function of the z-direction of the 3D periodic cell is applied in the x-direction on each of the particles

(1)

where, Lz is the length of the simulation box in the z direction, normal to the VL interface which forms along the xy plane. The factor 1/2 is due to the presence of two interfaces and angle brackets represent the ensemble average. The values were estimated from 20 000 frames per state point trajectory, while statistical uncertainties (95% confidence intervals) were obtained by dividing the production period of each simulation into five blocks.111 The cut-off of the LJ interactions influences the computed values of surface tension significantly, particularly in simulations of inhomogeneous systems such as VL interfaces, often leading to the incorrect estimation of such properties.112,113 In order to compensate for the truncation of dispersion forces, large values of the cut-off are typically used in combination with analytical corrections114,115 or lattice-based methods such as the LJ-PME49−52 to calculate long-range LJ interactions.116 In this study, the long-range LJ potential was cut-off using a radius of 2.3 nm, and the effect of using the LJ-PME49−52 was also investigated. The LJ-PME technique is computationally robust and has been previously employed successfully in the calculation of γVL.46,53,115 The goal was to investigate the effect of LJ-PME on the VL surface tension of n-alkanes, because previous works have shown that using this approach with force fields optimized for a cut-off approach leads to overestimated γVL values for several organic liquids.117 All the aforementioned LJ interaction treatments concern all models except MARTINI, where only the default simulation parameters described previously were employed. 2.4. Viscosity. As described in the Introduction section, the methods available for calculating viscosities of liquids from MD simulations are subdivided to (i) equilibrium and (ii) nonequilibrium. The former makes use of pressure fluctuations, either using Einstein’s relation or the GK formula72,74 ηGK =

V kBT

∫0

ax(z) = A cos(kz)

where, the amplitude A is the acceleration amplitude, which specifies the magnitude of the force applied and the wavenumber k k=

2π Lz

(4)

with a wavelength equal to the box length in the z direction, Lz. According to the Navier−Stokes equation,128 when the system reaches the steady state the generated oscillatory velocity profile is

∞

⟨Pαβ(t0)Pαβ(t + t0)⟩dt

(3)

(2)

ux(z) = u0 cos(kz)

where, V is the simulation box volume, T is the temperature, kB is the Boltzmann’s constant, and Pαβ denotes the off-diagonal elements of the pressure tensor. Angle brackets indicate an ensemble average over all time origins. Calculating the viscosity using the GK formulation is not trivial. Even for molecules of low molecular mass, a fair estimation of viscosity and its standard deviation requires a large number of independent simulations and long simulation times, rendering the calculation impractical if at all possible. Therefore in this study and for comparison reasons, only n-C12 viscosities were calculated using the TraPPE-UA force field and the GK method, following the strategy proposed by Zhang et al.118 by extracting shear viscosity estimates from 100 EMD simulations in the microcanonical ensemble (NVE), 2 ns long with a 0.5 fs time step where the pressure tensor was sampled every 5 fs. The NVE simulations were conducted in the temperature and pressure conditions of interest using configurations with simulation boxes set to the mean liquid density, ⟨ρ⟩, and total energy close to the mean value predicted from the TraPPE-UA force field for these conditions. Both ⟨ρ⟩ and ⟨E⟩ were obtained from NPT and NVT equilibration simulations,

(5)

where, u0 is the average velocity profile amplitude. Thus, viscosity is obtained using the following formula ηPP =

A ρ u0 k2

(6)

where, ρ is the equilibrium density. Viscosity depends on the acceleration amplitude, A, and optimal values should be small enough so that the perturbation does not cause shear rates and temperature drifts which disturb the equilibrium of the system significantly.74 In addition, utilizing a large system is optimal as it produces a longer wavelength.74,126 To this end, all systems were doubled in size along the z direction (Figure S1b). Values for A were selected by performing a series of 4 ns NVT simulations following the procedure proposed for LJ fluids by Vasquez et al.129 in the ranges of 0.005−0.20 nm ps−2 for n-C12 and n-C28, and of 0.0005−0.020 nm ps−2 for the model SX-70 wax, respectively. Test simulations were performed at the lowest and highest temperatures examined for each system, using conformations with the average liquid E


Article

The Journal of Physical Chemistry B density as calculated from the NPT equilibration simulation. For comparison, two force fields, namely, TraPPE-UA and LOPLS were examined, except for the model SX-70 where only TraPPE-UA was considered. Final acceleration amplitudes were selected from the plots of the computed viscosities from the regions, where ηPP is fairly constant with respect to A (Figures S2−S4),129 so that the linear Newtonian behavior is followed at the two extreme temperature end points for all systems. We have therefore used the following four acceleration amplitudes for n-C12, namely, A = 0.0075, 0.0100, 0.0125, and 0.0150 nm ps−2 and for n-C28: A = 0.005, 0.010, 0.015, and 0.020 nm ps−2. The latter acceleration amplitudes were also employed for the n-C12/n-C28 binary mixture. Lastly, for the model representing wax SX-70 the same procedure the acceleration amplitudes were: A = 0.0020, 0.0025, 0.0030, and 0.0035 nm ps−2. Production runs were performed for 4 ns in the NVT ensemble and analysis was performed over 20 000 frames collected from the last 2 ns of each NEMD-PP simulation. The obtained viscosities along with the uncertainties of each point were estimated using the block average method111 with respect to the acceleration amplitudes applied in the simulations, which are illustrated in Figures S5−S18. In all cases, the desired linear correlation exists in the region of the selected acceleration. At each temperature, a least-squares-fitting of the viscosity values was performed using the uncertainty of each point as the weight, and the shear viscosities for the “unperturbed” systems were retrieved by extrapolating the line to A = 0.126,127 The uncertainty for the extrapolated viscosities is equal to the 95% confidence interval of the root of mean square of deviations of the fits.

experimental (Expt.) for convenience. The average absolute deviation (AAD), that is, the absolute difference between reference and calculated values over the reference data, is only 1% for TraPPE-UA, following the correct trend with increasing temperature. On the other hand, the MARTINI force field follows the correct temperature dependence, but ultimately yields overestimated liquid density values with an AAD equal to 9%. This has also been reported previously131,132 and it is not surprising, as the MARTINI force field has been optimized to reproduce the free energies of vaporization, hydration, and partitioning between water and organic phases at ambient temperature and pressure conditions (300 K, 1 bar),86,91 which are commonly found in biological systems. Moreover, the fully atomistic L-OPLS and Lipid14 force fields are very successful in predicting liquid density of n-C12 at the lower end of the tested temperature range, yielding a value that is within 1 and 0.2%, respectively, from the experiment at 323.15 K. However, as the temperature increases to 573.15 K, both force fields result in significant deviations from the reference data, having a value that underestimates Expt. data by 24 and 23%, respectively. Interestingly, both AA force fields display the same temperature dependence in terms of liquid density prediction. Similar findings are observed in the case of pure n-C28 (Figure 2b and Table S8), where the TraPPE-UA force field produces liquid densities that closely follow Dutour’s experimental correlation values (0.1% uncertainty)133 reported by Koller et al.25 (also referred to as Expt.), being slightly overestimated by 2% on average and successfully reproduce the temperature dependence. The MARTINI force field again captures the correct liquid density trend as the temperature increases, however, the values are 5% too high as compared to the experiment. It is deduced that as the n-alkane chain length increases, the overestimation is not so pronounced and the results agree better with the Expt. values. Lipid14 shows excellent agreement within 0.4% from the reported liquid density of n-C28 at the lowest temperature examined (i.e., 398.15 K), while L-OPLS slightly underestimates it by 2.4%. As with n-C12, both force fields prove to capture liquid density less accurately with increasing temperatures, but deviations are less pronounced in this case, being approximately 9% lower than the experiment at 573.15 K. The same findings are applied in the case of the binary n-C12/n-C28 mixture, respectively (Figure 2c, Table S9). Overall, these results indicate that the TraPPE-UA model yields excellent agreement to the experimental liquid density predictions of the pure n-C12 and n-C28 (AAD of 1 and 2%, respectively), with a mean deviation from experimental data in close agreement with the previous simulation studies on nC1236 as well as on n-C10 and n-C16 (reported AAD ∼1%).87 MARTINI overestimates liquid density, particularly in the case of n-C12, but still produces the correct temperature dependence. L-OPLS and Lipid14 atomistic force fields show similar behavior, having excellent agreement at the lower temperature range, but as the temperature increases, greater deviations are observed, especially for n-C12. Therefore, for the prediction of liquid densities, TraPPE-UA is found to be the most accurate among the models tested; MARTINI offers the CG advantages in terms of computational efficiency, and in agreement with previous studies yields good predictions131 as well as appearing to be more effective with increasing chain length, while LOPLS and Lipid14 work well in the lower temperature range.

3. RESULTS AND DISCUSSION 3.1. Liquid Density, Structural, and Relaxation Properties. From the liquid density values of pure n-C12 illustrated in Figure 2a (Table S7), it is evident that the TraPPE-UA force field closely reproduces the Lemmon and Huber equation of state reference data (0.2% uncertainty)130 reported in the work of Koller et al.,25 which are referred to as

Figure 2. Liquid densities for (a) n-C12, (b) n-C28, and (c) the n-C12/ n-C28 binary mixture. The black line and squares corresponds to the reference values reported by Koller et al.25 while colored symbols correspond to the MD predictions made using the different force fields employed in this study: red diamonds correspond to TraPPEUA, green diamonds to MARTINI, blue triangles pointing up to LOPLS, and magenta triangles pointing down to Lipid14. Dotted lines with the respective colors serve as guides to the eye. F


Article


UA with a 2.3 nm LJ cut-off yields γVL results that are in excellent agreement with the experimental measurements, with an AAD equal to 4%. This is in line with the results from previous TraPPE-UA simulations with a 2.4 nm cut-off for nC10 in the temperature range 350−600 K, where the AAD from experiment was around 6%.10 In addition, VLE liquid density values are in close agreement with their NPT counterparts (Figure S21) and the predicted Tc ((654 ± 1) K) is also in very good agreement with the Expt. data (659 K) (Table S11). MARTINI is close in agreement, as it provides the second best γVL values with AAD of 6% from the experiment (Table S12), while predicting a Tc for n-C12 ((660 ± 2) K) slightly above the Expt. data (Table S11). It should be noted that MARTINI yields these values without any special treatment of the longrange LJ interactions, but only using the force field defaults. Both force fields display their largest deviation from the experiment at the high end of the temperatures examined (i.e., 548.15 and 573.15 K). On the other hand, both AA models produce poor γVL values. Specifically, L-OPLS and Lipid14 underestimate surface tension by 34 and 39% on average, respectively. The agreement with the experiment is slightly better for lower temperatures, but significantly worsens as the temperature increases, with the trend resembling the decay in liquid density predictions presented in the previous section. Both these force fields produce Tc values well below the Expt (Table S11). It should be stressed that as the temperature approaches the critical temperature predicted by these force fields (T > 0.9Tc), recommendations concerning the slab simulation box dimensions and the cut-off radius are no longer valid and one needs to proceed with caution.116,134−136 The VL interface was stable in all cases except at 573.15 K for LOPLS and Lipid14, where no interface was formed and surface tension could not be calculated. The slab method used for calculating the surface tension in this work has severe limitations at the critical temperature region for n-C12 and alternative methodologies that avoid an explicit interface should be applied.137,138 γVL values that fall to this region are selectively depicted in Figure 3 (empty symbols) and Table S12 (red color). The LJ-PME method was also utilized for the prediction of γVL for n-C12 using TraPPE-UA, L-OPLS, and Lipid14 (Figure 3b, Table S12). This method is computationally robust but has not been extensively tested in conjunction with force fields optimized for the cut-off approach.117 As far as TraPPE-UA is concerned, using LJ-PME leads to overestimated values by 10% on average, falling in close agreement with the results reported for n-C12 by Morrow and Harrison using the same approach in the temperature range of 461−560 K.53 The Tc calculated is slightly higher in value than when a 2.3 nm cutoff is used (Table S11). On the other hand, using LJ-PME with the L-OPLS and Lipid14 force fields still produces underestimated values compared to experiment, but improves agreement with experiment as they now deviate by 30 and 33% on average, respectively (Figure 3b, Table S12). L-OPLS results with the LJ-PME methods are also in line with Morrow and Harrison.53 Agreement is excellent at 323.15 K (1 and 3%, respectively), as with the long cut-off approach, but worsens as temperature increases. As far as n-C28 is concerned, the combination of TraPPE-UA with a 2.3 nm LJ cut-off again yields results that are in excellent agreement with the experiment, slightly overestimated by 2% on average from the experiment (Figure 4a, Table S13). This compares very well the AAD from the experiment for n-C20

The structural properties of the pure component systems were also analyzed. The predicted content of the transdihedrals by TraPPE-UA, L-OPLS, and Lipid14 is reduced with temperature. This is the case among all these force fields and systems examined, hence only TraPPE-UA calculated torsion probability distributions for n-C28 for the temperatures considered are depicted in Figure S19. MARTINI CG description inherently prohibits such conformational temperature dependence because there are no dihedral contributions to the potential energy for this force field. Therefore, for the TraPPE-UA, L-OPLS, and Lipid14 force fields, the n-alkane molecules examined adopt a less extended configuration with temperature increase. This is demonstrated by the decrease of the mean squared end-to-end vector, ⟨Ree2⟩, and mean squared radius of gyration, ⟨Rg2⟩, with temperature (Figure S20). As far as MARTINI is concerned, calculated ⟨Ree2⟩ and ⟨Rg2⟩ for n-C28 (n-C12 comprises only three beads and calculating these quantities was deemed unnecessary) were found to be temperature independent due to the inherent CG description (Figure S20c), rendering their comparison invalid. In terms of relaxation, the maximum orientation relaxation time, t0, (i.e., a dynamic property describing the time needed for the molecule orientation to decorrelate) followsas expectedthe same trend (Table S10) in all cases. 3.2. Surface Tension. VL slab simulations allowed the evaluation of surface tension, γVL, the estimation of VLE coexisting curves (Figure S21) and critical temperatures (Tc) predicted by each force field (Table S11). For n-C12 (Figure 3a, Table S12), it is deduced that the combination of TraPPE-

Figure 3. Experimental data25 and MD predictions from each force field for the surface tension of n-C12, utilizing (a) a long cut-off radius (2.3 nm) for the treatment of LJ potential in all cases except MARTINI and (b) the LJ-PME approach with TraPPE-UA, L-OPLS, and Lipid14 force fields. Symbols are the same as in Figure 1. Empty symbols illustrate points for L-OPLS and Lipid14 where T > 0.9Tc.116 Surface tension could not be determined with L-OPLS and Lipid14 at 573.15 K. G


Article


Figure 4. Experimental data25 and MD predictions from each force field for the surface tension of n-C28, utilizing (a) a long cut-off radius (2.3 nm) for the treatment of LJ potential in all cases except MARTINI and (b) the LJ-PME approach with TraPPE-UA, L-OPLS, and Lipid14 force fields.

Figure 5. Experimental data25 and MD predictions from each force field for the surface tension of the n-C12/n-C28 binary mixture, utilizing (a) a long cut-off radius (2.3 nm) for the treatment of LJ potential in all cases except MARTINI and (b) the LJ-PME approach with TraPPE-UA, L-OPLS, and Lipid14 force fields.

(∼0.5%) from Müller and Mejia.́ 10 As in the case of n-C12, MARTINI γVL values are only 4% on average from the experiment. L-OPLS and Lipid14 underestimate surface tension by 18 and 19% on average, respectively, performing better at lower temperatures showing significant deviations with increasing temperature. However, deviations are less pronounced compared to n-C12. The use of the LJ-PME method with TraPPE-UA yields overestimated γVL values by 15% on average, while LJ-PME with L-OPLS and Lipid14 improves agreement with the experiment, reducing AAD by 13 and 11%, respectively, showing excellent agreement at 398.15 K (AAD equal to 2 and 0.2%, respectively), which diminishes as the temperature increases (Figure 4b, Table S13). Similar trends are observed for the n-C12/n-C28 binary mixture, where again TraPPE-UA using a 2.3 nm LJ cut-off results in excellent agreement with AAD from the experiment, which is only 2% (Figure 5a, Table S14). MARTINI values are on average moderately lower by 4% from the experiment. LOPLS and Lipid14 underestimate surface tension with an AAD of 18 and 17%, respectively. LJ-PME TraPPE-UA simulations yield overestimated γVL values by 13% on average, while LJPME with L-OPLS and Lipid14 improves agreement with the experiment, reducing AAD to 13 and 10%, respectively (Figure 5b, Table S14). From these results, it is deduced that the implementation of TraPPE-UA with the use of large cut-offs produces γVL values with the best experimental agreement for the n-alkanes examined and their binary mixture. This is in line with previous studies where it has been reported that the cut-off radius value has a profound effect on the simulations for the TraPPE-UA potential,10 along with the slight tendency to overestimate surface tension by a few percent for different

temperatures and chain lengths.27 The MARTINI force field appears to be a very attractive alternative for VL surface tension predictions of n-alkanes because it offers the desired computational efficiency needed for systems of larger sizes, without the need for any special LJ treatment. It has been previously shown that MARTINI was very successful in the prediction of surface tension with respect to experiments for liquid−liquid interfaces of linear n-alkane−water systems, with deviations of about 20%.139 The present study of γVL for pure n-C12, n-C28, and their binary mixture shows excellent agreement with experiment, which affirms that the source of these deviations in linear n-alkane−water systems is the MARTINI water models.86,139,140 Lastly, both L-OPLS and Lipid14 produce poor γVL predictions, especially with increasing temperature, regardless of the long LJ interaction treatment, and their use should be restricted for the prediction of such interfacial properties at lower temperatures. 3.3. Viscosity. As shown in Figure 6a (Table S15) for nC12, the calculated viscosity values using the PP method at different temperatures follow the experimental trend in all cases. TraPPE-UA underestimates experimental viscosities with an AAD of 23% compared to the experiment, with its maximum deviation being 38% at 323.15 K. This underestimation is reduced to 8% at 573.15 K. The calculations for the viscosity of n-C12 following the GK formalism are presented in Figure S22 and Table S16 for the TraPPE-UA force field. The calculation protocol adopted118 resulted in converging viscosity values by increasing the number of independent simulations used for its estimation (Figure S23). The predicted viscosity underestimates the experimental values up to approx. 423.15 K (AAD is approx. 19%) and at higher temperatures overestimates them (AAD approx. 37%), and is H


Article


notable exception appears to be the predicted viscosity at 398.15 K, which is lower by 42% compared to the experiment, whereas in TraPPE-UA simulations it is underestimated by 30%. Both L-OPLS and Lipid14 force fields examined produce accurate viscosities for n-C28 with AADs equal to 6 and 0.2%, respectively. L-OPLS and Lipid14 yield overestimated viscosities at 398.15 K by 4 and 36%, respectively, while they perform best within the temperature range of 423.15− 498.15 K, with AADs of 1 and 6% from the experiment, respectively. As with the case of n-C12, in the temperature regime of 523.15−573.15 K, both models underestimate viscosity with AADs equal to 17 and 19%, respectively, but follow experimental results more closely. As shown from the results presented in Figure 6c (Table S18) for the n-C12/n-C28 binary mixture, the calculated viscosities reproduce the correct temperature dependence according to the experimental measurements. TraPPE-UA and MARTINI show the same behavior, the underestimated viscosities with AADs are equal to approximately 28 and 19%, respectively. TraPPE-UA is yields underestimated values by 35 and 30% at the lowest and highest temperatures examined, respectively. MARTINI performs slightly worse at lower temperatures, underestimating the viscosity by 41% at 373.15 K, but improves significantly as the temperature is increased beyond 473.15 K, yielding values that are lower by approximately 8% on average in the 523.15−573.15 K range. On the contrary, L-OPLS and Lipid14 show excellent agreement with experimental values, having AADs of approximately 9 and 4%, respectively, for the entire temperature range. Deviations for these force fields are more pronounced for lower temperatures (AADs equal to 10 and 27%, respectively, at 373.15 K), decrease significantly in the intermediate range (7 and 2%, respectively, at 423.15−498.15 K), and are increased again in the high temperature regime (19 and 21%, respectively, at 523.15−578.15 K). From the NEMD-PP viscosities calculated, one can draw several conclusions in terms of the performance of the force fields used, given that the simulation setups and acceleration amplitudes employed are identical for each system. Overall, TraPPE-UA and MARTINI are less accurate than the fully atomistic force fields, as L-OPLS and Lipid14 results are very close to the experimental values. The inadequacy of the TraPPE-UA model in the prediction of n-alkane dynamic properties such as self-diffusion coefficients and viscosities has been reported previously.70,79,81,83,87,141 Ungerer et al.9 as well as Dysthe et al.78 have pointed out the importance of a force field’s torsional potential coupling to viscosity, which can even be fine-tuned in order to improve results without compromising phase equilibrium properties.9 Even though this could be an important aspect for future attempts to reparametrize the TraPPE-UA force field for viscosity predictions, it is rather puzzling that the CG MARTINI force field, which completely lacks a torsional term, exhibits an analogous behavior. This similarity in viscosity prediction is attributed to the smoother potential energy surface that results from the neglect of atoms and diminished degrees of freedom in both force fields. Moreover, as far as the MARTINI force field is concerned, viscosity predictionseven though they are lower than experimental dataare within 17% on average from experimental measurements in all cases. Therefore, smoothening of the n-alkanes energy landscape by coarse graining with

Figure 6. Experimental and calculated shear viscosities using the NEMD-PP method for (a) n-C12, (b) n-C28, and (c) n-C12/n-C28.

in fair agreement with the NEMD-PP results. Moreover, the time decomposition method as proposed by Zhang et al.118 and applied in this study results in the correct viscosity− temperature dependence. The MARTINI force field exhibits a similar behavior using the PP method; however, the values are closer to their experimental counterparts as they are lowered by 12% on average, with the largest deviation found equal to 25% at 323.15 K and only 4.8% at 573.15 K. On the contrary, L-OPLS and Lipid14 exhibit a different behavior. Both force fields overestimate the viscosity at the lowest temperature, by 22 and 32%, respectively, while excellent agreement is observed for intermediate temperatures 348.15−498.15 K, where AAD is only 4 and 11% from the experiment. At temperatures in the range 523.15−573.15 K, both force fields show their highest deviations, underestimating viscosity with an AAD equal to 31 and 36%, respectively. The same trends are also observed for n-C28 (Figure 6b, Table S17). The zero shear extrapolated viscosity values using the PP method at different temperatures are overall following the experimental results, with TraPPE-UA and MARTINI force fields underestimating the experimental viscosities by 25 and 20% on average, respectively. Specifically, TraPPE-UA yields viscosity values that are underestimated by almost 30% at the lower temperature regime, but performs more accurately with increasing temperature with deviation of about 20% compared to the experiment. Our results show that as the chain length increases from n-C12 to n-C28, TraPPE-UA performs worse in reproducing viscosity using the PP method. The same is also true for MARTINI, which however seems to predict viscosity values that are in slightly better agreement with the experiment compared to TraPPE-UA for the entire temperature range, but particularly at higher temperatures. The I


Article

J

−3 (2) −17 (9) −10 (11) −5 (11) −18 (9) −20 (12) −12 (10)

The standard deviation is provided in parenthesis. The positive and negative percentage values indicate predictions above and below the reference data, respectively. The lowest values for each property and each of the three systems are indicated with bold numbers.

a

n-C28

−3 (2) −19 (9) −11 (10) 0.2 (20) −8 (7) −39 (25) −33 (26) −14 (22)

n-C12 n-C12/n-C28

−5 (2) −19 (12) −13 (9) −10 (6) −5 (2) −18 (10) −13 (10) −6 (10)

n-C28 n-C12

−8 (7) −34 (24) −30 (27) −10 (18) 6 (0.3) −5 (4)

n-C12/n-C28 n-C28

5 (0.7) −4 (4)

n-C12

9 (0.3) −6 (5)

2 (0.3) 2 (4) 13 (4) −27 (5)

n-C12/n-C28 n-C28 n-C12

1 (0.2) −4 (4) 10 (2) −23 (10)

property

ρ γVL γVL (LJ-PME) η

2 (0.3) 2 (4) 14 (3) −25 (6)

Lipid14 L-OPLS % deviation from reference data MARTINI TraPPE-UA

Table 3. Average Percentage Deviation between Molecular Simulation Values and the Reference Data25 for Liquid Density, VL Surface Tension, and Viscosity of the Pure nC12, n-C28, and Their n-C12/n-C28 Mixturea

MARTINI does not appear to result in effective time scales related to significant kinetics speed-up. This is in accordance to previous studies, which have shown that the CG speed-up of MARTINI is molecule dependent, and n-alkanes have a small speed-up compared to experiments or even a reduced speed compared to atomistic simulations.91,142 Our results are also in line with the works of Allen and Rowley,76 Payal et al.,81 and Ewen et al.,84 where it was shown that the AA representation is crucial in order to accurately capture viscosity because the detailed hydrogen atoms representation bears a strong effect on this property. Furthermore, even though L-OPLS and Lipid14 force fields have different torsional potential terms, they produce viscosities that overall follow the same trends in terms of temperature dependence and are close in value. Therefore, the present study is supportive of the importance of using an AA description for the description of viscosity. However, the exact role of the torsional potential requires further investigations. Unfortunately, using an AA representation increases the computational cost significantly, particularly as the system size increases. Moreover, the tested UA/CG and AA models produce liquid densities above and below the reference data, respectively, while the observed trends in viscosity are reversed, which is possibly related to the different conformational features or the free volume characteristics143 predicted for n-alkanes by each model, however, the latter assumption was not investigated any further. The use of the PP method appears to work very well for the prediction of n-alkane viscosities. Selecting an appropriate acceleration amplitude is very important and the computational cost of doing so is minimal. As it has been previously reported,126,127 the linear dependence of the PP calculated shear viscosity allows for linear extrapolation in order to predict the zero-shear viscosities. Utilizing the PP method in conjunction to the LOPLS force field overall produces the best agreement for the nalkane systems examined with experimental data; however, the L-OPLS model is computationally expensive which potentially hinders its implementation to larger systems. Lastly, the overall force field performance examined in terms of their percentage deviation from the reference data for each property considered, namely, liquid density, VL surface tension, and viscosity, is summarized in Table 3. The computational efficiency, which is defined as the nanoseconds per day (ns/day) ratio, obtained for a given force field over the MARTINI model (i.e., the one with the highest performance) is 1/320 for the AA force fields (L-OPLS and Lipid14) and 1/ 15 for the TraPPE-UA force field. 3.4. The SX-70 Wax Model. In light of the findings for the force fields examined in their performance towards the aforementioned properties, we modeled SX-70 wax (Figure S24a) using TraPPE-UA, L-OPLS, and MARTINI in the temperature range 473.15−523.15 K at 1 atm. The liquid density of SX-70 was experimentally obtained by volumetric measurements25 and as shown from the results presented in Figure 7a (Table S19), the calculated MARTINI values closely match the experimental line. It is also observed that the TraPPE-UA force field yields values that are lower from the experimental counterparts. Given the trends for the pure n-C12, n-C28, and their binary mixture, results for SX-70 is an indication of poor accuracy of the model for heavier n-alkanes (Figure S24b). On the other hand, this behavior is consistent with the study of Müller and Mejia,́ 10 where it was shown that

n-C12/n-C28



Article


fields developed for the study of biomolecular systems, namely, the L-OPLS, Lipid14, and MARTINI, on the reproduction of liquid densities, surface tension, and viscosities, against the well-established TraPPE-UA model as well as with the recently reported data by Koller et al.25 TraPPE-UA reproduces liquid densities most accurately and consistently over the entire temperature range. MARTINI overestimates liquid densities, particularly in the case of n-C12, but captures the correct temperature dependence. Both LOPLS and Lipid14 exhibit excellent agreement with experimental data for lower temperatures, but as the temperature increases, greater deviations are observed, especially for n-C12. For the special case of the SX-70 wax, TraPPE-UA predicts the lower liquid density values compared to the experiment, which is probably due to the tendency of UA models to produce lower liquid densities for longer n-alkanes.10 MARTINI over prediction of liquid density compared to TraPPE-UA results in excellent agreement with the experiment for the SX-70 system. As far as the VL surface tension is concerned, our simulations with TraPPE-UA using large cut-offs for the treatment of long LJ interactions yield the best agreement with the experiment, also in agreement with previous n-alkane studies.10 However, the MARTINI force field was found to perform equally well and is suggested as an alternative for such calculations. MARTINI offers enhanced computational efficiency and can be used for VL surface tension calculations, that is, without the need for any special LJ treatment. This was further corroborated in the case of the complex SX-70 n-alkane system, where MARTINI produced excellent results compared to the experiment. L-OPLS and Lipid14 produced poor predictions with increasing temperature, regardless of the LJ treatment. Our viscosity simulations using the PP method revealed that in agreement with previous studies,70,79,81,83,87,141 the TraPPEUA produces poor viscosity predictions in all n-alkane systems examined. MARTINI predictions are within 17% AAD from experimental measurements, therefore it is deduced that there is no effective time scale speed-up of this CG model for nalkane viscosity compared to experiments.91,142 MARTINI matches the behavior of TraPPE-UA, which is attributed to their coarser representation, a crucial factor in order to accurately capture viscosity.76,81,84 In any case, MARTINI is performing well for temperatures beyond 473.25 K, at a minimal computational cost, which was further demonstrated in the case of the SX-70 wax model. Lipid14 and L-OPLS, which produce very good results especially in the 423.15− 498.15 K range, outperform both TraPPE-UA and MARTINI. This is counterbalanced by their increased computational demands; therefore, their use is limited to relatively small system sizes. The GK method was tested on n-C12 with the TraPPE-UA model and was found to yield viscosity values similar to the NEMD-PP, suffering from the increased computational cost needed to overpass convergence problems. On the other hand, the PP method works very well for the prediction of n-alkane viscosities particularly with the L-OPLS force field, as it generally produces the best agreement with the experimental data. Overall, our simulations show that these widely used force fields are suitable candidates in the study of n-alkane properties, both in the pure and mixture states, by carefully considering the temperature range and properties of interest. There is however not a single outstanding force field which

Figure 7. (a) Experimental data and MD predictions for the liquid density of the SX-70 wax model system calculated with the TraPPEUA, L-OPLS, and MARTINI force fields. (b) Experimental data and MD predictions with TraPPE-UA and MARTINI force fields for surface tension. (c) Experimental data and MD predictions using the NEMD-PP method with L-OPLS and MARTINI for shear viscosity.

even though TraPPE-UA was able to reproduce liquid densities for n-C10 and n-C20, it predicted values below experimental measurements for n-C60 and n-C100. In any case, both force fields maintain the same trends observed previously, with MARTINI producing liquid densities that are above that of TraPPE-UA. L-OPLS consistently predicts liquid density values lower than that of the experiment, with an AAD equal to 9%. All models capture the SX-70 liquid density temperature dependence in the range considered. MARTINI yet again reproduces the γVL values in excellent agreement with the reported reference data, for all state points, being in AAD of 6% from experiments, while the TraPPE-UA calculated γVL values are slightly above Expt. with an AAD equal to 14% (Figure 7b, Table S20). NEMD-PP calculated viscosities with MARTINI and L-OPLS are also in excellent agreement with the experiment, with values that are in AAD equal to 2 and 6%, respectively (Figure 7c, Table S21). Therefore, it is deduced that MARTINI proves to be a promising alternative for surface tension and viscosity investigations of realistic complex n-alkane mixtures such as the SX-70 wax, coupling its good performance with the desired computational efficiency.

4. CONCLUSIONS The present work entails extensive MD simulations of n-C12, nC28, their binary mixture, and a model mixture of the commercially available hydrocarbon wax SX-70 over a broad temperature range from 323.15 to 573.15 K at ambient pressure. The goal was to evaluate the performance of force K


Article


(5) Israelachvili, J. N. Intermolecular and Surface Forces; 2nd ed.; Academic Press: London U.K., 1991. (6) Myers, D. Surfaces, Interfaces, and Colloids: Principles and Applications, 2nd ed.; Wiley-VCH; New York, U.S.A., 1999. (7) Diefendorf, A. F.; Freeman, K. H.; Wing, S. L.; Graham, H. V. Production of n-Alkyl Lipids in Living Plants and Implications for the Geologic Past. Geochim. Cosmochim. Acta 2011, 75, 7472−7485. (8) Tipple, B. J.; Berke, M. A.; Doman, C. E.; Khachaturyan, S.; Ehleringer, J. R. Leaf-Wax n-Alkanes Record the Plant-Water Environment at Leaf Flush. Proc. Natl. Acad. Sci. U.S.A. 2013, 110, 2659−2664. (9) Ungerer, P.; Nieto-Draghi, C.; Lachet, V.; Wender, A.; Di Lella, A.; Boutin, A.; Rousseau, B.; Fuchs, A. H. Molecular Simulation Applied to Fluid Properties in the Oil and Gas Industry. Mol. Simul. 2007, 33, 287−304. (10) Müller, E. A.; Mejía, A. Comparison of United-Atom Potentials for the Simulation of Vapor-Liquid Equilibria and Interfacial Properties of Long-Chain n-Alkanes up to n-C100. J. Phys. Chem. B 2011, 115, 12822−12834. (11) Jasper, J. J.; Kerr, E. R.; Gregorich, F. The Orthobaric Surface Tensions and Thermodynamic Properties of the Liquid Surfaces of the n-Alkanes, C5 to C28. J. Am. Chem. Soc. 1953, 75, 5252−5254. (12) Jasper, J. J.; Kring, E. V. The Isobaric Surface Tensions and Thermodynamic Properties of the Surfaces of a Series of n-Alkanes, C5 to C18, 1-Alkenes, C6 to C16, and of n-Decylcyclopentane, nDecylcyclohexane and n-Decylbenzene. J. Phys. Chem. 1955, 59, 1019−1021. (13) Doolittle, A. K.; Peterson, R. H. Preparation and Physical Properties of a Series of n-Alkanes. J. Am. Chem. Soc. 1951, 73, 2145− 2151. (14) Maeda, N.; Yaminsky, V. V. Surface Supercooling and Stability of n-Alkane Films. Phys. Rev. Lett. 2000, 84, 698−700. (15) Wu, Y.; Bamgbade, B.; Liu, K.; McHugh, M. A.; Baled, H.; Enick, R. M.; Burgess, W. A.; Tapriyal, D.; Morreale, B. D. Experimental Measurements and Equation of State Modeling of Liquid Densities for Long-Chain n-Alkanes at Pressures to 265 MPa and Temperatures to 523 K. Fluid Phase Equilib. 2011, 311, 17−24. (16) Sagdeev, D. I.; Fomina, M. G.; Mukhamedzyanov, G. K.; Abdulagatov, I. M. Experimental Study of the Density and Viscosity of n-Heptane at Temperatures from 298 K to 470 K and Pressure Up to 245 MPa. Int. J. Thermophys. 2013, 34, 1−33. (17) Dymond, J. H.; Robertson, J.; Isdale, J. D. Transport Properties of Nonelectrolyte Liquid MixturesIII. Viscosity Coefficients for nOctane, n-Dodecane, and Equimolar Mixtures of n-Octane + nDodecane and n-Hexane + n-Dodecane from 25 to 100°C at Pressures up to the Freezing Pressure or 500 MPa. Int. J. Thermophys. 1981, 2, 133−154. (18) Dymond, J. H.; Robertson, J.; Isdale, J. D. (p, ϱ, T) of Some Pure n-Alkanes and Binary Mixtures of n-Alkanes in the Range 298 to 373 K and 0.1 to 500 MPa. J. Chem. Thermodyn. 1982, 14, 51−59. (19) Rolo, L. I.; Caço, A. I.; Queimada, A. J.; Marrucho, I. M.; Coutinho, J. A. P. Surface Tension of Heptane, Decane, Hexadecane, Eicosane, and Some of Their Binary Mixtures. J. Chem. Eng. Data 2002, 47, 1442−1445. (20) Fröba, A. P.; Pellegrino, L. P.; Leipertz, A. Viscosity and Surface Tension of Saturated n-Pentane. Int. J. Thermophys. 2004, 25, 1323− 1337. (21) Queimada, A. J.; Caço, A. I.; Marrucho, I. M.; Coutinho, J. A. P. Surface Tension of Decane Binary and Ternary Mixtures with Eicosane, Docosane, and Tetracosane. J. Chem. Eng. Data 2005, 50, 1043−1046. (22) Queimada, A. J.; Marrucho, I. M.; Coutinho, J. A. P.; Stenby, E. H. Viscosity and Liquid Density of Asymmetric n-Alkane Mixtures: Measurement and Modeling. Int. J. Thermophys. 2005, 26, 47−61. (23) Mohsen-Nia, M.; Rasa, H.; Naghibi, S. F. Experimental and Theoretical Study of Surface Tension of n-Pentane, n-Heptane, and Some of Their Mixtures at Different Temperatures. J. Chem. Thermodyn. 2010, 42, 110−113.

captures the entire range of temperatures and properties satisfactorily, which is not surprising considering the different parametrization strategies and level of atomistic detail represented by all these force fields.

■

ASSOCIATED CONTENT

S Supporting Information *

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.jpcb.9b02840. Force field parameters for all species, representative nC12 VL surface tension and NEMD-PP viscosity starting conformations, NEMD-PP viscosity acceleration amplitude tests and calculation results, n-C12, n-C28, and nC12/n-C28 liquid densities, pure n-C12 and n-C28 mean squared end-to-end vectors, mean squared radii of gyration, C−C−C−C torsion probability distributions, maximum orientation relaxation times, VL equilibrium coexisting curves, critical temperatures, n-C12, n-C28, and n-C 12 /n-C28 calculated VL surface tensions and viscosities, and representative SX-70 wax model conformation with cumulative mole fraction probability distributions related to carbon number, SX-70 wax liquid densities, surface tension, and viscosities (PDF)

■

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. ORCID

Konstantinos D. Papavasileiou: 0000-0002-2322-7422 Loukas D. Peristeras: 0000-0002-3741-4077 Ioannis G. Economou: 0000-0002-2409-6831 Author Contributions

The manuscript was written through contributions of all authors. All the authors have given approval to the final version of the manuscript. Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS This publication is supported by the “Stavros Niarchos Foundation” Industrial Scholarship program and Scienomics SARL. The authors are thankful to the Greek Research & Technology Network (GRNET) for the computational time granted in the National High Performance Computing facilityARISunder project ID002043 and to the High Performance Computing center of Texas A&M University at Qatar for generous resource allocation.

■

REFERENCES

(1) Pedersen, K. S.; Christensen, P. L.; Shaikh, J. A. Phase Behavior of Petroleum Reservoir Fluids; 2nd ed.; CRC Press: Boca Raton, U.S.A., 2006. (2) Freund, M.; Mózes, G. Paraffin Products: Properties, Technologies, Applications; Elsevier Scientific Publishing Company: New York, U.S.A., 1982. (3) Jennings, D. W.; Weispfennig, K. Experimental Solubility Data of Various n-Alkane Waxes: Effects of Alkane Chain Length, Alkane Odd Versus Even Carbon Number Structures, and Solvent Chemistry on Solubility. Fluid Phase Equilib. 2005, 227, 27−35. (4) Nicholls, A.; Sharp, K. A.; Honig, B. Protein folding and association: Insights from the interfacial and thermodynamic properties of hydrocarbons. Proteins 1991, 11, 281−296. L


Article


(43) Nath, S. K.; Escobedo, F. A.; de Pablo, J. J. On the Simulation of Vapor-Liquid Equilibria for Alkanes. J. Chem. Phys. 1998, 108, 9905−9911. (44) Chilukoti, H. K.; Kikugawa, G.; Ohara, T. A Molecular Dynamics Study on Transport Properties and Structure at the LiquidVapor Interfaces of Alkanes. Int. J. Heat Mass Transfer 2013, 59, 144− 154. (45) Galicia-Andrés, E.; Medeiros, M. Vapour-Liquid Interfacial Properties of n-Alkanes. J. Mol. Liq. 2017, 248, 253−263. (46) Leonard, A. N.; Simmonett, A. C.; Pickard, F. C.; Huang, J.; Venable, R. M.; Klauda, J. B.; Brooks, B. R.; Pastor, R. W. Comparison of Additive and Polarizable Models with Explicit Treatment of LongRange Lennard-Jones Interactions Using Alkane Simulations. J. Chem. Theory Comput. 2018, 14, 948−958. (47) Klauda, J. B.; Venable, R. M.; Freites, J. A.; O’Connor, J. W.; Tobias, D. J.; Mondragon-Ramirez, C.; Vorobyov, I.; MacKerell, A. D.; Pastor, R. W. Update of the CHARMM All-Atom Additive Force Field for Lipids: Validation on Six Lipid Types. J. Phys. Chem. B 2010, 114, 7830−7843. (48) Lamoureux, G.; Roux, B. Modeling Induced Polarization with Classical Drude Oscillators: Theory and Molecular Dynamics Simulation Algorithm. J. Chem. Phys. 2003, 119, 3025−3039. (49) Isele-Holder, R. E.; Mitchell, W.; Ismail, A. E. Development and Application of a Particle-Particle Particle-Mesh Ewald Method for Dispersion Interactions. J. Chem. Phys. 2012, 137, 174107. (50) Wennberg, C. L.; Murtola, T.; Hess, B.; Lindahl, E. LennardJones Lattice Summation in Bilayer Simulations Has Critical Effects on Surface Tension and Lipid Properties. J. Chem. Theory Comput. 2013, 9, 3527−3537. (51) Isele-Holder, R. E.; Mitchell, W.; Hammond, J. R.; Kohlmeyer, A.; Ismail, A. E. Reconsidering Dispersion Potentials: Reduced Cutoffs in Mesh-Based Ewald Solvers Can Be Faster Than Truncation. J. Chem. Theory Comput. 2013, 9, 5412−5420. (52) Wennberg, C. L.; Murtola, T.; Páll, S.; Abraham, M. J.; Hess, B.; Lindahl, E. Direct-Space Corrections Enable Fast and Accurate Lorentz-Berthelot Combination Rule Lennard-Jones Lattice Summation. J. Chem. Theory Comput. 2015, 11, 5737−5746. (53) Morrow, B. H.; Harrison, J. A. Vapor-Liquid Equilibrium Simulations of Hydrocarbons Using Molecular Dynamics with LongRange Lennard-Jones Interactions. Energy Fuels 2019, 33, 848−858. (54) Siu, S. W. I.; Pluhackova, K.; Böckmann, R. A. Optimization of the OPLS-AA Force Field for Long Hydrocarbons. J. Chem. Theory Comput. 2012, 8, 1459−1470. (55) Best, R. B.; Zhu, X.; Shim, J.; Lopes, P. E. M.; Mittal, J.; Feig, M.; MacKerell, A. D. Optimization of the Additive CHARMM AllAtom Protein Force Field Targeting Improved Sampling of the Backbone ϕ, ψ and Side-Chain χ1 and χ2 Dihedral Angles. J. Chem. Theory Comput. 2012, 8, 3257−3273. (56) Edberg, R.; Morriss, G. P.; Evans, D. J. Rheology of n-Alkanes by Nonequilibrium Molecular Dynamics. J. Chem. Phys. 1987, 86, 4555−4570. (57) Morriss, G. P.; Daivis, P. J.; Evans, D. J. The Rheology of nAlkanes: Decane and Eicosane. J. Chem. Phys. 1991, 94, 7420−7433. (58) Daivis, P. J.; Evans, D. J.; Morriss, G. P. Computer-Simulation Study of the Comparative Rheology of Branched and Linear Alkanes. J. Chem. Phys. 1992, 97, 616−627. (59) Padilla, P.; Toxvaerd, S. Fluid n-Decane Undergoing Planar Couette Flow. J. Chem. Phys. 1992, 97, 7687−7694. (60) Berker, A.; Chynoweth, S.; Klomp, U. C.; Michopoulos, Y. Non-Equilibrium Molecular Dynamics (NEMD) Simulations and the Rheological Properties of Liquid n-Hexadecane. J. Chem. Soc., Faraday Trans. 1992, 88, 1719−1725. (61) Daivis, P. J.; Evans, D. J. Comparison of Constant-Pressure and Constant Volume Nonequilibrium Simulations of Sheared Model Decane. J. Chem. Phys. 1994, 100, 541−547. (62) Mundy, C. J.; Siepmann, J. I.; Klein, M. L. Decane Under Shear: A Molecular Dynamics Study Using Reversible NVT-SLLOD and NPT-SLLOD Algorithms. J. Chem. Phys. 1995, 103, 10192− 10200.

(24) Queimada, A. J.; Quiñones-Cisneros, S. E.; Marrucho, I. M.; Coutinho, J. A. P.; Stenby, E. H. Viscosity and Liquid Density of Asymmetric Hydrocarbon Mixtures. Int. J. Thermophys. 2003, 24, 1221−1239. (25) Koller, T. M.; Klein, T.; Giraudet, C.; Chen, J.; Kalantar, A.; van der Laan, G. P.; Rausch, M. H.; Fröba, A. P. Liquid Viscosity and Surface Tension of n-Dodecane, n-Octacosane, Their Mixtures, and a Wax between 323 and 573 K by Surface Light Scattering. J. Chem. Eng. Data 2017, 62, 3319−3333. (26) Alejandre, J.; Tildesley, D. J.; Chapela, G. A. Fluid Phase Equilibria Using Molecular Dynamics: The Surface Tension of Chlorine and Hexane. Mol. Phys. 1995, 85, 651−663. (27) Chen, B.; Siepmann, J. I.; Oh, K. J.; Klein, M. L. Simulating Vapor-Liquid Nucleation of n-Alkanes. J. Chem. Phys. 2002, 116, 4317−4329. (28) Nielsen, S. O.; Lopez, C. F.; Srinivas, G.; Klein, M. L. A Coarse Grain Model for n-Alkanes Parameterized from Surface Tension Data. J. Chem. Phys. 2003, 119, 7043−7049. (29) Pàmies, J. C.; McCabe, C.; Cummings, P. T.; Vega, L. F. Coexistence Densities of Methane and Propane by Canonical Molecular Dynamics and Gibbs Ensemble Monte Carlo Simulations. Mol. Simul. 2003, 29, 463−470. (30) Amat, M. A.; Rutledge, G. C. Liquid-Vapor Equilibria and Interfacial Properties of n-Alkanes and Perfluoroalkanes by Molecular Simulation. J. Chem. Phys. 2010, 132, 114704. (31) Harris, J. G. Liquid-Vapor Interfaces of Alkane Oligomers: Structure and Thermodynamics from Molecular Dynamics Simulations of Chemically Realistic Models. J. Phys. Chem. 1992, 96, 5077−5086. (32) Jorgensen, W. L.; Madura, J. D.; Swenson, C. J. Optimized Intermolecular Potential Functions for Liquid Hydrocarbons. J. Am. Chem. Soc. 1984, 106, 6638−6646. (33) Nicolas, J. P.; Smit, B. Molecular Dynamics Simulations of the Surface Tension of n-Hexane, n-Decane and n-Hexadecane. Mol. Phys. 2002, 100, 2471−2475. (34) Siepmann, J. I.; Karaborni, S.; Smit, B. Simulating the CriticalBehavior of Complex Fluids. Nature 1993, 365, 330−332. (35) Ibergay, C.; Ghoufi, A.; Goujon, F.; Ungerer, P.; Boutin, A.; Rousseau, B.; Malfreyt, P. Molecular Simulations of the n-Alkane Liquid-Vapor Interface: Interfacial Properties and Their Long Range Corrections. Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys. 2007, 75, 051602. (36) 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. (37) Gloor, G. J.; Jackson, G.; Blas, F. J.; de Miguel, E. Test-Area Simulation Method for the Direct Determination of the Interfacial Tension of Systems with Continuous or Discontinuous Potentials. J. Chem. Phys. 2005, 123, 134703. (38) Ismail, A. E.; Tsige, M.; Veld, P. J. I.; Grest, G. S. Surface Tension of Normal and Branched Alkanes. Mol. Phys. 2007, 105, 3155−3163. (39) Jorgensen, W. L.; Maxwell, D. S.; Tirado-Rives, J. Development and Testing of the OPLS All-Atom Force Field on Conformational Energetics and Properties of Organic Liquids. J. Am. Chem. Soc. 1996, 118, 11225−11236. (40) Smith, G. D.; Borodin, O.; Bedrov, D. A Revised Quantum Chemistry-Based Potential for Poly(Ethylene Oxide) and Its Oligomers in Aqueous Solution. J. Comput. Chem. 2002, 23, 1480− 1488. (41) Borodin, O.; Smith, G. D. Development of Many−Body Polarizable Force Fields for Li-Battery Components: 1. Ether, Alkane, and Carbonate-Based Solvents. J. Phys. Chem. B 2006, 110, 6279− 6292. (42) Mendoza, F. N.; López-Rendón, R.; López-Lemus, J.; Cruz, J.; Alejandre, J. Surface Tension of Hydrocarbon Chains at the LiquidVapour Interface. Mol. Phys. 2008, 106, 1055−1059. M


Article

The Journal of Physical Chemistry B (63) Mundy, C. J.; Siepmann, J. I.; Klein, M. L. Erratum: Decane under shear: A molecular dynamics study using reversible NVTSLLOD and NPT-SLLOD algorithms [J. Chem. Phys. 103, 10192 (1995)]. J. Chem. Phys. 1996, 104, 7797. (64) Cui, S. T.; Gupta, S. A.; Cummings, P. T.; Cochran, H. D. Molecular Dynamics Simulations of the Rheology of Normal Decane, Hexadecane, and Tetracosane. J. Chem. Phys. 1996, 105, 1214−1220. (65) McCabe, C.; Manke, C. W.; Cummings, P. T. Predicting the Newtonian Viscosity of Complex Fluids from High Strain Rate Molecular Simulations. J. Chem. Phys. 2002, 116, 3339−3342. (66) Baig, C.; Edwards, B. J.; Keffer, D. J.; Cochran, H. D. Rheological and Structural Studies of Liquid Decane, Hexadecane, and Tetracosane under Planar Elongational Flow Using Nonequilibrium Molecular-Dynamics Simulations. J. Chem. Phys. 2005, 122, 184906. (67) Marechal, G.; Ryckaert, J.-P. Atomic Versus Molecular Description of Transport Properties in Polyatomic Fluids: n-Butane as an Illustration. Chem. Phys. Lett. 1983, 101, 548−554. (68) Mundy, C. J.; Siepmann, J. I.; Klein, M. L. Calculation of the Shear Viscosity of Decane Using a Reversible Multiple Time-Step Algorithm. J. Chem. Phys. 1995, 102, 3376−3380. (69) Daivis, P. J.; Evans, D. J. Transport-Coefficients of Liquid Butane near the Boiling-Point by Equilibrium Molecular-Dynamics. J. Chem. Phys. 1995, 103, 4261−4265. (70) Mondello, M.; Grest, G. S. Viscosity Calculations of n-Alkanes by Equilibrium Molecular Dynamics. J. Chem. Phys. 1997, 106, 9327− 9336. (71) Kubo, R.; Toda, M.; Hashitsume, N. Statistical Physics II: Nonequilibrium Statistical Mechanics; Springer Berlin Heidelberg: Berlin, Heidelberg, 1985. (72) Allen, M. P.; Tildesley, D. J. Computer Simulation of Liquids; Oxford University Press: Oxford U.K., 1989. (73) Evans, D. J.; Morriss, G. P. Statistical Mechanics of Nonequilibrium Liquids; Academic Press: London U.K., 1990. (74) Hess, B. Determining the Shear Viscosity of Model Liquids from Molecular Dynamics Simulations. J. Chem. Phys. 2002, 116, 209−217. (75) Mundy, C. J.; Balasubramanian, S.; Bagchi, K.; Siepmann, J. I.; Klein, M. L. Equilibrium and Non-Equilibrium Simulation Studies of Fluid Alkanes in Bulk and at Interfaces. Faraday Discuss. 1996, 104, 17−36. (76) Allen, W.; Rowley, R. L. Predicting the Viscosity of Alkanes Using Nonequilibrium Molecular Dynamics: Evaluation of Intermolecular Potential Models. J. Chem. Phys. 1997, 106, 10273−10281. (77) Ryckaert, J.-P.; Bellemans, A. Molecular-Dynamics of Liquid Alkanes. Faraday Discuss. 1978, 66, 95−106. (78) Dysthe, D. K.; Fuchs, A. H.; Rousseau, B. Fluid Transport Properties by Equilibrium Molecular Dynamics. III. Evaluation of United Atom Interaction Potential Models for Pure Alkanes. J. Chem. Phys. 2000, 112, 7581−7590. (79) Gordon, P. A. Development of Intermolecular Potentials for Predicting Transport Properties of Hydrocarbons. J. Chem. Phys. 2006, 125, 014504. (80) Nieto-Draghi, C.; Ungerer, P.; Rousseau, B. Optimization of the Anisotropic United Atoms Intermolecular Potential for n-Alkanes: Improvement of Transport Properties. J. Chem. Phys. 2006, 125, 044517. (81) Payal, R. S.; Balasubramanian, S.; Rudra, I.; Tandon, K.; Mahlke, I.; Doyle, D.; Cracknell, R. Shear Viscosity of Linear Alkanes Through Molecular Simulations: Quantitative Tests for n-Decane and n-Hexadecane. Mol. Simul. 2012, 38, 1234−1241. (82) Tobias, D. J.; Tu, K.; Klein, M. L. Assessment of All-Atom Potentials for Modeling Membranes: Molecular Dynamics Simulations of Solid and Liquid Alkanes and Crystals of Phospholipid Fragments. J. Chim. Phys. Phys. Chim. Biol. 1997, 94, 1482−1502. (83) Makrodimitri, Z. Α.; Heller, A.; Koller, T. M.; Rausch, M. H.; Fleys, M. S. H.; Bos, A. N. R.; van der Laan, G. P.; Fröba, A. P.; Economou, I. G. Viscosity of Heavy n-Alkanes and Diffusion of Gases

Therein Based on Molecular Dynamics Simulations and Empirical Correlations. J. Chem. Thermodyn. 2015, 91, 101−107. (84) Ewen, J.; Gattinoni, C.; Thakkar, F.; Morgan, N.; Spikes, H.; Dini, D. A Comparison of Classical Force-Fields for Molecular Dynamics Simulations of Lubricants. Materials 2016, 9, 651. (85) Dickson, C. J.; Madej, B. D.; Skjevik, Å. A.; Betz, R. M.; Teigen, K.; Gould, I. R.; Walker, R. C. Lipid14: The Amber Lipid Force Field. J. Chem. Theory Comput. 2014, 10, 865−879. (86) Marrink, S. J.; Risselada, H. J.; Yefimov, S.; Tieleman, D. P.; de Vries, A. H. The Martini Force Field: Coarse Grained Model for Biomolecular Simulations. J. Phys. Chem. B 2007, 111, 7812−7824. (87) Moultos, O. A.; Tsimpanogiannis, I. N.; Panagiotopoulos, A. Z.; Trusler, J. P. M.; Economou, I. G. Atomistic Molecular Dynamics Simulations of Carbon Dioxide Diffusivity in n-Hexane, n-Decane, nHexadecane, Cyclohexane, and Squalane. J. Phys. Chem. B 2016, 120, 12890−12900. (88) Cornell, W. D.; Cieplak, P.; Bayly, C. I.; Gould, I. R.; Merz, K. M.; Ferguson, D. M.; Spellmeyer, D. C.; Fox, T.; Caldwell, J. W.; Kollman, P. A. A Second Generation Force Field for the Simulation of Proteins, Nucleic Acids, and Organic Molecules. J. Am. Chem. Soc. 1995, 117, 5179−5197. (89) Wang, J.; Wolf, R. M.; Caldwell, J. W.; Kollman, P. A.; Case, D. A. Development and Testing of a General Amber Force Field. J. Comput. Chem. 2004, 25, 1157−1174. (90) Papavasileiou, K. D.; Moultos, O. A.; Economou, I. G. Predictions of Water/Oil Interfacial Tension at Elevated Temperatures and Pressures: A Molecular Dynamics Simulation Study with Biomolecular Force Fields. Fluid Phase Equilib. 2018, 476, 30−38. (91) Marrink, S. J.; de Vries, A. H.; Mark, A. E. Coarse Grained Model for Semiquantitative Lipid Simulations. J. Phys. Chem. B 2004, 108, 750−760. (92) Baoukina, S.; Marrink, S. J.; Tieleman, D. P. Molecular Structure of Membrane Tethers. Biophys. J. 2012, 102, 1866−1871. (93) Rossi, G.; Monticelli, L.; Puisto, S. R.; Vattulainen, I.; AlaNissila, T. Coarse-Graining Polymers with the Martini Force-Field: Polystyrene as a Benchmark Case. Soft Matter 2011, 7, 698−708. (94) Baaden, M.; Marrink, S. J. Coarse-Grain Modelling of ProteinProtein Interactions. Curr. Opin. Struct. Biol. 2013, 23, 878−886. (95) Baoukina, S.; Monticelli, L.; Tieleman, D. P. Interaction of Pristine and Functionalized Carbon Nanotubes with Lipid Membranes. J. Phys. Chem. B 2013, 117, 12113−12123. (96) Scienomics. Maps Platform, 4.2; Paris, France. http://www. scienomics.com/ (accessed September 28, 2018). (97) Ramos, J.; Peristeras, L. D.; Theodorou, D. N. Monte Carlo Simulation of Short Chain Branched Polyolefins in the Molten State. Macromolecules 2007, 40, 9640−9650. (98) Bekker, H.; Berendsen, H. J. C.; Dijkstra, E. J.; Achterop, S.; von Drunen, R.; van der Spoel, D.; Sijbers, A.; Keegstra, H.; Reitsma, B.; Renardus, M. K. R., Gromacsa Parallel Computer for MolecularDynamics Simulations. In 4th International Conference on Computational Physics (PC 92); de Groot, R. A., Nadrchal, J., Eds.; World Scientific Publishing, 1993. (99) Berendsen, H. J. C.; van der Spoel, D.; van Drunen, R. Gromacs: A Message-Passing Parallel Molecular Dynamics Implementation. Comput. Phys. Commun. 1995, 91, 43−56. (100) Lindahl, E.; Hess, B.; van der Spoel, D. Gromacs 3.0: A Package for Molecular Simulation and Trajectory Analysis. J. Mol. Model. 2001, 7, 306−317. (101) van der Spoel, D.; Lindahl, E.; Hess, B.; Groenhof, G.; Mark, A. E.; Berendsen, H. J. C. Gromacs: Fast, Flexible, and Free. J. Comput. Chem. 2005, 26, 1701−1718. (102) Hess, B.; Kutzner, C.; van der Spoel, D.; Lindahl, E. Gromacs 4: Algorithms for Highly Efficient, Load-Balanced, and Scalable Molecular Simulation. J. Chem. Theory Comput. 2008, 4, 435−447. (103) Abraham, M. J.; Murtola, T.; Schulz, R.; Páll, S.; Smith, J. C.; Hess, B.; Lindahl, E. Gromacs: High Performance Molecular Simulations through Multi-Level Parallelism from Laptops to Supercomputers. SoftwareX 2015, 1−2, 19−25. N


Article

The Journal of Physical Chemistry B (104) van der Spoel, D.; van Maaren, P. J. The Origin of Layer Structure Artifacts in Simulations of Liquid Water. J. Chem. Theory Comput. 2006, 2, 1−11. (105) Darden, T.; York, D.; Pedersen, L. Particle Mesh Ewald: An N·Log(N) Method for Ewald Sums in Large Systems. J. Chem. Phys. 1993, 98, 10089−10092. (106) Salomon-Ferrer, R.; Götz, A. W.; Poole, D.; Le Grand, S.; Walker, R. C. Routine Microsecond Molecular Dynamics Simulations with Amber on GPUs. 2. Explicit Solvent Particle Mesh Ewald. J. Chem. Theory Comput. 2013, 9, 3878−3888. (107) Hess, B.; Bekker, H.; Berendsen, H. J. C.; Fraaije, J. G. E. M. LINCS: A Linear Constraint Solver for Molecular Simulations. J. Chem. Theory Comput. 1997, 18, 1463−1472. (108) de Jong, D. H.; Baoukina, S.; Ingólfsson, H. I.; Marrink, S. J. Martini Straight: Boosting Performance Using a Shorter Cutoff and Gpus. Comput. Phys. Commun. 2016, 199, 1−7. (109) Hoover, W. G. Canonical Dynamics - Equilibrium PhaseSpace Distributions. Phys. Rev. A: At., Mol., Opt. Phys. 1985, 31, 1695−1697. (110) Parrinello, M.; Rahman, A. Polymorphic Transitions in Single Crystals: A New Molecular Dynamics Method. J. Appl. Phys. 1981, 52, 7182−7190. (111) Flyvbjerg, H.; Petersen, H. G. Error-Estimates on Averages of Correlated Data. J. Chem. Phys. 1989, 91, 461−466. (112) Vega, C.; de Miguel, E. Surface Tension of the Most Popular Models of Water by Using the Test-Area Simulation Method. J. Chem. Phys. 2007, 126, 154707. (113) Goujon, F.; Malfreyt, P.; Boutin, A.; Fuchs, A. H. Direct Monte Carlo Simulations of the Equilibrium Properties of n-Pentane Liquid-Vapor Interface. J. Chem. Phys. 2002, 116, 8106−8117. (114) Janeček, J. Long Range Corrections in Inhomogeneous Simulations. J. Phys. Chem. B 2006, 110, 6264−6269. (115) Lundberg, L.; Edholm, O. Dispersion Corrections to the Surface Tension at Planar Surfaces. J. Chem. Theory Comput. 2016, 12, 4025−4032. (116) Ghoufi, A.; Malfreyt, P.; Tildesley, D. J. Computer Modelling of the Surface Tension of the Gas-Liquid and Liquid-Liquid Interface. Chem. Soc. Rev. 2016, 45, 1387−1409. (117) Fischer, N. M.; van Maaren, P. J.; Ditz, J. C.; Yildirim, A.; van der Spoel, D. Properties of Organic Liquids When Simulated with Long-Range Lennard-Jones Interactions. J. Chem. Theory Comput. 2015, 11, 2938−2944. (118) Zhang, Y.; Otani, A.; Maginn, E. J. Reliable Viscosity Calculation from Equilibrium Molecular Dynamics Simulations: A Time Decomposition Method. J. Chem. Theory Comput. 2015, 11, 3537−3546. (119) Baig, C.; Edwards, B. J.; Keffer, D. J.; Cochran, H. D. A Proper Approach for Nonequilibrium Molecular Dynamics Simulations of Planar Elongational Flow. J. Chem. Phys. 2005, 122, 114103. (120) Baig, C.; Edwards, B. J.; Keffer, D. J.; Cochran, H. D.; Harmandaris, V. A. Rheological and Structural Studies of Linear Polyethylene Melts under Planar Elongational Flow Using Nonequilibrium Molecular Dynamics Simulations. J. Chem. Phys. 2006, 124, 084902. (121) Ladd, A. J. C. Equations of Motion for Non-Equilibrium Molecular-Dynamics Simulations of Viscous-Flow in Molecular Fluids. Mol. Phys. 1984, 53, 459−463. (122) Evans, D. J.; Morriss, O. P. Non-Newtonian MolecularDynamics. Comput. Phys. Rep. 1984, 1, 297−343. (123) Mü ller-Plathe, F. Reversing the Perturbation in Nonequilibrium Molecular Dynamics: An Easy Way to Calculate the Shear Viscosity of Fluids. Phys. Rev. E: Stat. Phys., Plasmas, Fluids, Relat. Interdiscip. Top. 1999, 59, 4894−4898. (124) Gosling, E. M.; Mcdonald, I. R.; Singer, K. Calculation by Molecular-Dynamics of Shear Viscosity of a Simple Fluid. Mol. Phys. 1973, 26, 1475−1484. (125) Ciccotti, G.; Jacucci, G.; Mcdonald, I. R. Transport Properties of Molten Alkali-Halides. Phys. Rev. A: At., Mol., Opt. Phys. 1976, 13, 426−436.

(126) Zhao, L.; Wang, X.; Wang, L.; Sun, H. Prediction of Shear Viscosities Using Periodic Perturbation Method and OPLS Force Field. Fluid Phase Equilib. 2007, 260, 212−217. (127) Zhao, L. F.; Cheng, T.; Sun, H. On the Accuracy of Predicting Shear Viscosity of Molecular Liquids Using the Periodic Perturbation Method. J. Chem. Phys. 2008, 129, 144501. (128) Hansen, J.-P.; McDonald, I. R. Theory of Simple Liquids; 3rd ed.; Elsevier Academic Press: London U.K., 2006. (129) Vasquez, V. R.; Macedo, E. A.; Zabaloy, M. S. Lennard-Jones Viscosities in Wide Ranges of Temperature and Density: Fast Calculations Using a Steady-State Periodic Perturbation Method. Int. J. Thermophys. 2004, 25, 1799−1818. (130) Lemmon, E. W.; Huber, M. L. Thermodynamic Properties of n-Dodecane. Energy Fuels 2004, 18, 960−967. (131) Maerzke, K. A.; Siepmann, J. I. Transferable Potentials for Phase Equilibria−Coarse-Grain Description for Linear Alkanes. J. Phys. Chem. B 2011, 115, 3452−3465. (132) Cao, F.; Sun, H. Transferability and Nonbond Functional Form of Coarse Grained Force Field - Tested on Linear Alkanes. J. Chem. Theory Comput. 2015, 11, 4760−4769. (133) Dutour, S.; Lagourette, B.; Daridon, J. L. High-Pressure Speed of Sound, Density and Compressibility of Heavy Normal Paraffins: C28H58 and C36H74. J. Chem. Thermodyn. 2002, 34, 475−484. (134) Vrabec, J.; Kedia, G. K.; Fuchs, G.; Hasse, H. Comprehensive Study of the Vapour-Liquid Coexistence of the Truncated and Shifted Lennard-Jones Fluid Including Planar and Spherical Interface Properties. Mol. Phys. 2006, 104, 1509−1527. (135) Goujon, F.; Malfreyt, P.; Tildesley, D. J. The Gas-Liquid Surface Tension of Argon: A Reconciliation between Experiment and Simulation. J. Chem. Phys. 2014, 140, 244710. (136) Goujon, F.; Ghoufi, A.; Malfreyt, P.; Tildesley, D. J. Can We Approach the Gas-Liquid Critical Point Using Slab Simulations of Two Coexisting Phases? J. Chem. Phys. 2016, 145, 124702. (137) Binder, K. Monte Carlo Calculation of the Surface Tension for Two- and Three-Dimensional Lattice-Gas Models. Phys. Rev. A: At., Mol., Opt. Phys. 1982, 25, 1699−1709. (138) Potoff, J. J.; Panagiotopoulos, A. Z. Surface Tension of the Three-Dimensional Lennard-Jones Fluid from Histogram-Reweighting Monte Carlo Simulations. J. Chem. Phys. 2000, 112, 6411−6415. (139) Ndao, M.; Devémy, J.; Ghoufi, A.; Malfreyt, P. CoarseGraining the Liquid−Liquid Interfaces with the MARTINI Force Field: How Is the Interfacial Tension Reproduced? J. Chem. Theory Comput. 2015, 11, 3818−3828. (140) Ramazani, A.; Mandal, T.; Larson, R. G. Modeling the Hydrophobicity of Nanoparticles and Their Interaction with Lipids and Proteins. Langmuir 2016, 32, 13084−13094. (141) Mondello, M.; Grest, G. S.; Webb, E. B.; Peczak, P. Dynamics of n-Alkanes: Comparison to Rouse Model. J. Chem. Phys. 1998, 109, 798−805. (142) Marrink, S. J.; Tieleman, D. P. Perspective on the Martini Model. Chem. Soc. Rev. 2013, 42, 6801−6822. (143) Allal, A.; Boned, C.; Baylaucq, A. Free-Volume Viscosity Model for Fluids in the Dense and Gaseous States. Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys. 2001, 64, 011203.

O


Molecular Dynamics Simulation of Pure n-Alkanes and Their Mixtures

Recommend Documents