Computational predictions of interfacial tension, surface tension, and surfactant adsorption isotherms †

All-atom (AA) molecular dynamics (MD) simulations are employed to predict interfacial tensions (IFT) and surface tensions (ST) of both ionic and non-ionic surfactants. The general AMBER force field (GAFF) and variants are examined in terms of their performance in predicting accurate IFT/ST, g , values for chosen water models, together with the hydration free energy, D G hyd , and density, r , predictions for organic bulk phases. A strong correlation is observed between the quality of r and g predictions. Based on the results, the GAFF-LIPID force field, which provides improved r predictions is selected for simulating surfactant tail groups. Good g predictions are obtained with GAFF/GAFF-LIPID parameters and the TIP3P water model for IFT simulations at a water–triolein interface, and for GAFF/GAFF-LIPID parameters together with the OPC4 water model for ST simulations at a water–vacuum interface. Using a combined molecular dynamics-molecular thermodynamics theory (MD-MTT) framework, a mole fraction of C12E6 molecule of 1.477 (cid:2) 10 (cid:3) 6 (from the experimental critical micelle concentration, CMC) gives a simulated surface excess concentration, G MAX , of 76 C12E6 molecules at a 36 nm 2 water–vacuum surface (3.5 (cid:2) 10 (cid:3) 10 mol cm (cid:3) 2 ), which corresponds to a simulated ST of 35 mN m (cid:3) 1 . The results compare favourably with an experimental G MAX of C12E6 of 3.7 (cid:2) 10 (cid:3) 10 mol cm (cid:3) 2 (80 surfactants for a 36 nm 2 surface) and experimental ST of C12E6 of 32 mN m (cid:3) 1 at the CMC.


Introduction
2][3] Surfactants are employed in a range of applications including washing products, 4,5 pharmaceuticals, 6,7 and petroleum recovery. 8,9The presence of both hydrophilic and hydrophobic parts enables surfactants to lower interfacial tension (IFT), g, between two phases (where one phase is usually water or an aqueous phase). 10,11IFT reduction, (or in the case of a water-air interface -surface tension (ST) reduction), is one of the crucial factors in reflecting the capability of surfactant molecules in real-world applications. 124][15][16][17] Below this concentration, additional surfactant molecules added to the interface help reduce g.Therefore, the value of the CMC and the form of the surface adsorption isotherm, together with the inherent ability of certain surfactants to reduce g more than others, are crucial factors in choosing surfactants for ''real-world'' applications. 18omputer simulation has become a powerful technique for studying and making quantitative predictions of natural phenomena. 19Several previous studies have calculated g for surfactants at both water-oil and water-vacuum interfaces using all-atom (AA) molecular dynamics (MD), [20][21][22][23] and the impact of both ionic and non-ionic surfactants on lowering water-oil IFT or water-vacuum ST has been investigated.5][26][27][28][29][30][31][32][33] While experimental techniques used for g measurements (e.g.pendant drop tensiometry) provide g with changing bulk concentration, C, [34][35][36][37] simulations provide g measurement based on controlling the number of surfactants at the interface, i.e. a controlled surface concentration (G). 18This difference adds difficulty in providing quantitative comparisons between simulation and experimental results.Nevertheless, the comparison can still be made by evaluating adsorption isotherms G(C) experimentally by various methods including colorimetry, chemical and radiochemical techniques, spectroscopic techniques, 38,39 and by employing the Gibbs adsorption isotherm equation 40 where G is surface excess concentration, k B is the Boltzmann constant, T is absolute temperature, g is IFT or ST, C is bulk concentration of surfactants.n = 1 for a non-dissociating adsorbent (e.g.non-ionic surfactant) and n = 2 for a dissociating adsorbent (e.g.ionic surfactant).Recently, Sresht and Jusufi suggested a molecular dynamics simulation-molecular thermodynamics theory (MD-MTT) framework to fully predict the adsorption isotherm graph G(C) using molecular dynamics simulation, and (hence) also provide predictions of G MAX and CMC. 41 key component of all AA MD studies is provided by the molecular force field employed within the simulation.While some literature studies focus on the simulation of g for one category of surfactants at the interface by employing a specific force field (e.g.SDK force field for polyethylene glycol alkyl ethers), [41][42][43][44] developments of force fields for predicting g of various categories of surfactants remain an active research area in high-throughput screening of surfactant formulations, where lower g at CMC (g CMC ) and lower CMC are targeted for higher efficiency and lower cost.45 In the current study, we test the performance of general AMBER force field (GAFF), and variants thereof, in the prediction of g for different surfactants at interfaces.20 We show that the performance of the force field in predicting g at a water-air interface is very sensitive to the water model used, and also to the ability of the force field to predict the correct molecular flexibility and density (r) for hydrocarbon chains.We show that the relatively new OPC4 46 (4-point, rigid optimal point charge) water model works well in combination with the GAFF/GAFF-LIPID in predicting ST for realistic G MAX values, and we show that GAFF/GAFF-LIPID works well in combination with a TIP3P (transferable intermolecular potential 3 points) water model for predicting the IFT for a (typical ''real world'' interface) watertriolein system.

