Impact of chosen force fields and applied load on thin film lubrication

The rapid development of molecular dynamics (MD) simulations, as well as classical and reactive atomic potentials, has enabled tribologists to gain new insights into lubrication performance at the fundamental level. However, the impact of adopted potentials on the rheological properties and tribological performance of hydrocarbons has not been researched adequately. This extensive study analyzed the effects of surface structure, applied load, and force field (FF) on the thin film lubrication of hexadecane. The lubricant film became more solid-like as the applied load increased. In particular, with increasing applied load, there was an increase in the velocity slip, shear viscosity, and friction. The degree of ordering structure also changed with the applied load but rather insignificantly. It was also significantly dependent on the surface structure. The chosen FFs significantly influenced the lubrication performance, rheological properties, and molecular structure. The adaptive intermolecular reactive empirical bond order (AIREBO) potential resulted in more significant liquid-like behaviors, and the smallest velocity slip, degree of ordering structure, and shear stress were compared using the optimized potential for liquid simulations of united atoms (OPLS-UAs), condensed-phase optimized molecular potential for atomic simulation studies (COMPASS), and ReaxFF. Generally, classical potentials, such as OPLS-UA and COMPASS, exhibit more solid-like behavior than reactive potentials do. Furthermore, owing to the solid-like behavior, the lubricant temperatures obtained from OPLS-UA and COMPASS were much lower than those obtained from AIREBO and ReaxFF. The increase in shear stress, as well as the decrease in velocity slip with an increase in the surface potential parameter ζ, remained conserved for all chosen FFs, thus indicating that the proposed surface potential parameter ζ for the COMPASS FF can be verified for a wide range of atomic models.


Introduction
Hexadecane lubricant, which is widely used in many tribological systems, has been analyzed in several theoretical studies over the last two decades, particularly in studies on steel surfaces [1][2][3][4][5][6][7][8][9][10]. Under confined conditions, thin hexadecane film behaves differently from its bulk form. Theoretical investigations revealed the formation of layering and slip phenomena at solid-liquid interfaces [2, 6−8], which depend on many variables such as surface roughness [2,6,7], lubricant additives [4], surface structure and film thickness [2], as well as applied pressure [10]. Moreover, an increase in the solid-liquid adhesive strength significantly reduced the slip phenomenon at the interface [4,11,12], whereas the slip occurred internally between the molecular layers within the lubricant [11]. Zheng et al. [10] determined the superior anti-wear performance of the long alkane chains and concluded that they behaved more solid-like than other shorter alkanes.
Although iron oxides are commonly found on mild steel strip surfaces during cold rolling, studies on the lubrication performance of hexadecane on steel surfaces at the atomic level remain inconclusive.
Although several studies were conducted on the confined shear of hydrocarbons between metal surfaces, only a few of them considered the metal oxide surfaces in molecular dynamics (MD) simulations [1,2,4]. To characterize the adhesive strength between alkane and metal/metal oxide surfaces, Manias et al. [13] and Steven et al. [14] proposed a factor k, which varies in a wide range, to present the different interaction strengths between the alkane and the confined surfaces. Berro et al. [15] and Savio et al. [1,2] adopted this factor, which is based on the energy difference per CH 2 group obtained from experimental measurements of alkanes adsorbed on Ni [16], Ag [17], and Au [18] surfaces, to compare the tribological behavior of alkanes between different metal (Au, Fe) and metal oxide (CuO, Fe 2 O 3 ) surfaces. Ta et al. [19] presented a more practical adhesive strength between hydrocarbons with iron and iron oxide surfaces by deriving an interfacial force field (FF) from ab initio calculations.
In MD simulations, different conventional and reactive potentials have been applied to hydrocarbons. However, studies are yet to be carried out to ascertain their influence on tribological performance and degree of layering order. In a few notable attempts, the non-equilibrium MD simulations indicated that the condensed-phase optimized molecular potential for atomic simulation studies (COMPASS [20]) FF predicts a higher zero-shear viscosity of hexadecane than assisted model building with energy refinement (AMBER) 96 [21,22]. A benchmark was also set in terms of density and viscosity by Ewen et al. [23], who showed that united-atom (UA) FFs can accurately reproduce experimental density, whereas viscosity was significantly underpredicted compared to the all-atom (AA) model. They explained that the significant underestimation of viscosity by UA FFs suggested that the molecules moved past one another more easily, which indicates that the coefficient of friction (CoF) calculated for confined systems using UA FF will also be under-predicted. Although classical FFs can accurately model the physical processes of molecular systems under certain conditions, their potential functions are unable to reproduce the chemical reactions of these systems. Additionally, these FFs are generally poorly transferable when the properties of interest were not directly parameterized.
In this study, we considered the effects of the applied load and interfacial crystalline structure of iron compounds on the molecular structure and lubricability of hydrocarbon tribofilms. This study also aims to analyze the effect of available conventional and reactive potentials in predicting the tribological properties of hexadecane. The resulting article is organized in the following order. First, Section 2 presents the implemented methodology and relevant parameters for MD modeling. After that, Section 3 discusses the influence of loading pressure and crystalline structure on the rheological properties of hexadecane under confined shear conditions, as well as the impact of chosen FFs on the ordering structure and tribological properties of hydrocarbons. Finally, a conclusion is presented in Section 4.

