Molecular Simulations of Surfactant Adsorption on Iron Oxide from Hydrocarbon Solvents

The performance of lubricant additives, such as organic friction modifiers (OFMs), depends critically on their ability to adsorb onto the surfaces of moving components and form protective self-assembled layers (SAMs). Therefore, understanding the relationship between the concentration of the additive in the base oil and the resulting surface coverage is extremely important for lubricant formulations, as well as many other surfactant applications. Here, we use molecular dynamics (MD) simulations to study the adsorption isotherms of three different OFMs, stearic acid (SA), glycerol monoostearate (GMS), and glycerol monooleate (GMO), onto a hematite surface from hydrocarbon solvents, n-hexadecane and poly-α-olefin (PAO). First, we calculate the


Introduction
The adsorption and self-assembly of surfactants at solid-liquid interfaces plays a key role in a wide range of technological processes, from nanoparticle stabilisation to ore flotation. 1 In tribology, the adsorption of polar lubricant additives from nonpolar base oils onto solid surfaces is critical to their effective operation. 2 For example, organic friction modifiers (OFMs), such as stearic acid (SA), adsorb and form thin (∼2 nm) self-assembled monolayers (SAMs) on steel surfaces, which significantly reduce friction and wear. 3 Similarly, inorganic friction modifiers, such as molybdenum dialkyldithiocarbamate (MoDTC), adsorb and decompose to form thin (∼1 nm), low-friction lamellar molybdenum disulphide films on steel. 4 Moreover, antiwear additives, such as zinc dialkyldithiophosphate (ZDDP), adsorb, form SAMs, 5 and decompose on steel surfaces to form thick (∼100 nm) protective tribofilms, composed mostly of metal phosphates. 6 Modern lubricant formulations contain at least ten different types of additive dissolved in a blend of base oils. 7 Combinations of different additives often behave antagonistically, that is the additive blend gives poorer performance than any of the individual additives alone at equal concentration. 8 Prominent examples of antagonistic lubricant additive combinations include overbased calcium sulphonate detergents with ZDDP, 9 succinimide dispersants with ZDDP, 10 and MoDTC with ester OFMs. 11 Many of these antagonistic interactions are due to competitive adsorption processes, whereby the most surface-active component occupies the majority of surface sites and inhibits the action of the other additives. 8 Therefore, it is important to understand the relative adsorption strength of the various lubricant additives from base oils on steel surfaces.
OFMs are one of the earliest and most widely studied types of lubricant additive. 12 Hardy and Doubleday 13 were the first to suggest that OFM friction reduction was due to the formation of vertically-oriented surfactant monolayers on the sliding surfaces. They found that fatty acids with longer chain length (C 4 -C 12 ) showed lower boundary friction and that fatty acids gave lower friction than fatty alcohols 13 . Bowden and Leben 14 showed that the friction of SA (C 18 ) monolayers deposited on steel surfaces with the Langmuir-Blodgett technique was initially the same as that for SAMs formed from a hydrocarbon base oil. Fowkes 15 showed that OFM adsorption strength from hydrocarbon solution onto steel varied in the order carboxylic acid > amine > alcohol, and Groszek 16 showed that this correlated with antiwear performance. Rowe 17 also found that adsorption strength and antiwear performance were correlated for OFMs, but with a slightly different order (amine > carboxylic acid > alcohol). Using this data, Rowe developed a mathematical model to calculate thermodynamic quantities describing adsorption from experimentally-obtained wear rates. 17 Spikes and Cameron 18 showed that adsorption and lubricant failure both had a similar dependence on fatty amine concentration and temperature, suggesting that lubricant failure could be due to desorption of the OFM molecules. Jahanmir and Beltzer 19 studied the adsorption and friction behaviour of a wide range of carboxylic acid, ester, alcohol, amine, and amide OFMs with different alkyl tailgroup lengths, branching, and unsaturation. They found that additives that adsorbed more strongly generally gave lower friction. 19 Jahanmir and Beltzer 20 also built upon the model due to Rowe 17 to describe the relationship between the friction coefficient and adsorption energy for a two-component lubricant (i.e. base oil and one additive), assuming that the minimum friction coefficient observed at high OFM concentration was due to the formation of a close-packed SAM. 20 Recent depletion isotherm and polarised neutron reflectometry (PNR) experiments by Wood et al. 21 for the adsorption from n-dodecane onto iron oxide surfaces have shown that OFMs with linear, saturated tailgroups (SA) do form monolayers with a high surface coverage, similar to that expected for a close-packed monolayer, at high concentration. Conversely, OFMs with kinked, Z-unsaturated tailgroups, such as oleic acid (OA), form monolayers with much lower surface coverage. This observation explained the higher friction coefficient for OA than SA reported in nanoscale atomic force microscopy (AFM) by Lundgren et al. 22 and microscale surface forces apparatus (SFA) experiments by Ruths et al. 23 The same observations were made at the macroscale by Campen et al. 24 for OA and glycerol mono-oleate (GMO) compared to SA and glycerol mono-stearate (GMS) using steel substrates and an n-hexadecane solvent. Using a range of OFMs with different headgroups and tailgroups, Fry et al. 25 showed through a combination of quartz crystal microbalance (QCM), ellipsometry, and tribology experiments that OFMs which formed higher coverage monolayers gave lower initial friction coefficients.
In recent years, molecular simulations have provided additional insights into the nanoscale behaviour of a wide range of lubricant additives, including OFMs. 26,27 The adsorption energy of lubricant additives on metal oxide surfaces can be accurately determined with first principles methods such as density functional theory (DFT). [28][29][30] However, due to its very high computational cost, DFT is typically performed in vacuum at zero temperature, and can only capture the enthalpy contributions to the adsorption free energy. For many surfactants at solid-liquid interfaces, entropy plays a central role in the adsorption process. 31 These entropic contributions can be calculated using classical molecular dynamics (MD) simulations, 32 which can be used to access much larger time and length scales than DFT.
Nevertheless, capturing the adsorption and self-assembly processes from solution using realistic surfactant concentrations represents a considerable computational challenge. 27 Until recently, a certain surface coverage of surfactant molecules molecules on the solid surfaces was usually assumed in MD simulations of OFMs. 26,27 However, Chia et al. 33 and Jaishankar et al. 34 used MD simulations to study the adsorption of surfactants on iron oxide surfaces from hydrocarbon solvents. Chia et al. 33 studied the competitive adsorption of the fuel additive 2-[(3-(dimethylamino)propyl)amino]methyl-4-dodecyl-6-methyl-phenol and ethanol from i-octane onto the α-Fe 2 O 3 (0001) surface at different temperatures. They showed that, while the free energy of adsorption of the additive was approximately two times larger than that for a single ethanol molecule, a single additive needs to displace up to five ethanol molecules to strongly adsorb onto the surface. Jaishankar et al. 34 studied SA adsorption from n-heptane and n-hexadecane onto Fe 3 O 4 (111), Fe 3 O 4 (001), and hydroxylated Fe 3 O 4 (001) surfaces at room temperature. They found that adsorption was stronger on the hydroxylated surface and adsorption from n-hexadecane was much slower than from n-heptane due to stronger alignment of the molecules on the surface. They used a combined approach using MD and molecular thermodynamic theory (MTT) to calculate adsorption isotherms, which agreed well with experimental data obtained with QCM. 34 The theoretical MTT model for the adsorption of non-ionic surfactants, originally developed by Nikas et al., 35 was applied in combination with MD for the first time by Sresht et al. 36 to study surface tensions at a water-air interface. Since then, it has been used to study adsorption in many systems, such as sodium dodecyl sulphate (SDS) on graphene from water, 37 SDS on α-Al 2 O 3 nanoparticles from water, 38 and organomolybdenum friction modifiers on borondoped diamond-like-carbon (DLC) from n-heptane. 39 Chia et al. 33 and Jaishankar et al. 34 extracted the potential of mean force (PMF) along the reaction coordinate using umbrella sampling (US) and the weighted histogram analysis method (WHAM). 40,41 Even using modern simulation software and high performance computing (HPC) architectures, this method requires hundreds of nanoseconds of simulation time to obtain converged PMFs. For example, Jaishankar et al. 34 used an overall simulation time of more than 1 µs to obtain the PMF describing the adsorption of SA onto hematite from n-hexadecane. A potentially more efficient method to obtain the PMF for surfactants at solid-liquid interfaces is the adaptive biasing force (ABF) algorithm. 42,43 The ABF method has been used to study adsorption at solid-liquid interfaces for a range of surfactants from water, such as SDS on graphene, 37 chlorogenic acid on polydimethylsiloxane, 44 lipoic acid on silver, 45 peptides on gold, 46 proteins on hydroxyapatite 47 and graphene, 48 ethylene glycols on calcite 49 and hematite, 50 and a wide range of organic molecules on carbon nanotubes, 51 graphene, 52 and graphene oxide. 53 In this study, we compare the adsorption of different OFMs from hydrocarbon solvents onto iron oxide surfaces. More specifically, ABF-MD is used to study the adsorption of SA, GMS and GMO from n-hexadecane and hydrogenated 1-decene trimer onto the α-Fe 2 O 3 (010) surface. We show that the OFM headgroup and tailgroup structure, as well as the solvent, significantly affect the adsorption behaviour. Once the adsorption PMF is obtained, the hard-disk area occupied by the molecules on the surface due to steric effects is estimated using separate annealing MD simulations. Subsequently, the relation between the bulk concentration of the OFMs in the solvent and the surface coverage, i.e. the adsorption isotherm, is determined using the MTT model. 34,35,54 Finally, the Jahanmir and Beltzer model, 20 together with High Frequency Reciprocating Rig (HFRR) friction measurements, are used to experimentally validate the calculated adsorption energies. 20,34