Simulation details
In this study, the GAFF parameter sets 20 and semi-empirical with bond charge correction (AM1-BCC) charge model 47,48 were employed for the AA MD simulations, using the simple harmonic functional form shown in eqn (2), as taken from the Antechamber package 49 obtained from the AmberTools package. 50,51Here, r eq and y eq are equilibration bond length and bond angle, n or n d is multiplicity, o is dihedral angle, g and o d are phase angles, K r , K y , V n and k d are bond, angle and torsional force constants respectively, s ij and e ij are Lennard-Jones parameters, and q i and q j are partial electronic charges.Force fields were generated for the structures shown in Fig. 1.The triolein molecule, which is found in naturally occurring products such as olive oil, was chosen as a typical ''oil molecule''.The sodium dodecyl sulphate (SDS), sodium 4-dodecylbenzene sulfonate (SDBS), sodium 4-(dodecan-2-yl)benzene sulfonate (SD2BS) and sodium laureth-3 sulfate (SLE3S) were chosen as typical ionic surfactants, and polyethylene glycol alkyl ethers (C12E3, C12E5 and C15E7) were chosen as non-ionic surfactants.
The generated GAFF topology and coordinate files were converted to GROMACS format by the ACPYPE script. 52The TIP3P model was used as a compatible water model to be used with GAFF parameter sets. 53We note that GAFF parameters and the TIP3P water model were successfully employed previously by Wilson and co-workers to study self-assembly within several lyotropic systems, e.g.Sunset Yellow (SSY), 54 and a number of ionic cyanine dyes, [55][56][57][58] and this force field combination therefore provides a suitable starting point for this study.In addition, (as described below -see Section 3), alternative 3-site and 4-site water models were also used and the GAFF force was amended/enhanced for some calculations.Finally, the OPLS-AA force field was used for some of the test calculations, to calibrate and understand the influence of different force field factors on IFT.OPLSA-AA force fields were generated using the LigGenPar server, [59][60][61] with the 1.14*CM1A charge model. 60imulations and analysis of results made use of the GRO-MACS 2021.1 package. 62Molecular structures and positions were viewed in Visual Molecular Dynamics (VMD). 63A leap-frog integrator was used for all simulations 64 using a time step of 2 fs, and the Linear Constraints Solver (LINCS) constraint algorithm 65 was applied to all molecular bonds in both equilibration and production runs.Long-range electrostatic interactions were modified by the Particle Mesh Ewald (PME) method, 66 and the long-range dispersion corrections for energy and pressure were applied to long-range van der Waals interactions, having applied a cutoff distance of 1.2 nm for shortrange interactions.

Interfacial tension
For IFT measurements, a simulation box was constructed as follows: water molecules were randomly inserted into a 6 nm Â 6 nm Â 3 nm box; oil molecules were randomly inserted into a 6 nm Â 6 nm Â 6 nm box (the number of molecules was chosen according to their experimental r); and surfactant molecules (number of molecules ranging from 10 to 100) were inserted into a 6 nm Â 6 nm Â 3 nm box such that x, y positions were chosen randomly but they were orientated with their principal axis along the z-axis of the simulation box.These ''sub boxes'' were combined (in a sandwich geometry) along the z direction, with the order of water sub box + oil sub box + water sub box, or water sub box + surfactants sub box + oil sub box + surfactants sub box + water sub box, depending on whether surfactants were required.It should be noted that the head groups of the surfactants were always arranged to point towards the water box, as demonstrated in the initial configuration snapshot of Fig. 2.
After energy minimization, with a steepest descent algorithm, 67 a 200 ps NP N AT pre-equilibration was carried out with a V-rescale thermostat and a Berendsen barostat, 68,69 followed by a 50 ns NP N AT equilibration run with a Nose ][72] A further production run was carried by for 200 ns as a NP N AT ensemble simulation using a Nose ´-Hoover thermostat (with a time constant of 1.0 ps) and a Parrinello-Rahman barostat (with a time constant of 2.0 ps).Here, P N and A are normal pressure and surface area of x-y plane respectively.In both the equilibration and production runs, a temperature of 298.15K and a pressure of 1 bar were maintained.
The NP N AT ensemble runs in GROMACS were carried out using a semi-isotropic pressure coupling scheme, with a reference pressure of 1.0 bar and compressibilities set to 0 bar À1 in the x-y plane and to 4.5 Â 10 À5 bar À1 in z direction respectively.IFT, g, was calculated from the pressure tensor elements (P XX , P YY and P ZZ ) using the equation 73 where L Z is the simulation box length in the z direction, which is normal to x-y plane.The factor of 1 2 arises from the presence of two interfaces in the simulation box and h. ..i t denotes an ensemble average over simulation time.Final values of g were calculated as ensemble averages over 200 ns simulation.Errors for g were obtained from the standard error calculated for 20 Â 10 ns simulation blocks from the 200 ns simulation runs.
We note in passing that to obtain reproducible surface tensions we always use the procedure of producing an initial surfactant layer with head groups pointing towards the water.Following the equilibration steps described above leads to a liquid layer of surfactant with considerable liquid disorder.In principle, this procedure is not mandatory, i.e. we could use a random arrangement of surfactants in a layer or start from surfactants dissolved initially in the aqueous phase.However, although both of these cases tend toward a configuration with the surfactant head groups pointing towards the water, equilibration is very slow (on the time scales used for atomistic simulation).So, from a practical point of view, the procedure outlined above is necessary.