Methodology
The MD simulation was performed using an AA model, in which the methyl groups of CH 2 and CH 3 were presented with the explicit atoms. COMPASS FF was applied to the lubricant [20], whereas Buckingham FFs were applied to the FeO and Fe 2 O 3 surfaces. The embedded-atom method/Finnis-Sinclair potential (EAM/FS) was applied to the Fe (100) surface [24]. The confined shear process, as illustrated in Fig. 1, was implemented using the large-scale atomic/ molecular massively parallel simulator (LAMMPS) code [25], while the initial configuration was constructed in Material Studio 7.0 (MS 7.0). The parameters for the Buckingham potential, as reported by Guillot et al., are presented in Table S1 in the Electronic Supplementary Material (ESM) [26,27]. The interactions between surface atoms with fluid atoms, as well as between different lubricant molecules, were determined using a Lenard-Jones 9-6 (L-J 9-6) potential with a cut-off distance of 12.5 Å combined with a long-range Columbic interaction. This long-range Columbic interaction adopted the particle-particle particle-mesh (PPPM) method with a desired relative error of 10 -4 in forces, and the slab technique with a ratio of 3, which is the ratio between the total extended volume used in the 2D approximation and the volume of the simulation domain [28]. The combination rules for this L-J 9-6 potential were applied to the interaction between the lubricant atoms, as well as between the lubricant and surface, which are expressed as Eqs. (S1) and (S2) in the ESM. To investigate the role of the applied load, COMPASS FF was used to describe the atomic interactions between hexadecane molecules. Theoretical MD simulations indicated that this FF applied to several substances and was commonly used to investigate the adsorption energy between organic molecules on metal oxide surfaces [29−32]. The functional forms and expressions for this FF are described in Table S2 in the ESM while the corresponding parameter values can be found in Ref. [20].
Interfacial interaction is a critical factor that affects the structure of organic molecules on the surface, as well as its rheological properties under confined conditions. Based on the adsorption energy between ethane and the iron surface, measured by Govender et al., the interaction between alkanes and iron was derived from ab initio calculations, as shown in Fig.  S1 in the ESM [33]. The procedure for calibrating this interaction parameter was implemented using their (Govender et al. [33]) work with L-J 9-6 potential. This parameterization was carried out at an absolute zero temperature. The ethane molecule and two upper iron layers were relaxed, whereas the rest were fixed. To obtain adsorption energy values consistent with Govender's result, the well-depth energy of iron ( Fe ) for L-J 9-6 interactions between the surface and lubricant was maintained in the range of 0.02-0.05 eV. The dependence of the adhesive strength between the ethane and iron surface on the well-depth energy  Fe is presented in Fig. S2 in the ESM, which shows that the  Fe value of 0.0425 eV appropriately represents the adsorption energy,   ad 0.23 eV, E obtained from previous ab initio calculations. The non-bond interaction parameters between the iron oxide surfaces and hexadecane ( Table S3 in the ESM) were obtained from the COMPASS FF and semi-ionic oxide model proposed by Zhao et al. [34].
As shown in Fig. 1, the molecular model was constructed with hexadecane allocated between the FeO (100) surfaces with a system domain of approximately 35 Å × 35 Å × 60 Å. Periodic boundary conditions were also applied in the lateral directions. To avoid self-interaction with its periodic images, the simulation domain was assigned a width twice the molecular chain length of hexadecane. This chosen system domain was similar to those presented by Berro et al. [35] and Jabbarzadeh et al. [36]. Each surface, cleaved from its regular crystalline lattice, had a thickness of approximately 10 Å, whereas the initial thickness of the lubricant was 40 Å. Although the Fe 2 O 3 (001) surface has three chemically distinct terminations, in this study, it was set up by a single iron termination because this termination had the most stable surface configuration [37−39]. Details of the domain sizes and number of hexadecane molecules for each confined shear model are presented in Table S4 in the ESM.
The established model underwent tribological conditions, in which an applied pressure was introduced to the upper surface. To monitor its influence on the rheological and structural properties of hexadecane, this applied pressure varied in the range of 100− 500 MPa. This MD simulation includes three main stages. The model was relaxed for 0.5 ns to neglect poor atom contacts during the first relaxation stage. Subsequently, a dynamic stage was introduced in the next 0.5 ns with an applied pressure on the top region to compress the model, while a fixed constraint was applied at the bottom region. To continuously conduct the frictional heat generated owing to compression and achieve an equilibrium condition, the temperature of the system was kept thermostatic at 300 K at the thermostatic layers on the surfaces. The NVT ensemble with a damping coefficient of 100 fs was used for this ensemble. In contrast, the NVE ensemble was applied to lubricant, as it was subjected to a constantly applied pressure. Finally, in the confined process, the surfaces were moved in opposite directions along the x-direction at a speed of 10 m/s for 10 ns, whereas the constraints of loading pressure, NVT ensemble for the surfaces, and NVE ensemble for the lubricant were retained. A time step of 1.0 fs was applied to this MD simulation.
To analyze the influence of the selected FFs, different potentials, from the simplest UA model to an advanced reactive FF, were used. The UA potential employed in this work was the optimized potential for liquid simulations of united atoms (OPLS-UA) derived by Jorgensen and Tirado-Rives. [40], whereas COMPASS [20], adaptive intermolecular reactive empirical bond order (AIREBO) [41], and ReaxFF [42] were used for the AA models. Note that the applied FFs for the surfaces and their interactions with hydrocarbons were similar in all cases to be solely considered for the role of the lubricant model [22]. A comparison between these FFs was carried out based on the analysis of molecular structure, velocity gradient, shear stress, and slip ratio.

