The Water-Alkane Interface at Various NaCl Salt Concentrations: A Molecular Dynamics Study of the Readily Available Force Fields

In this study, classical molecular dynamic simulations have been used to examine the molecular properties of the water-alkane interface at various NaCl salt concentrations (up to 3.0 mol/kg). A variety of different force field combinations have been compared against experimental surface/interfacial tension values for the water-vapour, decane-vapour and water-decane interfaces. Six different force fields for water (SPC, SPC/E, TIP3P, TIP3Pcharmm, TIP4P & TIP4P2005), and three further force fields for alkane (TraPPE-UA, CGenFF & OPLS) have been compared to experimental data. CGenFF, OPLS-AA and TraPPE-UA all accurately reproduce the interfacial properties of decane. The TIP4P2005 (four-point) water model is shown to be the most accurate water model for predicting the interfacial properties of water. The SPC/E water model is the best three-point parameterisation of water for this purpose. The CGenFF and TraPPE parameterisations of oil accurately reproduce the interfacial tension with water using either the TIP4P2005 or SPC/E water model. The salinity dependence on surface/interfacial tension is accurately captured using the Smith & Dang parameterisation of NaCl. We observe that the Smith & Dang model slightly overestimates the surface/interfacial tensions at higher salinities (>1.5 mol/kg). This is ascribed to an overestimation of the ion exclusion at the interface.


Results
The Alkane-Vapour Interface. Figure 1 presents the calculated surface tension values of decane over a range of temperatures. The results show that the three tested force fields accurately reproduce the overall trend of surface tension at various temperatures.
Whilst we have examined decane solely, the considered force fields are widely transferable. Both OPLS-AA and TraPPE-UA have previously been shown to accurately predict the surface tension values of ethane through to hexadecane 33,34 . To the author's knowledge, there have been no surface tension calculations of alkanes using the CGenFF force field.
Both OPLS-AA and TraPPE-UA systematically underestimate the surface tension of decane at all temperatures, OPLS more-so than TraPPE-UA. CGenFF overestimates surface tension at lower temperatures, and underestimates surface tension at higher temperatures. Our results concur with the work of Ismail et al., who show that OPLS-AA consistently underestimates the surface tension of various saturated n-alkanes 33 and with Mendoza et al., who present that TraPPE-UA is generally in good agreement with experimental data 34  Previous work by Nicolas & Smit has shown that the accurate calculation of surface tension for saturated alkanes is linearly proportional to the liquid density predicted in the simulations 36 . Values for the liquid and vapour density of decane have been calculated by fitting the interfacial density profile to a tanh function (see Equation 6 in Methods). The liquid density (ρ l ) of decane at 293.15 K is presented in Table 1, along with values for the interfacial width (Δ). Note that the vapour density (ρ v ) at 293.15 K is calculated to be zero, and is therefore not presented.
Our results show that both TraPPE-UA and OPLS-AA marginally underestimate the liquid density of decane at 293.15 K, and consequently, both underestimate the surface tension. Conversely, CGenFF overestimates the liquid density, and therefore the surface tension of decane. The accurate reproduction of liquid-vapour densities depends on the treatment of long-range dispersion forces 37 . Ismail et al. show that Lennard-Jones cutoffs greater than 16 Å are required to yield results in agreement with explicit long range methods (evaluated in reciprocal space) 33 , whilst the simulations of Mendoza et al. and Lopez et al. also show the utility of using Ewald-summations to treat these long range dispersion forces 34,38 . Note that our simulations satisfy both these criteria. Figure 2 presents the interfacial density profiles (both intrinsic and non-intrinsic) of the decane-vapour interface at 293.15 K. The density profiles show that all three force fields act remarkably similarly, even though CGenFF and OPLS-AA are all-atom representations of decane, and TraPPE-UA is a united-atom representation. The intrinsic density profiles present four discrete peaks of decane before the density becomes bulk-like.
The Water-Vapour Interface. The calculated surface tension values of pure water are presented in Table 2.
It can be seen that almost all water models underestimate the experimentally observed surface tension of water by at least 15%. The only exception to this is TIP4P2005, which underestimates the surface tension of water by 7%. Our results compare favourably to the study of Vega & De Miguel 25 . Like their work, we find that the accurate reproduction of surface tension follows the trend: The TIP3Pc model, which introduces LJ sites on hydrogen atoms, performs better than the original TIP3P parameterisation, but still underestimates the surface tension of water by 23%. We observe that, generally, the four-point water models perform better than the three-point models, in keeping with the observations of Vega & De Miguel 25 . Of all the three-point water models, SPC/E most accurately reproduces the experimental surface tension of water, whilst of the four-point water models, the TIP4P2005 parameterisation performs the best.
The density profiles of TIP4P2005 and SPC/E water are presented in Fig. 3. Analysing the intrinsic density profiles, it can be seen that water forms three separate density peaks at the water-vapour interface. This information is lost due to the capillary wave fluctuations present in the simulation when examining the non-intrinsic density profile solely. Also notable, is that most water models appear to overestimate the interfacial width of water   13 , see Table 2. The TIP4P2005 and SPC/E water model most accurately predict the experimentally observed interfacial width of water. Overall, the TIP4P2005 is best placed to model water-vapour interfaces. Of the three-point water models, SPC/E is the most favourable.
The Water-Vapour Interface at Various NaCl Concentrations. Proceeding, we have chosen three water models to examine the trends in surface tension as a function of salinity. The TIP4P2005 and SPC/E water models were examined as they represent the best four-point and three-point water models respectively. The TIP3Pc water model was also examined, as it is frequently simulated in conjunction with secondary organic phases 6,7 . The variation of surface tension with NaCl concentration is presented in Fig. 4. Firstly, we note that the overall trend of surface tension variation with NaCl concentration is well captured by the classical MD simulations. The simulations are able to accurately predict the experimental surface tensions of   13 . NaCl solutions up to 1.5 mol/kg, however, the simulations slightly overestimate the surface tension of NaCl electrolyte at higher concentrations, deviating more drastically as concentrations increase beyond 2.0 mol/kg. Similar behaviour was noted by D' Auria & Tobias, who observed that the surface tension of the TIP3P water model was overestimated when used in conjunction with the Smith & Dang parameterisation at NaCl concentrations of 1.2 mol/kg and 6.17 mol/kg 39 . Secondly, we observe that the simulated surface tensions of all three water models agree with each other at each concentration. This implies that the accurate modelling of interfaces containing electrolytes depends more so on the parameterisation of the ion, rather than the water model. This conclusion concurs with the work of Neyt et al., who presented that the salinity dependence of the electrolyte-vapour interface varied primarily with ion parameterisation, rather than the water model, whilst the utilised water model sets the accuracy of the surface tension calculation without additional NaCl 40 . Notably, Neyt et al. presented that the TIP4P2005 water model used in conjunction with the OPLS parameterisation of NaCl, typically returned surface tensions within 0.2% of the experimental trend 40 . Table 3 presents the interfacial tension values of the decane-water interface, without salts, for various force field combinations as predicted by the MD simulations. We have chiefly examined TraPPE-UA and CGenFF due to their compatibility (i.e. consistent use of Lennard-Jones mixing rules) with ClayFF; which makes these force fields more suitable for use in complex three-phase simulations. Priority has also been placed on testing the TIP4P2005 and SPC/E water models, as these most accurately model the interfacial properties of pure water. TIP3Pc has also been examined as it is the default three-point water model used in conjunction with CGenFF 8 .