Surface tension
Two 6 nm Â 6 nm Â 6 nm vacuum sub-boxes, two surfactant sub-boxes and a 6 nm Â 6 nm Â 6 nm water sub-box were combined in z direction (with the order of vacuum + surfactants + water + surfactants + vacuum box) to form the initial simulation box.After minimisation, a 200 ps NVT pre-equilibration was run with the V-rescale thermostat, followed by a 50 ns NVT equilibration with the Nose ´-Hoover thermostat.The production run was carried out by a 200 ns NVT ensemble simulation with the Nose ´-Hoover thermostat (with a time constant of 1.0 ps).In both the equilibration and production runs, a temperature of 298.15K was maintained.Data analysis for ST simulations employed the same procedure described in Section 2.2.

Density of bulk phases
We calculate the bulk density for several systems: methanol, 1,2-dimethoxyethane (DME), dodecane, pentadecane and C12E5.For these simulations, molecules were randomly inserted into a 7 nm Â 7 nm Â 7 nm box according to its experimental r. 74Long-range electrostatic interactions were modified by the Particle Mesh Ewald (PME) method, and the long-range dispersion corrections for energy and pressure were This journal is © the Owner Societies 2024 applied to long-range van der Waals interactions, both with a cut-off distance of 1.2 nm.After energy minimisation, a 100 ps NVT equilibration was run with the V-rescale thermostat, followed by a 1 ns NPT equilibration with the V-rescale thermostat and Berendsen barostat.The MD production run was carried out using a 10 ns NPT ensemble simulation with a Nose ´-Hoover thermostat (with a time constant of 1.0 ps) and a Parrinello-Rahman barostat (with a time constant of 2.0 ps).In both equilibration and production runs, a temperature of 298.15K and a pressure of 1 bar were maintained.

Free energy of hydration
The target molecule was inserted into a rhombic dodecahedron box with a distance between the solute and the box edge of 1.2 nm, followed by a solvation step using explicit water molecules in the box.After energy minimisation and equilibration, the production run was carried by a 5 ns NPT ensemble simulation at each of 14 l states.The l values were 0.0, 0.1, 0.2, 0.3, 0.4, 0.45, 0.5, 0.55, 0.6, 0.65, 0.7, 0.8, 0.9 and 1.0 (l = 0 for the decoupled state and l = 1 for the fully coupled state).Here, both van der Waals and electrostatic interactions were turned off when l = 0 (soft core interactions were applied to avoid singularities, where soft core-a = 1.0, soft-core-s = 0.3 and soft core power = 1), 75 and all interactions were turned on when l = 1.At each stage, both van der Waals and electrostatic interactions were gradually turned on together.In both the equilibration and production runs, a temperature of 298.15K and a pressure of 1 bar were maintained, if applicable.In addition, an SD-integrator (a leap-frog stochastic dynamics integrator) was used with a time step of 2 fs in MD simulations, and the LINCS constraint algorithm was applied to all molecular bonds in both the equilibration and production run.The free energy of each l state was calculated using the Bennett acceptance ratio (BAR) method, 76 and the final free energy of hydration (DG hyd ) was calculated by the difference of free energy obtained between the states l = 0 and l = 1 states.

Free energy of adsorption
A 6 nm Â 6 nm Â 6 nm vacuum sub box and a 6 nm Â 6 nm Â 6 nm water sub box (with a surfactant molecule) were combined along the z-direction to form the initial simulation box.After energy minimisation, a 200 ps NVT pre-equilibration was performed with the V-rescale thermostat, followed by a 50 ns NVT equilibration with the Nose ´-Hoover thermostat.In equilibration, a temperature of 298.15K was maintained.
To obtain the free energy of adsorption (DG ads ), umbrella sampling was employed to calculate the potential of mean force (U PMF ).U PMF is compiled from average forces calculated as a function of the intermolecular coordinate between two interacting species using where R A and R B are the initial and final intermolecular coordinates, f (R) is the force in terms of the intermolecular coordinate and h. ..i t donates an ensemble average over time for individual simulations. 77,78e employed centre-of-mass (COM) pulling, 79 between the centre of mass of a surfactant molecule and the centre of mass of a water slab.COM windows were selected with a neighbouring distance of 0.15 nm between them (25 windows in total).Each window was equilibrated for 1 ns followed by a 10 ns production run to obtain a range of configurations in which the COM of the pulling species was constrained.
To plot U PMF as a function of full and continuous intermolecular coordinates, the energy values in adjacent windows were reassembled to produce a continuous function using a weighted histogram analysis method (WHAM). 79The value of DG ads was calculated by the difference between the highest and lowest values in the U PMF curve.The statistical error raised from the WHAM, due to limited counting (i.e.finite time steps in simulations), was conveniently calculated by bootstrap analysis. 80Results and discussion