Influences of the applied load
When the tribosystem was subjected to shear forces and applied pressure, the lubricant density increased gradually with the simulation time. Figure S3 in the ESM shows the time evolution of the atomic mass density, which oscillates during the first 100 ps owing to the instant impact of the applied load and becomes more stable after 11 ns. This figure also reveals that although an increase in applied pressure induces a subsequent increase in oscillation amplitude at the beginning of the loading process, the oscillation amplitude becomes more stable during the sliding process.
The average lubricant density in Fig. 2(a) indicates a linear correlation between the lubricant density and the applied load, as a higher pressure induces an increase in the density of confined hexadecane. This observation agrees with previous investigations by Martini et al. [8,10]. The reduction in film thickness due to the applied pressure ( Fig. 2(b)) could be the reason for the density increase in Fig. 2(a). However, a slight difference appeared among the averaged lubricant densities on different surfaces. This discrepancy may be due, in part, to the difference in the domain size of each surface using the same dimensions, number of lubricant molecules, and applied load for this comparative analysis. Additionally, the difference in the lattice constant for each surface material resulted in an apparent difference in surface dimensions (Table S4 in the ESM). However, the main source for this discrepancy is the difference in the ordering structure depending on the surface type.
Under relative sliding, kinetic energy was transferred into the fluid via the adhesive forces between the surface and lubricant [43,44]. This resulted in the formation of layers or ordering behavior of molecules adjacent to the surfaces. Figure 3 indicates that  although this layering density profile oscillates with the largest amplitude at the solid-fluid interfaces, its intensity gradually attenuates toward the middle regime of the lubricant film. This observation is consistent with Refs. [1,2,7,10,36,45,46]. The thickness of the film decreases with increasing applied pressure, which results in a decrease in the number of layers. In particular, when the applied load was increased beyond 300 MPa, the number of ordering layers was reduced from eight to seven for the FeO,  One of the interfacial phenomena that has attacted considerable interest in thin-film lubrication investigations is the interfacial slip. This slip was determined as the velocity difference between the effective speed applied on the surface and an apparent one calculated from the lubricant layer at solidluid interfaces. The velocity profiles across the film thickness for each surface represented in Fig. 4 demonstrate that most of the velocity slip is observed for the FeO surface, whereas it is insignificant for To quantify the slip phenomenon, a dimensionless parameter s is measured based on the obtained apparent app ( )