The Water-Alkane Interface at Various NaCl Concentrations.
In addition to the usual TIP4P2005 water model interacting with TraPPE-UA, we have examined the parameters as described by Ashbaugh et al. 41 , denoted TIP4P2005*. The TIP4P2005* water model uses custom Lennard-Jones interactions between the oxygen atoms of water and the TraPPE-UA beads. These custom interactions have been finely tuned to accurately reproduce the hydration free energy of alkanes 41 . All other parameters of the TIP4P2005* water model are identical to the original TIP4P2005 parameterisation.
We note that the MD simulations are able to predict the interfacial tension of the water-decane interface within 10% of the empirically observed values. The simulations of this interface are notably more accurate than the simulations of the water-vapour interface. The combination of different force fields appears more robust for the water-decane interface, compared to the water-vapour interface. This is likely due to the absence of a vapour phase in the simulations, which when inaccurately modelled (in particular, the vapour density), can lead to deviations in the predicted surface tension, as noted by Nicolas & Smit 36 .    sets agree with each other at each NaCl concentration. This again highlights that the ion parameterisation may be more important for the accurate modelling of the water-vapour/water-decane interface.

Discussion
The presented results show the quality with which MD simulations can capture the bulk properties and structural properties of the liquid-liquid interface. One key advantage of such simulations, compared to many experiments, is the unparalleled resolution with which one can interpret the bulk phenomena on display. In this instance, we can use the MD simulations to further our understanding of why the surface tension at the water-vapour/ water-alkane interface increases with increasing NaCl concentration. Typically, this behaviour is explained due to an exclusion of ions from the liquid-vapour interface 3 . This directly alters the surface tension between the two phases as described by the Gibbs adsorption isotherm: where Γ i is the surface excess and μ i is the chemical potential of species i∈{ Na, Cl }. In this study, the surface excess of Na + and Cl − ions across the electrolyte-water interface has been calculated as: . Z Gibbs is the location of the Gibbs dividing plane, which is priorly calculated by minimising Γ water , such that: Note that equation 3 has been used to centre the non-intrinsic density profiles presented in the results section. Figure 6 presents the negative total surface excess of ions (Γ = Γ + Γ NaCl Na Cl ) at the water-vapour interface calculated from all simulations as a function of NaCl concentration. We observe that, compared to the experimental data of Ali et al., the simulations marginally overestimate the negative total surface excess of NaCl, −Γ NaCl , especially at higher concentrations 43 . In turn, this elucidates why the surface tension values are   40,44 . Both studies concluded that current models of polarizable force fields (including both polarizable water and/or ions), are not mature enough to capture the salinity dependence on surface tension. These studies tested both the Drude oscillator model and the electronic continuum correction model (whereby classical point charges are shifted by a screening factor to account for polarization). Work by Jiang & Panagiotopolous further show that theelectronic continuum correction model fails to accurately predict the interfacial properties of electrolytes with polarizable models 45 . More recently, work by Jiang et al. presents that the BK3 polarizable model is able to accurately capture the interfacial properties of NaCl electrolyte solutions 46 . However, it is worth noting that such polarizable models typically run 5-10 times slower compared to the classical parameterizations used in this study, and that the implementation of the BK3 model currently only runs on a single processor using GROMACS (the software suite used in this study).
Also noteworthy is that the variation in interfacial tension predicted at the interface between water and decane is remarkably similar to that at the water-vapour interface. This is apparent in Fig. 7, which compares the change in surface tension at the water-vapour interface against the change in interfacial tension at the water-alkane interface (using the TIP4P2005 water model and the TraPPE-UA oil model). The overall change in interfacial tension is therefore due to the properties primarily within the water phase. This result is somewhat unsurprising as the dielectric properties of decane and water vapour are much similar when compared to bulk liquid water. The surface tensions and interfacial tensions scale linearly as the ions are excluded from their respected interfaces in equal value, and the ions are not miscible in the vapour/decane phase respectively.
To conclude; in this study, the interfacial properties of water and decane have been examined at various NaCl concentrations using classical molecular dynamics simulations. By choosing an appropriate set of interaction parameters (force fields), one can obtain remarkable agreement between model and experimental observation. In particular, the TIP4P2005 water model is best placed to examine the interfacial properties of water. The SPC/E parameterisation is the best three-point water model to interpret the interfacial behaviour of water. CGenFF, OPLS-AA and TraPPE-UA all accurately reproduce the interfacial properties of decane. In combination, CGenFF and TraPPE-UA are compatible with various water models, and are able to accurately predict the interfacial tension of the water-decane interface. The salinity dependence on surface/interfacial tension is well captured using the Smith & Dang parameterisation of NaCl. We observe that the model slightly overestimates the surface/interfacial tensions at higher salinities. This is due to an overestimation of the ion exclusion at the interface. The simulations further suggest that the salinity dependence on surface/interfacial tension is dictated by the parameterisation of the salt ions. Future work will examine this hypothesis, and will model different ion parameterizations. Future work will also examine the role of divalent cations at the water-vapour/water-alkane interface, whereby, charge inversion may play a determining role on the behaviour of the liquid-liquid interface.