Initial testing and validation of force fields of surfactants at water-oil interfaces
In initial tests, the IFT of a water-triolein interface shows a simulated value of 34.03 AE 0.93 mN m À1 using GAFF parameters and the TIP3P water model, which compares well with the experimental value of 32 mN m À1 . 88his interface was then used in initial test simulations to evaluate the IFT of the surfactants (i.e.SDS, SDBS, SD2BS, SLE3S, C12E3, C12E5 and C15E7) as a function of the number of surfactant molecules at the water-triolein interface within the symmetrical slab geometry shown in Fig. 2. The results are presented in Fig. 3, where initially both triolein and surfactant molecules are simulated with GAFF parameters.
We note that at high G the interface starts to distort from the planar geometry as it becomes oversaturated with surfactant.This leads to an ''apparent rise'' in IFT, as the increase in the area of the interface is not being accounted for.Hence, to make a comparison with experimental results, we, where available, compare the calculated IFT with the experimental value of IFT measured at G MAX .
Unexpectedly high IFT values were obtained from simulations at high G for both ionic and non-ionic surfactants at the water-triolein interface in Fig. 3.For example, the simulated IFTs of 70 SDS and 80 C12E5 molecules at a 36 nm 2 crosssection water-triolein interface (approximate values expected at G MAX ) show values of 21.12 AE 0.92 mN m À1 and 17.26 AE 1.27 mN m À1 respectively, whereas the experimental IFT CMC of ionic and non-ionic surfactants at a water-oil interface can approach 5 mN m À1 in the case where the oils are triglycerides with similar structure to triolein 34,82,83,[89][90][91][92][93][94] (noting that the experimental values of G MAX of SDS and C12E5 are 3.2 Â

Improving predictions of IFT in GAFF
Here, C12E5 (Fig. 1), a non-ionic alcohol ethoxylate was chosen to perform the initial tests and improvements of the GAFF parameter set in combination with the TIP3P water model.The three key parts of alcohol ethoxylates consist of a hydrophobic hydrocarbon tail group, together with a head group composed of a terminal hydrophilic alcohol and an extended ethylene oxide (EO) chain.7][98] Accordingly, we consider the performance of GAFF for representative fragment molecules: methanol, 1,2-dimethoxyethane (DME) and dodecane. 42,86,99,100t is known already that GAFF underestimates DG hyd of alcohol and DME, with a simulated DG hyd for methanol of À3.49AE 0.02 kcal mol À1 (experimental data of À5.10 AE 0.60 kcal mol À1 ), and a simulated DG hyd for DME of À3.10 AE 0.03 kcal mol À1 (experimental data of À4.84 AE 0.60 kcal mol À1 ). 75,86lso, GAFF is known to provide relatively poor predictions of r for long hydrocarbon chains.Here, long hydrocarbon chains are slightly too stiff in GAFF and, consequently, r values are too high with GAFF parameters from bulk phase simulations. 100he r of dodecane molecules simulated with GAFF shows a value of 824 kg m À3 , 98 compared to the experimental value of 750 kg m À3 . 87ere, it is interesting to compare the results from tests of the OPLS-AA force field.According to Akinshina et al., 99 OPLS-AA parameters successfully reproduced the hydration of the short EO chains of the chromonic amphiphile 2,3,6,7,10,11hexa-(1,4,7-trioxa-octyl)-triphenylene (TP6EO2M) in aqueous solution, and hence the chromonic aggregation behaviour of TP6EO2M molecules in bulk water (i.e.single molecule-cross section stacks).In contrast, GAFF parameters underestimate the hydration of EO chains, resulting in a high aggregation free energy for TP6EO2M molecules and compact disordered clusters.Calculation of DG hyd of DME with GAFF and OPLS-AA parameters using the methodology of Section 2.5 shows that the OPLS-AA parameters estimate enhanced hydration of DME with a value of À23.31 kJ mol À1 , compared to GAFF parameters with a value of À15.97 kJ mol À1 (experimental value À20.25 kJ mol À1 ). 75FT simulations were performed (with GAFF and OPLS-AA parameters) for a water-octane interface system and a wateroctane interface system with 80 C12E5 molecules at each of the two 36 nm 2 interfaces (Table 1), 35 in order to test the relationship between DG hyd and calculated IFTs.OPLS-AA parameters estimate a much lower IFT for the 36 nm 2 water-80 C12E5octane system (with both TIP3P and SPC/E water models) compared to the GAFF parameters, noting that the OPLS-AA and GAFF parameters simulate similar IFTs arising from the pure water-octane interface.It could reasonably be hypothesised that the force field, which simulates a more negative DG hyd for the head groups of surfactants, would estimate a lower IFT for those surfactants at a water-oil interface.101 However, OPLS-AA parameters are not capable of reducing the IFT of the water-70 SDS-octane system to a much lower value than GAFF parameters: a value of 25.00 AE 0.54 mN m À1 is obtained from GAFF and a value of 25.71 AE 0.47 mN m À1 is obtained for OPLS-AA with the SPC/E water model.
Although DG hyd might be one important factor influencing the accurate determination of IFTs in the MD, further validations are still needed.GAFF and OPLS-AA have different functional forms and parameter sets. 20,22The difference between the two force fields might not provide a clear conclusion about the relationship between DG hyd and IFT.However, it is possible to investigate the relationship between DG hyd and IFT by careful modification of the GAFF force field.
Here, three separate modifications of GAFF parameters are performed (Table 2).First, following the approach of Fennell et al., the partial atomic charge of the EO chains in GAFF is scaled by a factor of 1.1 to correctly calculate DG hyd of DME. 86econd, following Boyd, 74 the e of ''atomtype os'' (ether oxygen of EO chains) in GAFF is changed from 0.71128 kJ mol À1 to 1.60656 kJ mol À1 , which improves the calculations of DG hyd for DME.Third, GAFF-DC parameters are applied to the alcohol functional group to correctly simulate DG hyd of methanol. 86t was initially hypothesised that a modified GAFF, which simulates correctly DG hyd of either DME or methanol, would estimate a lower IFT for the water-C12E5-octane interface. 44,101owever, both the second and third modifications of GAFF parameters surprisingly lead to the prediction of higher IFTs (Table 2).As a result, we additionally determined the influence of force field modifications on both r and DG hyd for methanol, DME and dodecane (Table 3).OPLS-AA parameters simulate enhanced hydration for DME, compared to GAFF parameters, although both OPLS-AA and GAFF underestimate DG hyd of methanol (i.e.alcohol head group of non-ionic surfactant).However, the results indicate that the lower IFT value of the 36 nm 2 water-80 C12E5-octane system obtained from OPLS-AA is probably correlated with the lower r of methanol and DME predicted with OPLS-AA.GAFF-DC, which simulates a higher r for methanol, produces a higher IFT for the 36 nm 2 water-80 C12E5-octane compared to GAFF parameters.In addition, it should be noted that both GAFF and OPLS-AA simulate a r for dodecane that is too high compared to the experimental values, as has been noted elsewhere. 100,102