and effective ( eff  ) shear rates, as expressed in the following formula: The apparent shear rate is defined as app    Figure 5 shows the linear relationship between this slip parameter and applied pressure. In addition, this slip also depends strongly on the crystalline structure in the order Fe 2 O 3 (001) < Fe (100) < Fe 2 O 3 (012) < FeO (100). Fillot et al. [48] proposed that the slip varies linearly with the applied load. Therefore, a linear function was fit for the obtained MD results as expressed in the following equation: where p is the applied pressure. The fitting results for the different surfaces presented in Fig. 5    indicate that the onset of the slip could start at a small applied load for FeO, Fe, and Fe 2 O 3 (012) surfaces, although this phenomenon only occurs when an applied load larger than 0.06 MPa is imposed on the Fe 2 O 3 (001) surface. This observation suggests that a practical implication to this phenomenon is that the thin film lubrication of hydrocarbon could be optimized by utilizing the smoother surfaces more than Fe 2 O 3 (001), which is generally considered as the most stable facet of α-Fe 2 O 3 [49,50]. It is also commonly found on steel surfaces of the high load and elevated temperature [51,52]. Moreover, it also predicts that the full velocity slip at an applied load of 0.5 GPa is for the FeO (100) Fig. 5(b) shows a contrasting tendency in that the effective shear rate decreases almost linearly with an increase in the applied load. Friction is a critical component commonly considered in tribological investigations because it reflects the lubricity of the lubricant. In the current MD simulation, the shear stress was calculated as a fraction of frictional force, which was evaluated by the sum of the lateral forces on the surface exerted by the lubricant via the pairwise interfacial interactions over the surface area. This shear stress value was time-averaged during the last 5 ns for all cases, and the linear relationship between the calculated shear stress and an applied load is presented in Fig. 6. However, the figure shows a remarkable discrepancy in the CoF value for each surface material. Particularly, the respective CoF values of 0.009, 0.055, 0.027, and 0.072 were found for FeO,    (012) surfaces were smoother, therefore, they could generate a lower CoF value than that of Fe 2 O 3 (001). This result is in good agreement with that determined by Berro et al. [4] and Savio et al. [1].
Although FeO (100) also has a smooth surface similar to Fe (100), we observed that the CoF value of FeO (100) is lower than that of the pure iron surface. This is an interesting observation and an explanation for this phenomenon could be the commensurability of hexadecane with different surfaces. Thomson and Robbins [53] and Savio et al. [1] stated that a greater commensurability between a surface and fluid could lead to a greater enhancement of the interfacial slip. The current observation agrees with their studies because the gap between surface atoms on the FeO (100) surface was shorter than that on the Fe surface. In fact, gaps of 2.866 and 2.1865 Å were determined for Fe (100) and FeO (100) surfaces, respectively.
Lubricant viscosity is a critical rheological property in thin-film lubrication. Under the confined gap, this physical component is different from the bulk phase owing to the ordering rearrangement of the lubricant molecules and the change in density. The effective viscosity  eff of hexadecane was calculated to observe the effect of the adhesive strength and surface properties on this rheological property. The obtained results are presented in Fig. 7, which shows that the viscosity of hexadecane increases significantly with the applied load, especially in the case of FeO (100) and Fe 2 O 3 (012) surfaces. Furthermore, this effective viscosity exhibits a large increase compared with the bulk value, thus implying that the hexadecane lubricant behaves more solid-like under extreme conditions. The results also show that the solid-like behavior is more pronounced on FeO and Fe 2 O 3 (012) surfaces in comparison with Fe and Fe 2 O 3 (001). The dependence of  eff on eff  or p can be described as a power function in the following expression: where x is the effective shear rate or applied load, 0  is the effective viscosity at zero shear rate or zero-applied load, a and b are the parameters. The fitting results are plotted in Fig. 7.