Methods
The datasets generated and analysed during the current study are available from the corresponding authors on reasonable request.
All simulations were calculated using GROMACS 5.1.4 47 with an electrostatic and van der Waals cutoff radii of 2.0 nm. Long range electrostatics were calculated using a Particle-Mesh-Ewald (PME) summation with grid spacings of 0.1 nm. The PME summation used a spline interpolation order of 4, and long-range electrostatic interactions were accurate to within 99.999%. All simulations were initialised with an energy minimisation calculation to minimise any unphysical atomic overlaps. This was achieved using a steepest descents algorithm, which was terminated once the maximum force on any one atom was less than 100 kJ/(mol nm).
All simulations were subsequently equilibrated for 1 ns. Simulations of liquid-vapour phase surface tensions were calculated in the canonical ensemble (constant particle number -N, constant volume -V and constant temperature T) at 293.15 K, using a V-rescale thermostat set to rescale system temperatures every 0.5 ps. Liquid-liquid interfacial tension simulations were calculated in the isothermal-isobaric ensemble (constant particle number -N, constant volume -V and constant pressure -P) at 293.15 K and 1.01325 bar. NPT simulations were equilibrated using a V-rescale thermostat with a temperature coupling constant of 0.5 ps. Pressure coupling was achieved using a Berendsen barostat, with a pressure coupling constant of 1 ps.
Following equilibration, all simulations were run for a production period of 10 ns. During the production period, all simulations used a Nosé-Hoover thermostat with a temperature coupling constant of 0.5 ps and a Nosé-Hoover chain length of 1. NPT simulations used an isotropic Parrinello-Rahman barostat during the production period, with a pressure coupling constant of 1 ps. Unless otherwise stated, all production simulations have been calculated at 293.15 K and 1.01325 bar.
Three different force field parameterisations for decane have been tested in this study: TraPPE-UA 30 , OPLS-AA 31,48 and CGenFF [8][9][10] . TraPPE-UA is a united atom force field, where CH 3 -and -CH 2 -groups are modelled as soft spheres. OPLS-AA and CGenFF are both all atom models, where carbon and hydrogen atoms are modelled explicitly. Six different water models were examined in the water-vapour interface: SPC 27  System Setup. Three different systems are presented in this study. In the first, a 5 × 5 × 5 nm 3 box of 392 decane molecules is inserted into 5 × 5 × 20 nm 3 simulation box. In the second, a 5× 5× 5 nm 3 box of 4139 water molecules is inserted into 5 × 5× 20 nm 3 simulation box. The number of decane and water molecules in each film is calculated to match the bulk density of each solvent at 293.15 K and 1.01325 bar, 0.727 g/cm 3 for decane and 0.9982 g/cm 3 for water respectively 51 . In the third system, the decane and water films are combined in a 5 × 5× 10 nm 3 simulation box. Each system was generated using the Packmol software package 52 . Systems involving water were further examined at various NaCl concentrations. Systems were setup in terms of NaCl molal concentration, up to a maximum of 3.00 mol/kg, in increments of 0.50 mol/kg. This was achieved by replacing water molecules with the relevant amount of sodium and chloride ions. The number of water molecules and ions present in each simulation is presented in Table 4. Additionally, 0.20 mol/kg NaCl solution was examined. Whilst simulation results are presented in molal concentration (mol/kg), experimental results are often presented in terms of molar concentration (mol/L). The conversion from molar concentration to molal concentration is presented in Table 4, using data extrapolated from the CRC Handbook of Chemistry & Physics 53 .
Analyses. Thermodynamic data from each simulation were output every 1 ps. Final values for thermodynamic quantities were averaged over all 10 ns, and errors were calculated using a block-averaging method, with each block averaging over a 1 ns timeframe. In all figures, error bars are presented to ±2 standard errors of the mean (a confidence interval of 95%).
The surface tension across an interface has been calculated using the diagonal components of the local pressure tensor: where L z is the length of the simulation in z (the direction normal to the interface) and p N (z) & p T (z) represent the normal and tangential components of the pressure tensor with respect to the interface:  The diagonal components of the local pressure tensor (p xx , p yy & p zz ) have been calculated using the Irving-Kirkwood formalism 54 .
Density profiles across the interface have been calculated following a two-stage process. Firstly, the simulation trajectory is centred about the centre of mass of the primary solvent phase in each simulation (typically water). This reduces artefacts caused by the collective drift of the interface throughout the simulation. The primary phase is calculated by clustering all molecules in the system. The largest cluster is then selected as the primary (liquid) phase. Any molecule further than 0.35 nm from the primary phase is excluded from the centre of mass calculation, and therefore does not affect the centering of the system. A cutoff of 0.35 nm was selected as this corresponds to the first minimum in the radial distribution function (RDF) of water 55,56 . Consequently, water molecules in the vapour phase are excluded in the centre of mass calculation for the bulk liquid water phase, and therefore do not artificially shift the resulting density profiles normal to the interface. After the system has been centred, the density profile is calculated across the interface using a bin size of 0.01 nm. Where applicable, the density profile has been fit to the equation: where ρ l is the liquid density of the primary phase, ρ v is the vapour density of the primary phase, z 0 is the location of the interface, and Δ is the interfacial width. The density profiles of Na + and Cl − ions have subsequently been calculated relative to the definition of the water interface. The density profiles calculated using the above methodology are subject to capillary waves due to thermal fluctuations. Recently, computation techniques have been able to resolve the intrinsic density profile of the liquid-vapour and liquid-liquid interface 21 . That is, the interface between two phases excluding the contribution of thermal fluctuations (capillary waves), which act as to smear the density profiles across the interface. The amplitude of these capillary waves scales as: 2 where the maximum wave vector, q, depends upon the size of the simulated interface. The calculation of the intrinsic density profile has been evaluated by offsetting the amplitude of the thermal fluctuations (ξ) from the interface: where index i sums over all atoms of phase i, and z is the position of the local non-intrinsic interface. The intrinsic density profiles have been calculated using the ITIM method 21 as presented by Sega et al. 57 , using a probe sphere radius of 0.2 Å. Within the ITIM method a probe sphere is moved along test lines perpendicular to the plane of the fluid-vapour or fluid-fluid interface. Once the probe sphere touches the first atom within of the phase of interest, this molecule is marked as being interfacial. This process is repeated over the entire interfacial area in the simulation. The intrinsic density profile is then calculated using the offset (ξ) of the marked interfacial atoms.