IFT predictions using a combination of GAFF and GAFF-LIPID force fields
From the above, correctly capturing r for pure components plays a key role in predicting IFT simulations successfully.4][105] Thus, GAFF-LIPID, 100 which makes better predictions for r, the heat of vaporisation and the isothermal compressibility of long carbon chain compounds, can be suggested for the tail group parameters of surfactants in the IFT simulation.In addition, Lennard-Jones parameters (s and e) of ''atomtype c3'' described in GAFF-LIPID are utilised in the PEO chains of non-ionic surfactants (C12E3, C12E5 and C15E7) to simulate a reduced r for the heads.However, the use of GAFF-LIPID ''atomtype c3'' for SLE3S PEO chains does not result in a significant difference compared to the original GAFF, which could be due to the existence of the sulfate group.Therefore, we chose the original GAFF parameters to simulate PEO chains for SLE3S (see ESI, † for details).GAFF-LIPID parameters are also tested to describe the acyl chains in the triolein molecule.For the phenyl and sulfonate functional groups involved in the SDBS and SD2BS molecules, GAFF parameters are capable of simulating the correct intercolumn distances and free energies of binding of a molecule to a stack of Sunset Yellow 54 and a thiacyanine dye, 55 a system consisting of aromatic rings and sulfonate functional groups.Thus, the original GAFF parameters are applied directly to ionic surfactant head groups.As discussed below, this combination of GAFF/GAFF-LIPID parameters provides improved IFT predictions for surfactants at the water-triolein interface (Fig. 4).
Because of the lack of experimental data for some surfactants at the water-triolein interface, a water-surfactant-octane system is also employed for validation purposes.IFT simulation a IFT CMC of water-C12E4-hexadecane interface at a G MAX of 3.6 Â 10 À10 mol cm À2 . 82C12E4 surfactant shows a same experimental G MAX as C12E5, both with a value of 3.6 Â 10 À10 mol cm À2 (i.e.80 surfactants at a 36 nm 2 interface). 82In addition, the water-hexadecane interface shows a similar experimental IFT to the water-octane interface, where the IFT of the water-hexadecane interface is 55.2 mN m À1 and the IFT of the water-octane interface is 52.5 mN m À1 . 81 ). 82IFT calculations for the 36 nm 2 water-70 SDS-octane system shows a value of 18.45 AE 0.48 mN m À1 , when GAFF-LIPID is used for SDS tail groups with the SPC/E water model, compared to a value of 25.00 AE 0.54 mN m À1 simulated from GAFF and SPC/E (experimental value of 9.7 mN m À1 ). 83 should be noted that the SPC/E water model with GAFF parameters provides a more accurate IFT for the water-octane interface than the TIP3P water model (Table 1).In contrast (see Section 3.4), the TIP3P water model with GAFF/GAFF LIPID parameters is best for the water-triolein interface (Table 4).The reason why different oil interfaces require different water models for optimal IFT results is yet to be fully understood.However, it has been previously noted that subtle changes to water models can have a significant effect in both structural and dynamic properties. 106

IFT predictions with GAFF/GAFF-LIPID for different water models
Different water models (SPC/E, TIP4P, TIP4P/2005, OPC4 in addition to TIP3P) are tested with GAFF/GAFF-LIPID parameters at the water-triolein interface (Table 4).The results suggest that the TIP3P water model gives a better performance in combination with this force field compared to the other tested water models.The combination of GAFF/GAFF-LIPID and TIP3P provides a value of 31.96AE 0.92 mN m À1 (experimental value of 32 mN m À1 ). 88Noting that this is in line with a good performance from GAFF-derived parameters with the TIP3P water model in simulating IFT for the water-toluene interface. 107alculations for pentadecane in combination with different water models again show that IFT and DG hyd predictions are not strongly correlated (Table 5).Here, the IFT values for the water-pentadecane interface decrease from 71.07 AE 0.83 mN m À1 (GAFF and SPC/E) to 55.05 AE 0.44 mN m À1 (GAFF-LIPID and SPC/E) using the more flexible GAFF-LIPID chains (experimental value of 54.9 mN m À1 ); 81 noting that for the same systems the values of DG hyd increase from 22.48 AE 0.58 kJ mol À1 to 30.13 AE 0.78 kJ mol À1 (experimental value 17.30 kJ mol À1 ). 75GAFF-LIPID also simulates reasonably good values for IFT of a 36 nm 2 water-80 C12E5-octane interface and r for the surfactant molecule C12E5 (Table 6).Here, however   This journal is © the Owner Societies 2024 DG hyd for C12E5 is underestimated.][110] Finally, it is worth noting that, through visualization of dense surfactant layers at the interface, it is evident that GAFF parameters give significantly stiffer surfactant tails in comparison to GAFF-LIPID, resulting in lower areas per surfactant and higher order parameters in comparison to experiments; 100 and, consequently, higher simulated IFT values are obtained for the same G.In these cases, using GAFF chains, more surfactant molecules can ''fit'' at a crowded interface because of closer molecular packing.For example, in simulations of watersurfactants-triolein interfaces, ''frozen islands'' of surfactants and holes in the surfactant layer interface have been observed with GAFF parameters, whereas a liquid state of surfactants at the interface is always observed with GAFF-LIPID chains (Fig. 6).