Influence of force fields
The average density of C 16 H 34 during the confined shear simulation was calculated, and the results were presented in Table 1. The data reveal that classical FFs of OPLS UA and COMPASS result in higher densities than those obtained from the reactive FFs of AIREBO and REAXFF. COMPASS exhibits the best agreement with experimental values of 0.77 g/cc [54], whereas other FFs show larger deviations up to 7%.  It should be noted that a larger number of surfaces have been used in this case to improve the reliability of our observations on the effect of used FFs on the hydrocarbon lubricating performance. The obtained results are presented in Figs. 8−10 for different surfaces, which indicate that the variation in densities across the film thickness for the UA model are much higher than those obtained from the AA models. Remarkably, in the middle regime of the lubricant, a significant decrease in density peaks is observed for AA models; in contrast, the decrease is insignificant for the UA model. These observations reveal in general, that the molecular structure of hydrocarbon or layering property, in particular, depends on a chosen FF.
These results indicate that the FFs adopted to influence the tribological performance of the thin film.       In Ref. [22], we discovered that this slip phenomenon was mitigated by an increase in the surface potential parameter (ζ (kcal/mol/Å)), which is a general formula used to characterize the influence of surface properties, including surface energy, corrugation force, and commensurability height. Additionally, the ζ values increased in the order: Fe (110) < FeO (100) < Fe (100) < Fe2O3 (012) < Fe (111) < FeO (110) < Fe (111)(010) < Fe2O3 (001) < FeO (111), with respective values of 1. 28, 2.90, 8.32, 10.45, 18.59, 20.3, 77.36, 104.92, and 154.4 kcal/mol/Å [22]. To verify the transferability of our theory to different atomic models, the dependence of shear stress τ xz (MPa) and velocity slip ratio s on ζ was analyzed, and the results are presented in Fig. 14. The results show that the current trends are consistent with previous observations in which τ xz increases with ζ, whereas s decreases with an increase in ζ.
However, the shear stresses vary significantly between different FFs, with the smallest values within 20 MPa for AIREBO FF over various surfaces, whereas COMPASS FF exhibits results that are two to four times larger than AIREBO. Given the same surface properties and interfacial interactions between the surfaces and lubricants, the model of the lubricant could likely be the root cause of this difference. The less ordered structure, larger velocity gradient, and smaller shear stress reveal that AIREBO behaves more liquid-like than other FFs. In contrast, despite the OPLS-UA model producing smaller shear stresses than COMPASS and ReaxFF, the stronger layering structure (Figs. 8−10) and small velocity gradients  verify that this kind of FF is more solid-like than the others. The error bar for ReaxFF is larger than that of other force fields owing to the larger deviation in the viscous friction within the lubricant.
With an increase in ζ, the velocity slip when the lubricant is confined decreases. A similar concern has been aimed at shear stress, which is a critical property in lubrication. The determined values in Fig. 14(a) indicate an increase in shear stress with ζ. Although the figure exhibits some fluctuation, the relationship presented is approximately linear. This observation is not only found in previous studies on COMPASS FF, but is also applicable to other atomic models in the current work. This indicates that our proposed ζ function can be considered as  a general formula for characterizing the influence of surface properties, including surface energy, corrugation force, and commensurability height, on the tribological performance of lubricants at the atomic level.
The average lubricant temperature was also assessed and is presented in Fig. 14(c), which shows that the temperature generated within the lubricant by frictional heat is an independent function of the surface properties on smooth surfaces. As illustrated in Fig. 1, the thermostat with the NVT ensemble was only applied to the inner layer of the surface, such that temperature was induced from kinetic energy owing to sliding motion. The lowest temperature range of 230-250 K is determined for the OPLS-UA model. The higher 280−285 K values were evaluated for COMPASS FF. Reactive potentials result in larger temperature ranges of 295−320 K and 320−330 K for AIREBO and ReaxFF, respectively. The lower lubricant temperature observed for the OPLS-UA model is a direct consequence of the velocity slip at the interface because this slip reduces the kinetic energy transferred from the surface to the lubricant. The less transferred momentum from using the OPLS-UA model reduces the temperature of the lubricant to become smaller than other FFs. In essence, the atomic model provides a stronger layering structure and a smaller velocity gradient that will decrease frictional heat and lubricant temperature.

Conclusions
MD simulations were implemented to investigate the effects of applied load, chosen atomic models, and surface structure by adopting different crystalline structures of iron and iron oxides. The COMPASS force field with an AA model was used for hexadecane to study the influence of the applied load on the degree of ordering, shear stress, and velocity slip, which was addressed; whereas EAM and Buckingham potentials were applied to the surfaces. In addition, different potentials have been applied as lubricants to validate the transferability of our previously proposed surface potential parameter, which characterizes the influence of surface energy, corrugation force, and commensurability. The simulation results revealed as following: 1) The lubricant film became more solid-like with an increased applied load. In particular, there was an increase in the velocity slip, shear viscosity, and friction as the applied load increased. The degree of ordering structure also changes with the applied load but rather insignificantly. It is also remarkably dependent on the surface structure.
2) The chosen FFs significantly influence the performance of lubrication, rheological property, and molecular structure. With the smallest velocity slip, degree of ordering structure, and shear stress, AIREBO potential results in more liquid-like behavior compared with OPLS-UA, COMPASS, and ReaxFF. In general, classical potentials such as OPLS-UA and COMPASS behave more solid-like than reactive ones. Furthermore, the lubricant temperatures obtained from OPLS-UA and COMPASS are much lower than those of AIREBO and ReaxFF owing to their solid-like behavior.
3) The increase in shear stress and decrease in velocity slip with an increase in parameter ζ are still conserved for all chosen FFs, thus indicating that the proposed surface potential parameter ζ for COMPASS FF can be validated for a wide range of atomic models.