System Setup
The studied OFMs are SA, a carboxylic acid containing a saturated C 18 tailgroup, GMS, a glyceride mono-ester containing a saturated C 18 tailgroup, and GMO, a glyceride ester containing a C 18 tailgroup with a Z-double bond in the middle of the chain. While SA is commonly used as a prototypical OFM, it is no longer used in commercial formulations due to its acidic nature, and has been replaced with amines, amides, and esters, such as GMS and GMO. 12 The molecular structures of the OFMs are shown in Fig. 1a.
There has been some debate in the literature as to whether GMO and GMS adsorb  and reduce friction as intact glyceride esters or are hydrolysed to form the corresponding carboxylic acid that adsorbs, self-assembles, and reduces friction. 12,[55][56][57] Recent experiments have shown that the initial friction of GMO solutions is much lower than OA on silica surfaces. 25 Therefore, in this study, we consider the adsorption of intact GMO and GMS molecules.
For the hydrocarbon solvent, n-hexadecane was compared to hydrogenated 1-decene trimer. While n-hexadecane is commonly used as a model base oil in adsorption and boundary friction experiments, 24,25 hydrogenated 1-decene trimer is the primary constituent of commercial PAO4 base oil, a poly-α-olefin with a kinematic viscosity of 4 cSt at 100°C. 58 The molecular structures of the base oils are shown in Fig. 1b.
All of the systems for the PMF adsorption calculation were constructed using the Materials and Processes Simulations (MAPS) platform from Scienomics SARL, whereas the systems to estimate the hard-disk area on the surface were constructed using proprietary scripting.
When exposed to air, steel forms passivating iron oxide surfaces, such as hematite (α- 59 Two slabs of α-Fe 2 O 3 with dimensions (xyz) of approximately 55Å × 55Å × 12 A were used as the substrates, with the (010) plane exposed to the liquid phase. Periodic boundary conditions were imposed in the x-and y-directions. Initially, the two hematite slabs were separated by 120Å of vacuum space. For the adsorption calculations, a single OFM molecule was placed at the centre of the simulation cell, and either 700 n-hexadecane or 373 PAO molecules were randomly distributed between the slabs using the amorphous builder plugin in MAPS, see Figure 2. This gave a similar total number of C atoms (∼ 11200) for both systems. Two systems were used for the hard-disk area calculations (see Supporting Information). In the first set-up (dimerisation simulations), two OFM molecules were placed close to the surface of one of the hematite slabs, using the same system size and number of solvent molecules as in the adsorption calculations. In the second set-up (annealing simulations), one hematite slab was covered a with pre-formed monolayer of OFM molecules with different surface coverages. The number of OFM molecules included varied from 106 to 198 and the system size was approximately 55Å × 55Å × 100Å. This yields initial surface coverages, Γ = 3.5-6.5 nm −2 The base oil and OFM molecules were modeled using the all-atom long chain-optimised potential for liquid simulations (L-OPLS-AA) 61  It is noteworthy that recent DFT calculations have shown that strong molecule-surface interactions can occur for OFMs on iron oxide surfaces in the absence of solvent molecules. 29 DFT calculations have also been used to parametrize classical force fields to more accurately represent amide-surface chemisorption. 67 However, this task has not yet been executed for the OFM molecules of interest in this study.
Simulation procedure MD simulations were performed using the Large-scale Atomic/Molecular Massively Parallel Simulator (LAMMPS) software. 68 The velocity-Verlet algorithm 69 was used to integrate the MD equations of motion. A time step of 1 fs was used for all of the simulations and all bonds involving H atoms were fixed using the SHAKE algorithm. 70 Long-range electrostatic interactions were calculated using the slab implementation of the particle-particle, particle- mesh (PPPM) algorithm with a relative force accuracy of 10 -4 . 71 The PMFs were extracted using the ABF method 42,43 with the collective variables (colvars) module implemented in LAMMPS. 72 Before starting the PMF production runs, the systems were energy minimized using the conjugate gradient algorithm, with all of the hematite atoms being fixed in place. Subsequently, the system was equilibrated at 300 K and 101 kPa for 2 ns. The outer layer of atoms in the bottom wall of hematite were fixed in all directions, whereas the outer layer of the top wall was fixed only in the x-and y-directions, but allowed to move in the z-direction during the equilibration phase. These outer layers consisted of all of the Fe and O atoms that were 1Å or less from the bottom and the top of the simulation box respectively (440 atoms in both cases).
The pressure was controlled during the equilibration phase by applying a constant force to the outer layer of atoms in the top wall. After equilibration, the outer layers of both hematite slabs were fixed in place for the subsequent PMF production simulations. At the start of the equilibration phase, the velocities of the solvent, OFM and non-fixed hematite atoms were initialized using with a Gaussian distribution of velocities corresponding to a temperature of 300 K. The temperature was then controlled using the Nosè-Hoover 73,74 thermostat with a time relaxation constant of 0.1 ps. The thermostat was applied to all of the non-fixed atoms.
The PMF was calculated using the ABF method, an importance sampling algorithm developed by Pohorille and Darve. 42,75 The algorithm is based on the computation of the mean force along a certain collective variable (colvar), which is divided into small bins. The mean force is then cancelled out by an equal and opposite biasing force. This allows the system to move freely along that colvar, overcoming energy barriers and recovering ergodicity.
The mean force profiles converge once the relevant range of the colvar has been explored and sufficient samples have been accumulated within each bin. The convergence of this method is subservient to the efficient exploration of the phase space orthogonal to the chosen colvar.
In our case, this corresponds to the different orientations and conformations of the OFM molecules for each distance to the surface. Finally, the PMF is calculated from the mean force profile, since the average force exerted by the system onto the colvar is the gradient of the free energy profile along the colvar with a sign change: where A is the free energy profile as a function of the colvar ξ and F is the force exerted by the system onto the colvar. 42,75 In principle, ξ is a continuous variable and equation 1 applies for every value ξ 0 . In practice, the colvar is divided into small bins of width ∆ξ and samples are aquired for each bin during the simulation to compute the mean force in the interval (ξ 0 , ξ 0 + ∆ξ). Then, the mean value of those samples for each bin is used to compute the gradient and, through integration, the PMF. The review by review by Comer et al. 76 provides a more detailed description of the ABF method, together with recommendations and good practices.

Adsorption PMF
To estimate the PMF that captures the adsorption process, the colvar was defined as the projection onto the z-axis of the distance between the center of mass (COM) of the OFM polar headgroup and the COM of the topmost layer of atoms of the bottom hematite surface.
To accelerate convergence, the relevant range of the collective variable can be divided into several regions or "windows" to which the ABF algorithm is applied simultaneously. This approach is known as the stratification scheme. 77 The length of these windows should not be too small, as restricting the range of the colvar too much in a single simulation may lead to a slowdown in convergence or even a loss of ergodicity. 76 In other words, the simulation should visit the same configurations and in the same proportion as it would do using a much larger window, so that the artificial limit to the colvar explored in each window does not introduce a bias in the sampling.
In this case, the window closest to the surface should be large enough so as to avoid that the tail remains intermingled in the first maximum of the solvent layering while the colvar is being explored. Otherwise all the samples acquired for the first window would correspond to the same tail position, biasing the exploration of the phase space. Due to these considerations, the calculations for SA and GMS were run with three windows: from 1.8-2.3 to 20Å, 20 to 30Å and 30 to 40Å, since their saturated linear tailgroups are more extended than the kinked Z-unsaturated tailgroups in GMO. Therefore, four windows were used for GMO: 2.3 to 10Å, 10 to 20Å, 20 to 30Å and 30 to 40Å. The initial configurations for each window were obtained using steered molecular dynamics (SMD), 78 as implemented in the LAMMPS colvars package. 72 The force constant of the applied harmonic biasing potential was 10 kcal mol -1Å-2 , and the equilibrium distance was moved at a velocity of 4 m s −1 from the initial position of the molecule towards the surface. One initial condition was saved for each window and was equilibrated for 2 ns before applying the ABF algorithm to each window. The chosen width of the colvar bins was ∆ξ = 0.05Å.
The ABF algorithm was applied to each window until convergence was reached. As for other importance sampling techniques, the assessment of convergence is not trivial. 76 Several necessary conditions for convergence were monitored to determine a reasonable trade-off between computational cost and precision. These included the uniform increase of the samples across bins, the absence of further changes to the PMF, and the continuity of the mean force between different windows (see the Supporting Information). The mean force converged more slowly in the windows that were close to the surface due to the existence of several orientations and conformations that correspond to different values of the mean force. The total simulation times required to obtain converged PMFs varied between 165 ns and 750 ns for windows adjacent to the surface and between 75 ns and 150 ns for the others. The calculations using n-hexadecane converged more quickly than those with PAO4. This is consistent with the observations of Jaishankar et al., 34 who found that calculations converged more quickly for n-heptane than n-hexadecane, due to shorter relaxation times of the former.

Hard-disk area
To obtain a simulated adsorption isotherm using MTT, it is also necessary to estimate the area that each adsorbed molecule occupies on the surface. As explained by Nikas et al., 35 several difficulties arise concerning the estimation of the hard-disk area of each molecule on the surface due to steric effects. For instance, this area depends on the molecular conformations accessible to the adsorbed surfactant, which will depend on the surface coverage. As a result, the effective cross-sectional area will be coverage-dependent. In practice, it is usually sufficient to estimate the hard-disk area at high concentration and thus high coverage, since this is where the steric repulsion of adsorbed molecules limits further increases to the surface coverage.
OFMs are known to form vertically-orientated SAMs on iron oxide surfaces. 12,13 Therefore, the surfactant headgroup size will mostly determine the hard-disk area of the adsorbed molecule. Two approaches were used to estimate the hard-disk area occupied by the OFMs on the hematite surface. In the first approach, the hard-disk radius was estimated using the ABF method to calculate the PMF for the lateral interactions between OFMs. Two OFMs were adsorbed on the hematite surface and the projection onto the x-y plane of the distance between the headgroups was used as colvar. Subsequently, the cross-sectional area occupied by each OFM on the hematite surface was estimated using the position of the PMF minimum, following the approach used by Jaishankar et al. 79 for SA at the n-heptane-iron oxide interface and Sresht et al. 36 for polyethylene glycol at the air-water interface. The results obtained using this first method can be found in the Supporting Information. In the second approach, annealing simulations of preformed OFM monolayers were performed. A similar approach to that described by Zhang et al. 80 was used here. Different numbers of OFM molecules were placed close to the surface to form monolayers with a range of initial surface coverages. A subsequent equilibration of 2 ns allowed the molecules to adsorb onto the surface. Then the temperature was increased instantaneously from 300 K to 400 K and the system was equilibrated for 16 ns, which allowed the unbound molecules to desorb and escape the surfactant monolayer. Subsequently the temperature was reduced back to 300 The surfactants in the bulk are considered to be below the critical micelle concentra-tion (CMC), 82 thus the bulk chemical potential of the surfactants µ b is that of an ideal solution with no contributions coming from the interactions between surfactants. Imposing thermodynamic equilibrium µ σ = µ b , the surface coverage Γ and the bulk molar fraction X are related in a transcendental equation that can be solved using the bisection or Newton method. For more details, see the original publication 35 or the derivation from Sresht et al. 54 The final result reads: where X is the bulk mole fraction of monomer surfactant (assumed here to be the total concentration of surfactant without distinguishing between monomers, dimers, or micelles), Γ is the surface coverage, a is the hard-disk area, and ∆µ 0 ≡ µ σ,0 − µ b,0 is the free energy of adsorption at infinite dilution, i.e. with only one molecule. And B is the coefficient associated to attractive pairwise interactions between surfactant molecules, which won't be taken into account in our computational results. ∆µ 0 is calculated from the transfer free energy or difference in free energy between the bulk region and the adsorbed state. If the interface region is broad and not well defined, then the minimum of the PMF can be used directly for ∆µ 0 , considering that the surfactants are located at the Gibbs dividing surface, as proposed by Sresht et al. 54 Here, however, we follow the approach used by Jaishankar et al., 34 where the surfactants are not considered to be in a 2D plane above the hematite surface, but inside a 3D interface, where the different adsorbed states are taken into account. In this case, the expression for ∆µ 0 reads: where L = z 2 − z 1 , z 2 marks the point at which the PMF deviates from the bulk value, z 1 is the upper bound of the free energy well and z 0 is the lower bound of the free energy well. v s is the molecular volume of the solvent and ∆G(z) is the adsorption PMF.
Jaishankar et al. 34 considered an additional term that takes into account the entropic penalty of increasing the coverage as the chains are not free to rotate. This contribution has been analytically estimated and increases the potential energy of the surfactants at the surface 83 and is included here as well in the MTT model: Here N c = (n c + 1)/3.6 is what Nagarajan et al. 83 consider as an orientationally free segment In hydrocarbon solvents, carboxylic acids such as SA 84 dimerise, while glycerol esters, such as GMO 85 and GMS 86 can form reverse micelles. Previous MD simulations have shown that SA monomers and dimers exhibit equal affinity for iron oxide surfaces. 34 As micellation is not taken into account here, the results will be representative of systems below the CMC of the surfactants in the hydrocarbon solvent. 82 To our knowledge, the CMC has not been quantified for the systems studied here. Future work could use the same techniques adopted here to study the adsorption of reverse micelles on solid surfaces. 56,87 Experiments To validate the calculations of the adsorption free energy, friction was used as a proxy for surface coverage using the Jahanmir and Beltzer model. 20 The obtained surface coverage can then be related to the adsorption energy. Friction measurements were performed using a highfrequency reciprocating rig (HFRR) from PCS Instruments (Acton, UK). 88 The standard specimens were used, namely a 6 mm diameter ball (upper) loaded on a 10 mm diameter disk (lower). Both of the specimens were made from AISI E-52100 steel. A 2 N load was applied to the upper specimen, which followed a reciprocating movement of 10 Hz with a stroke length of 0.1 mm, giving a maximum velocity of 2 mm s −1 . The corresponding maximum Hertzian pressure was 0.8 GPa, whereas the mean was 0.5 GPa. High pressure and low speed were used to ensure that boundary lubrication conditions were maintained over most of the stroke.
This is important since it is here that the friction behaviour depends strongly on the OFM concentration and consequently the coverage. 12,24 Ambient temperature (300 K) was used and each test was conducted for 1 h to reach steady state friction. The friction coefficient was measured for pure base oil and for OFM concentrations between 10 -6 M and 0.02 M, for a total number of 7 HFRR tests. At the beginning of each test and once the contact was flooded, the ball and disc were held stationary and no pressure was applied during 5 minutes to allow for the formation of the OFM film, which is similar to the procedure used by Fry et al. 25 The two specimens were cleaned prior to each test with acetone and toluene.
The base oil used was PAO4 from ExxonMobil Chemical (Houston, USA) and the OFMs were GMO and GMS supplied by Sigma-Aldrich (Madrid, Spain) with a high purity grade (> 99%). The steady-state friction coefficient was determined averaging over the last 10 min of the experiment. Three repetitions were used for each OFM and concentration to reduce the noise in the HFRR data.
To validate the calculation of the adsorption energy at infinite dilution, the Jahanmir and Beltzer model was used. 20 This model relates the measured friction from the HFRR experiments with the fractional coverage. The model assumes that the contact is in boundary lubrication conditions (so friction is given by asperity contact and not the viscosity of the base oil) and that the friction coefficient is a linear function of the fractional coverage: where f is the measured friction coefficient, θ ≡ Γ/Γ max is the fractional coverage, f a is the friction coefficient at maximum additive coverage and f b is the friction coefficient with no additive present. Although we considered θ = 1 to determine f a from the HFRR measurements, complete coverage is not achieved in practice according to our computational adsorption isotherm results. This observation is in agreement with previous MTT-MD studies. 34 Once f a and f b are determined, the Jahanmir and Beltzer model can be applied to the friction data to obtain a relation between the OFM bulk concentration and the fractional surface coverage, to which an adsorption isotherm model can be fitted. Jaishankar et al. 34 used the Langmuir model for this purpose. However, according to Jahanmir and Beltzer,20 the Temkin isotherm model is better suited for these OFMs than the Langmuir isotherm model since it is able to account for inhomogeneous surfaces and lateral interactions. In the Temkin isotherm model the adsorption energy increases linearly with coverage (making the adsorption process less favorable): where ∆G 0 is the adsorption energy at θ = 0, α is a positive constant depending on the interactions of the molecules on the surface, and θ refers to the fractional coverage. As Representing the fractional coverage as a function of ln(X), α can be found from the slope m, whereas ∆G 0 is found from the intercept n as shown in equations 8 and 9. ∆G ads = −αn (9) Results and discussion Adsorption For both solvents, the bulk density in the centre of the film was in good agreement with experimental values. We obtained approximately 742±5 g/L for n-hexadecane and 778±20 g/L for PAO4, which is below the experimental values by 2 % and 4 %, respectively. At the end of the equilibration phase, the separation between the two walls was 12 nm for nhexadecane and 11 nm for PAO4. The interaction of the hematite surface with the base oil gave in both cases rise to a layering of the solvent molecules, as can be seen from the mass density profiles of the solvent molecules after equilibration in Fig. 3. The layering next to the wall is stronger for linear n-hexadecane than branched PAO4, in agreement with previous MD simulations. 89 The layers are slightly thicker for PAO4 than n-hexadecane.  be determined by the interactions of the OFM molecule with both the surface and the solvent.
All calculated adsorption PMFs (Figs. 4a, 4b, 4c) exhibit at least one minimum close to the surface. The PMF profiles also show energy barriers that correspond to the interaction of the tail with the layering of the solvent molecules close to the surface. There are only minor differences between PMFs from n-hexadecane and PAO4. The simulation time required to reach convergence of the PMF profile is much greater for PAO4 than n-hexadecane, due to its longer relaxation times. 34 The results for SA are in reasonable agreement with the ones from Jaishankar et al. 34 for the adsorption of SA from n-hexadecane onto Fe 3 O 4 (111). As the SA molecule approaches the α-Fe 2 O 3 (010) surface, the free energy oscillates with a period (∼0.5 A) that corresponds to the thickness of the solvent layer (Fig. 3). At OFM-surface distances of around 2Å, 1.5Å, 1.0Å, and 0.5Å, there are free energy barriers that correspond to the OFM moving through the solvent layers. The tailgroup is preferentially located in the maxima of the solvent mass density profile, and a small energy barrier has to be crossed for the tailgroup to move between these maxima. The height of these free energy barriers increases as the OFM gets closer to the surface, because the solvent layering is more pronounced.
The final barrier, corresponding to moving through the last solvent layer, is approximately  combinations. For SA, the barrier height is larger for PAO4 than n-hexadecane, despite the stronger layering for the latter (Fig. 3). This implies that the solvent layers are more difficult for SA to penetrate for branched alkanes solvents than linear ones. In practical terms, this means that SA adsorption will be slower from PAO4 than n-hexadecane at 300 K. 34 The improved ability for the saturated linear OFM molecule to penetrate the n-alkane of similar chain length could be the molecular mechanism of the chain matching phenomena observed experimentally. 90 For GMO, the free energy barrier heights are similar for both solvents. In between these barriers, there are local minima in the free energy profiles corresponding to outer-sphere complexes. These have also been observed in previous MD simulations of oxalic acid adsorption on rutile from water 91 and napthenic acid adsorption on calcite from water. 92 For SA, the headgroup was defined as the carboxyl group −COOH. Only one innersphere adsorption mode was observed in the PMF at 2.3Å above the surface, see Fig. 4a.
The −OH group of the carboxylic acid was oriented towards either one oxygen atom or a bridge between two oxygen atoms of the hematite surface with a OH···O distance of ∼ 1.

Hard-disk area
The hard-disk area was determined using annealing MD simulations. 80 Various monolayers of different coverages were created for the three OFMs, namely from 3.5 nm -2 to 6.5 nm -2 for SA and from 3.5 nm -2 to 5 nm -2 for GMO and GMS. The final coverage, or alternatively the final area per molecule, can be obtained by counting the number of molecules adsorbed at the end of the simulation. 81 The threshold distance to consider a molecule as desorbed was 4Å, measured from the oxygen of the hydroxyl group to the COM of the first layer of surface atoms. Increasing the threshold distance above 3.5Å did not change the results significantly, whereas decreasing it below this value underestimated the number of adsorbed molecules at low coverages, where all the molecules are assumed to be adsorbed, as shown in the Supporting Information. It is important to note that these results only give information regarding the minimum area due to steric effects that can be occupied by molecules adsorbed on the surface in a thermodynamically stable monolayer. It is not necessarily a realistic depiction of the maximum coverage that can be kinetically accessed in realistic timescales, since the OFM monolayers are preformed. Jaishankar et al. 34 suggested that the maximum accessible coverage can be illustrated both using the MTT and the random sequential adsorption (RSA) model. 95 The RSA model describes the irreversible adsorption of hard-disks on a uniform surface, which predicts a maximum coverage below the closed-packed maximum one, in which all of the surface would be covered. This is a result of irreversible adsorption and absence of restructuring between adsorption events, creating a jamming limit that was in agreement with the results by Jaishankar et al. 34 In reality the value of the effective hard-disk area occupied by each molecule on the surface  They explained these findings in terms of a reduced surface packing efficiency due to the presence of the kinked Z-alkene. However, in the annealing simulations presented here no significant difference is observed between both, suggesting that the observed differences in surface coverage by Wood et al. 21 are instead due to kinetic effects due to steric hindrance, rather than surface packing efficiency.

Adsorption isotherms from MTT
The main input parameters for the MTT model are the adsorption energy, the hard-disk area, and the number of orientationally-free carbon atoms n c in the tail. The value of the harddisk area for each OFM is taken as A hp from Table 1 n-hexadecane. 34 Moreover, the hard-disk area (Table 1) is also consistent with previous estimates of the molecular area of SA. 97,98 Consequently, we suggest that the rate at which coverage increases with OFM concen-      obtain an experimental estimate of the adsorption energy, see Figure 9. This approach yields and adsorption energy of -24.9 kJ mol −1 for GMO (α = 13.1) and -33.2 kJ mol −1 for GMS (α = 15.5). These values are the same order of magnitude as the computational results calculated through the adsorption PMF of each OFM, see Table 2. The agreement between the computational and experimental results is much better for GMS than for GMO and the α value found for both OFMs is the same. One would expect to find different α values if different chain packing efficiencies lead to reduced dispersion forces between adsorbed molecules in the case of GMO and thus lower coverage. This was the case with SA and isosearic acid reported by Jahanmir and Beltzer,20 in which the presence of the methyl side group decreases packing efficiency and thus dispersion forces. They found a different slope m and therefore a different α for both surfactants. The same authors, however, ascribed the differences in coverage between SA, elaidic acid and oleic acid to changes in chain dispersion interaction even though they found similar values of α for the three surfactants.
Overall, the differences between our computational adsorption energy (with a single molecule) and the experimental estimation fitting the Temkin adsorption isotherm for GMO can be understood through the steric hindrance for molecules adsorbing into partially formed monolayers. The kink in the Z-unsaturated tailgroup of GMO will lead to larger kinetic barriers for molecules to pass through the film and adsorb on the hematite surface. As well as affecting the kinetics of adsorption, this may also influence the thermodynamics. The steric barrier could also lead to a reduction in the adsorption energy for GMO (as predicted from the HFRR experiments) since on average, molecules are trapped further from the surface.
Another factor worth bearing in mind is that the PMF calculation is completely agnostic to both the interaction between molecules in solution. It is possible that reverse micelle formation is more favourable in GMO than GMS, 85,86 reducing the availability of molecules to adsorb on the hematite surface.  Table 2.

Conclusions
In this work, we have used MD simulations together with the MTT model to predict the adsorption isotherm of three OFMs in hydrocarbon solvents onto an iron oxide surface. First, we calculated the adsorption PMF for the three OFMs using MD and the ABF method, which proved to be more efficient than other comparable enhanced sampling approaches.
The computational results showed that the adsorption strength increased in the order SA < GMS < GMO, suggesting that the adsorption energy increases with the number of functional groups that interact with the surface. The stronger adsorption energy for GMO than GMS is due to the additional Z-alkene function group in the tailgroup. The different characteristics of the PMF have been interpreted in terms of headgroup-surface and tailgroup-solvent interactions. Both n-hexadecane and PAO4 give similar results, with the calculations using PAO4 or GMS taking longer to converge, due to slower transitions of the OFM tailgroup between solvent layers. Subsequently, the hard-disk surface area has been estimated using annealing simulations of pre-formed monolayer films at different surface coverages. The results show that the hard-disk area increases in the order: SA < GMS ∼ GMO. The adsorption energy and hard-disk area are used to calculate the adsorption isotherm using the MTT model.
The amount adsorbed at low coverage is mainly controlled by the adsorption energy and the maximum coverage by the hard-disk area. Lateral interactions between adsorbed molecules control how fast the adsorption isotherm reaches its asymptotic value. While one experimental data-set for SA can be adequately described with only the adsorption energy and hard-disk area, two other data-sets require the lateral interactions to be accounted for. The adsorption isotherm for GMS and GMO are similar, in contrast with previous experimental results regarding the packing of saturated and unsaturated surfactants. We propose that these differences in coverage that are found experimentally, are mainly due to differences in the kinetic barriers to the formation of high-coverage monolayers, with differences in ag- for the OFM with a Z-unsaturated tailgroup (GMO). We have shown that it is possible to predict the adsorption isotherm of saturated OFMs in silico. It is expected that the differences for the unsaturated OFMs can be accounted for by accounting for the decrease in adsorption energy on already partially formed monolayers.