Surfactants at a water-vacuum interface
Despite the excellent performance of TIP3P for water-triolein IFTs, TIP3P water is found to provide poor predictions for ST at a water-vacuum interface.We see this when a variety of different water models are tested (Table 7), with several models failing to correctly predict the ST of pure water or the reduction in ST arising from surfactants.However, the 4-point rigid optimal point charge (OPC4) water model, which is designed to provide a better description of the electrostatic potential around a water molecule, 46 is found to dramatically improve ST predictions using GAFF-derived force fields (Fig. 5).As noted above, the reason why different interfaces require different water models for optimal IFT results or ST results is yet to be fully understood.We have also included the ST simulated from C12E6 (Fig. 1) for prediction of the adsorption isotherm in Section 3.7.ST results for the water-vacuum interface with the OPC4 water model show a value of 70.66 AE 0.20 mN m À1 (experimental value of 72 mN m À1 ), 111 and for a 36 nm 2 water-100 SDS-vacuum shows a value of 34.51 AE 0.43 mN m À1 (experimental value of 39.5 mN m À1 ). 112   A further comparison is carried out between the simulated and experimental STs of the non-ionic surfactant C12E3 at a water-vacuum interface using the promising OPC4 water model.The experimental ST CMC is approximately 30 mN m À1 and G MAX is 4.5 Â 10 À10 mol cm À2 , 113 compared to a simulated ST of 33.32 AE 0.37 mN m À1 at this G MAX , which corresponds to approximately 100 molecules of C12E3 at a 36 nm 2 interface.

Maximum packing of surfactants
Good predictions of G MAX from simulation would ideally play an important role in screening new surfactants before synthesis.However, it should be noticed that within the framework of atomistic simulations, it is relatively easy for surfactants to become overpacked at the interface and the surface layer in the simulation box to become metastable, compared to multibody dissipative particle dynamics (M-DPD) simulations. 114Surfactants prefer to stay at an interface in AA MD even when the number of molecules leads to G MAX being exceeded.The experimental G MAX will occur at the CMC, when the chemical potential of the monomers at the interface is the same as that in solution and micelles.However, we have no good way to determine this equilibrium within AA simulations.However, we expect AA simulations to offer some information about G MAX , arising from packing considerations at the interface.
Head group (i.e.-SO 4 À ) density graphs of four surface concentrations of SDS molecules (i.e.50 SDS, 70 SDS, 80 SDS and 130 SDS) at a 36 nm 2 water-triolein interface are presented in Fig. 6.For the maximum packing of 70 SDS molecules, consistent heights are observed for surfactant head groups at the two interfaces.However, above this packing, variations from these smooth distributions start to be seen, eventually leading to major distortions of the surfactant layer.This suggests that it might be reasonable to estimate G MAX of surfactants from qualitative observation of a flat monolayer of surfactants.For example, it is straightforward to visualise that the monolayer of 70 SDS molecules at a 36 nm 2 interface is relatively flat compared to 80 SDS molecules.However, this visualisation can be difficult to achieve in practice for many systems.We note that for 130 molecules of SDS (bottom frame in Fig. 6), which is well above G MAX , major distortions of the interface occur, with the system essentially creating additional interface to adjust to a large number of surfactants.In this case, the calculation of g via eqn (3) immediately fails because additional interfacial area has been created.
Interestingly, by applying qualitative observation of the monolayer of surfactants, the simulated G MAX and IFT at this G MAX of SDBS at the water-triolein interface are 100 molecules at a 36 nm 2 interface and 10.06 AE 1.45 mN m À1 respectively, whereas the simulated G MAX and IFT at this G MAX of SD2BS are 80 molecules at a 36 nm 2 interface and 3.69 AE 1.61 mN m À1 respectively.The difference between these two molecules might be explained by the lower density of dodecan-2-yl chains at the interface, compared to dodecyl. 115,116This observation is consistent with experimental value where experimental G MAX and ST CMC of SD2BS at water-air interface (''hard river'' water, 303.15 K) are 4.16 Â 10 À10 mol cm À2 and 36.4 mN m À1 respectively, experimental G MAX and ST CMC of SD4BS (''hard river'' water, 303.15 K) are 3.44 Â 10 À10 mol cm À2 and 28.2 mN m À1 respectively and experimental G MAX and ST CMC of SD6BS (''hard river'' water, 303.15 K) are 3.15 Â 10 À10 mol cm À2 and 27.5 mN m À1 respectively. 117

Adsorption isotherms of surfactants
Using an MD-MTT framework (a combination of MD and molecular thermodynamics theory), 41,118 it is, in principle, feasible to fully determine adsorption isotherms for non-ionic surfactants at a water-vacuum interface by simulation.AA MD determines three parameters: (i) the radius r of the surfactant, (ii) the second virial coefficient B for the interaction of surfactants at the interface and (iii) the free energy of adsorption DG ads of the surfactant.By fitting the ST graph of C12E6 molecules at a water-vacuum interface, parameters r and B can be deduced.The relationship between surface pressure and surface concentration of non-ionic surfactants follows the equation 41,118 where P is the surface pressure, G is the surface concentration of surfactants, Z (= Ga) is the packing fraction of surfactants at the interface, a (= pr 2 ) is the excluded area of a surfactant and B is the pairwise coefficient for short-range van der Waals interaction between surfactants. 41After the determination of the parameters r and B, the relationship between the mole fraction of surfactants in bulk water and the surface where X is the mole fraction.
Here, GAFF/GAFF-LIPID parameters with the OPC4 water model, and the OPLS-AA parameters with the SPC/E water model (as used by Sresht et al.) 41 are used to simulate the ST of the C12E6 molecules at the water-vacuum interface.The fitted data from the GAFF/GAFF-LIPID and OPC4 model provides a r value of 0.279 nm and B value of 0 nm 2 , and the OPLS-AA and SPC/E model provides a r value of 0.272 nm and B value of À0.364 nm 2 (Fig. 7).We note that neither graph in Fig. 7 exactly followed the theoretical functional form suggested.However, the start and end points at G = 0 nm À2 and G = G MAX are quite close and the fitted curve for the OPLS-AA and SPC/E model quite closely follows the simulation data.It should be noticed that the fitting data assumes a prior knowledge of an experimental G MAX of surfactants (i.e.2.22 nm À2 for C12E6), 119 and simulated IFT values that exceed G MAX are discarded.To fully predict the parameters r and B from the simulations, qualitative observations of a surfactant monolayer (Fig. 6) could be used to determine an approximate G MAX of C12E6 surfactants from simulations.Then, a calculated G MAX could be determined from the adsorption isotherm (Fig. 10) and compared with the previous G MAX value.If they are found to be different, then the fitting procedure in Fig. 7 could be iterated to calculate a new G MAX .
As part of an MD-MTT framework, it is necessary to know the change in the free energy for a molecule that is transferred from bulk water to the water-vacuum interface.The free energy curves for GAFF/GAFF-LIPID and OPC4, and GAFF/GAFF-LIPID and SPC/E are shown in Fig. 8 in the form of potentials of mean force, and the calculated DG ads for different force fields and water model combinations are shown in Table 8.Interestingly, the results in Table 8 suggest that GAFF/GAFF-LIPID parameters (for surfactants) with the SPC/E water model reproduce the experimental DG ads most closely.This is even though the combination of GAFF-LIPID and the SPC/E water produces a DG hyd for pentadecane that is more positive than the experiment, worse than GAFF and SPC/E, and worse than the excellent prediction for GAFF and TIP3P (Table 5).Here, of course, DG hyd is not the only contribution to DG ads , and most of the experimental DG ads arises from a lowering of the ST. 40,95This explains why GAFF/GAFF-LIPID parameters, which we see can provide accurate predictions of the lowering of ST, provide good agreement with the experimental DG ads .
Based on simulation results from Fig. 7 and Table 8, r, B and DG ads of C12E6 molecule (Table 9) are obtained with GAFF/ GAFF-LIPID parameters and the OPC4 water model (for r and B)  and the SPC/E water model (for DG ads ).As discussed in Section 3.3, the reasons why different physical measurements perform better with different optional water models are yet to be fully explored.Though we know that even small changes in the water models (charge and shape etc.) can lead to measurable differences in radial distribution functions and this is likely to have an effect on surface properties. 106rom the calculated adsorption isotherm graph (Fig. 9), the experimental mole fraction of C12E6 of 1.477 Â 10 À6 gives a simulated G MAX of 76 C12E6 molecules at a 36 nm 2 watervacuum interface, which corresponds to a simulated ST of 35 mN m À1 (Fig. 10).9][120] The agreement between the simulation data and the experimental data also validates the use of GAFF-LIPID ''atomtype c3'' for non-ionic surfactant PEO chains.
With further validation, the MD-MTT framework could be extended to predict both ionic surfactants and non-ionic surfactants at the water-oil and water-vacuum interface, with the help of the force fields refined in this paper. 121,122In combination with CMCs calculated directly from dissipative particle dynamics (DPD) simulations, [123][124][125][126][127] or slightly more expensive approaches such as many-body dissipative particle dynamics (M-DPD) 114 or coarse-grained molecular dynamics, 128 this provides a way to fully predict G MAX , g CMC , CMC and the adsorption isotherm of surfactants directly from simulation.
Based on the simulations presented in this study, we summarise the most suitable force field combinations in Table 10 for each AA MD simulation system mentioned above.

Conclusions
AA MD simulations have been used to predict accurately the value of g and the reduction of g due to surfactants at the watertriolein and water-vacuum interfaces.In this study, we tested the performance of GAFF, and variants thereof, in the    This journal is © the Owner Societies 2024 prediction of g for different surfactants at interfaces. 20Correct predictions of r for pure components are important to ensure accurate prediction of g.Thus, GAFF-LIPID, 100 which makes better predictions for r and isothermal compressibility of long carbon chain compounds, provides much better predictions for g, in most systems, than GAFF.Based on the results from this study we highlight recommended force field/water model combinations for predictions of IFT, ST, DG hyd and DG ads from AA MD simulation systems as summarised in Table 10.The TIP3P water model was used to provide good IFT predictions for water with a series of oils and was found to be compatible with GAFF and GAFF-derived parameter sets. 53However, TIP3P water provides poor predictions for ST of water-vacuum interfaces.A variety of different water models were tested to improve these predictions, with several models failing to correctly predict the ST of pure water or the reduction in ST arising from different surfactants.However, the 4-point rigid optimal point charge (OPC4) water model, which provides a better description of the electrostatic potential around a water molecule, is found to dramatically improve ST predictions using GAFF-derived force fields. 46 major issue related to AA MD predictions of g values arises from difficulties in predicting the maximum surface concentration, G MAX , to find the maximum reduction in ST for a given surfactant.With the help of a MD-MTT framework and an experimental CMC, simulations successfully predict a ST CMC for C12E6 of 35 mN m À1 at G MAX with GAFF/GAFF-LIPID parameters, which compares favourably to the experimental value of 32 mN m À1 .[118][119][120] Finally, it is worth highlighting some of the slightly unsatisfactory findings arising from this work.TIP3P, the best water model for the water-triolein interface (with the GAFF/GAFF LIPID force field), was much poorer than OPC4 in predicting surface tensions at the water-vacuum interface.Moreover, the various water models have differing performance in relation to the predictions of DG ads and DG hyd .Here, it is worth recalling that all these water models are effective two-body potentials, in which higher body terms have been combined into the effective pair potential.Although this is reasonable as an approximation for bulk water, the interactions vary at an interface or when the dielectric environment changes.Hence, it is difficult to predict which water model will behave best in combination with different oil and surfactant force fields at different interfaces.An interesting extension to the current work may arise from the use of polarisable force fields, 129 where the ability of the force field to adjust to changes in the molecular environment and relative permittivity potentially provide a way to address these issues.

Fig. 3
Fig. 3 Simulated interfacial tensions (IFT) for different surfactants at a water-triolein interface, with GAFF parameters and the TIP3P water model.Lines connecting points act as a guide for the eye, noting that at very high values of surface packing the interface becomes oversaturated with surfactant, leading to a nonphysical increase in IFT.

of a 36
nm 2 water-80 C12E5-octane interface shows a value of 8.83 AE 0.35 mN m À1 , with GAFF-LIPID for both C12E5 tail group and head group with the SPC/E water model, compared to a value of 25.21 AE 0.80 mN m À1 simulated from GAFF and SPC/E (experimental value of water-C12E4-hexadecane interface of 3.1 mN m À1

Fig. 5
Fig. 5 Simulated surface tensions (ST) of different surfactants at a watervacuum interface, with GAFF/GAFF-LIPID parameters and the OPC4 water model.In the top panel, the blue vertical line represents the experimental G MAX for SDS.In the lower panel the blue vertical line represents the experimental G MAX for C12E3, and the red-line represents the experimental G MAX for C15E7.The latter is coincident with the value of G MAX for C12E5 and C12E6.

Fig. 7
Fig. 7 Simulated and fitted surface tensions (ST) data of different number of C12E6 molecules at water-vacuum interface, with GAFF/GAFF-LIPID parameters and the OPC4 water model, which gives r value of 0.279 nm and B value of 0 nm 2 (top) and OPLS-AA parameters and the SPC/E water model, which gives r value of 0.272 nm and B value of À0.364 nm 2 (bottom).

Fig. 8
Fig. 8 Potential of mean force curve for pulling a C12E6 molecule from the water-vacuum interface to bulk water, with GAFF/GAFF-LIPID parameters and the OPC4 water model (top), and GAFF/GAFF-LIPID parameters and the SPC/E water model (bottom).

Fig. 9
Fig. 9 Calculated adsorption isotherm of C12E6 molecules.r and B are obtained from simulations with GAFF/GAFF-LIPID parameters and the OPC4 water model, and DG ads is obtained from the simulation with GAFF/GAFF-LIPID parameters and the SPC/E water model.

Fig. 10
Fig. 10 Simulated and experimental surface tensions (ST) against mole fraction of C12E6 molecules.r and B are obtained from simulations with GAFF/GAFF-LIPID parameters and the OPC4 water model, and DG ads is obtained from the simulation with GAFF/GAFF-LIPID parameters and the SPC/E water model.GAFF represents the original GAFF parameters for C12E6, and GAFF/GAFF-LIPID represents GAFF-LIPID parameters for both tail groups and head groups of C12E6.

Table 3
Densities (r) and free energies of hydration (DG hyd ) of methanol, DME and dodecane, with GAFF, OPLS-AA and GAFF-DC parameters using the TIP3P water model À3 DG hyd /kJ mol À1 r/kg m À3 DG hyd /kJ mol À1 r/kg m À3 DG hyd /kJ mol À1

Table 4
Interfacial tensions (IFT) of the 36 nm 2 water-triolein interface, with GAFF/GAFF-LIPID parameters for triolein using different water models

Table 7
Surface tensions (ST) of a 36 nm 2 water-vacuum interface (second column) and a 36 nm 2 water-100 SDS-vacuum interface (G of 2.78 nm À2 ) (third column), with GAFF/GAFF-LIPID parameters (for SDS) using different water models

Table 9
Radius (r), second virial coefficient (B) and free energy of adsorption (DG ads ) of C12E6 molecule.r and B are obtained from simulations with GAFF/GAFF-LIPID parameters and the OPC4 water model, and DG ads is obtained from the simulation with GAFF/GAFF-LIPID parameters and the SPC/E water model r/nm B/nm 2 DG ads /kJ mol À1