The prediction of thermophysical properties at extreme conditions is an important applicationofmolecularsimulations. Thequalityofthesepredictionsprimarilydepends on the quality of the employed force ﬁeld. In this work, a systematic comparison of clas- sical transferable force ﬁelds for the prediction of diﬀerent thermophysical properties of alkanes at extreme conditions, as they are encountered in tribological applications, was carried out using molecular dynamics simulations. Nine transferable force ﬁelds from three diﬀerent classes were considered (all-atom, united-atom, and coarse-grained force ﬁelds). Three linear alkanes (n-decane, n-icosane, and n-triacontane) and two branched alkanes (1-decene trimer and squalane) were studied. Simulations were carried out in a pressure range between 0.1 and 400 MPa at 373.15 K. For each state point, density, viscosity, and self-diﬀusion coeﬃcient were sampled, and the results were compared to experimental data. The Potoﬀ force ﬁeld yielded the best results. 1 Introduction Thermophysical properties of lubricants at extreme pressures are highly important for the design of processes in mechanical engineering.1Especially the density and the transport properties determine the performance of the lubricant in the lubrication gap. The pressure in tribological systems with small contact areas are often higher than 1,000 MPa.2For such conditions, experimental data are rare due to high costs and complexity of the measure- ments.3Hence, reliable predictive models are needed. Molecular simulation is an attractive method for this. Molecular dynamics (MD) simulations have been used often for the prediction of ther- mophysical bulk properties of lubricants,4–11wetting, and interfacial properties on solid walls,12–15as well as for the simulation of tribological contact processes.16–23Simulations of lubricants usually aim at predicting transport properties, in particular the viscosity, using either equilibrium molecular dynamics (EMD)24–28or non-equilibrium molecular dynamics (NEMD) simulations.7,29–32Both methods are challenging: EMD simulations require a large computational eﬀort to compute the viscosity, especially for higher viscosities due to the slow convergence of the auto correlation function (ACF).33In NEMD simulations, the shear rate has to be small enough to reach the Newtonian regime which can require very long simulation times due to the bad signal-to-noise ratio at small shear rates.34As the viscosity becomes very large at high pressures,35the prediction of transport properties at high pressures is a challenging task with high computational costs. Special methods have been applied for tackling this task, such as the time decomposition method.33 MD simulations of lubricants require force ﬁelds that model the inter- and intramolecular interactions. Especially for the quantitative prediction of bulk properties, the choice of an adequate force ﬁeld is crucial. There is a large number of diﬀerent force ﬁelds that can be used for modeling lubricants.36–38The choice of the force ﬁeld is not trivial as diﬀerent force ﬁelds have diﬀerent strengths and weaknesses, and comparisons on equal terms are usually diﬃcult from the information that is available in the literature alone. The most 2 important attribute of the force ﬁeld is its the ability to predict the target properties with the desired accuracy. Another attribute is the type of the force ﬁeld, i.e. the approach that is taken to model the atomistic architectures and their interactions. This is decisive for the computational cost but also relevant for the choice of the simulation tool, as not all programs can handle all types of force ﬁelds. Three diﬀerent basic modeling levels are distinguished in the present work (cf. Fig. 1): all-atom (AA), united-atom (UA), and coarse-grained (CG) force ﬁelds. united-atom coarse-grainedall-atom Figure 1: Sketch of a squalane molecule modeled by the diﬀerent types of force ﬁelds. In AA force ﬁelds, every atom is assigned to an individual interaction site. They are also often called ’explicit hydrogen’ models. In contrast, UA force ﬁelds assign multiple atoms to a single interaction site. Usually, methylene and methyl groups are modeled by a single interaction site. The modeling by CG force ﬁeld goes a step further and chain segments of multiple heavy atoms are combined to a single interaction site. For AA force ﬁelds, two subtypes can be further distinguished: (i) reactive AA force ﬁelds, which explicitly model chemical bonds (that can form and break during the simulation), and (ii) classical AA force ﬁelds, which are built using intramolecular potentials that do not break, e.g. a harmonic potential. Nine force ﬁelds were applied in the present work for the modeling of ﬁve linear and branched alkanes, which are typical components of poly- -oleﬁnes (PAO) or used as simple model lubricants. PAOs are an important group of lubricants.39The substances include the three linear alkanes n-decane (C 10H22), n-icosane (C 20H42), and n-triacontane (C 30H62) 3 as well as two isomers of n-triacontane: 1-decene trimer and squalane (both C 30H62). The chemical structures of the ﬁve substances investigated in this work are depicted in Fig. 2. In the following, the substances are abbreviated as n-C10, n-C20, n-C30, TRI, and SQU (cf. Fig. 2). Hence, the study comprises alkanes with diﬀerent chain lengths as n-decane (n-C10) n-icosane (n-C20) n-triacontane (n-C30) 1-decene trimer (TRI) squalane (SQU)Name Structural formula Chemical formula C10H22 C20H42 C30H62 Figure 2: Overview of substances included in this study. well as diﬀerent structure for a given molecular formula. The molecule 1-decene trimer is the simplest (long branched) representative of poly-(1-decene) and it is an important component of PAO base oils.39Squalane (SQU) is frequently used as a model lubricant in experiments29,40and also an important substance in cosmetics.41The predictions for the density, the viscosity, and the self-diﬀusion coeﬃcient of these substances were systematically compared to experimental data. Moreover, the computational costs of the diﬀerent force ﬁelds were assessed. The study covers pressures up to 400 MPa, the temperature is always 373.15 K. The properties of lubricants at high pressure are of particular interest due to the high pressure conditions occurring in tribological applications.2Furthermore, the linear alkanes are common components of diesel fuels, whose properties at high pressure are also of high practical interest.42 Comparisons of diﬀerent force ﬁelds have been published for a wide variety of substances. They address the prediction of phase equilibria43,44as well as bulk properties like the den- sity,25,26,37,45–48isothermal compressibility,37,44and the speed of sound.45In some of the works, the ability of the force ﬁelds to predict thermophysical properties at elevated tem- peratures and pressures was explicitly addressed.47,49,50Studies explicitly addressing the 4 properties of lubricants are sparse. Lin et al. calculated the density and viscosity of a polyol ester (POE) lubricant at ambient pressure with three diﬀerent force ﬁelds and compared them to experimental data.48Ewen et al. compared ten UA and AA force ﬁelds in their ability to predict the density and the viscosity of n-hexadecane up to 200 MPa.26The vis- cosity and the self-diﬀusion coeﬃcient of branched alkanes as model lubricants have been studied by Kondratyuk et al.25A broad systematic comparison of force ﬁelds of diﬀerent types regarding their ability to predict properties that are important for lubrication, i.e. density and viscosity, up to high pressures was, however, not available up to now. This gap is closed by the present study. The diﬀerent force ﬁelds lead to diﬀerent computational costs of the molecular simula- tions. For example, the number of interaction sites per molecule varies for diﬀerent force ﬁelds and has a strong inﬂuence. The computational costs normally scale linearly with the number of interaction sites for systems without electrostatic interactions.51With AA force ﬁelds, the number of interaction sites is approximately three times higher compared to UA force ﬁelds for alkanes. As CG force ﬁelds combine multiple UA sites to one site, the computational are further reduced. Another advantage of coarser modeling of the molecular structure consists of preventing high frequency oscillations (e.g. CH stretching bonds), which enables a larger time step. Furthermore, the type of interaction potentials to be evaluated strongly inﬂuences the computational costs, for example multibody or electrostatic interac- tions increase the computational time considerably.52,53Overall, multiple factors inﬂuence the computational cost of a simulation with a given force ﬁeld and the costs can not be pre- dicted a priori. Presently, no systematic comparison of the computational costs of diﬀerent transferable force ﬁeld types is available. The following text is organized as follows: ﬁrst, the methodology of the study is de- scribed which includes details of the conducted simulations, the sampling of the observables, information on the applied force ﬁelds, and the methods applied for the evaluation of the results. Then, the results for the three examined properties, the density, the viscosity, and 5 the self-diﬀusion coeﬃcient, for all studied substances are presented and discussed. There- after, the computational costs of simulations with diﬀerent force ﬁelds are compared. Finally, conclusions are drawn. Methodology Simulation details and evaluation Equilibrium molecular dynamics simulations were carried out to calculate the density, the viscosity, and the self-diﬀusion coeﬃcient of the studied pure substances. A cubic simulation box with periodic boundary conditions was used. The number of particles was chosen such that there were at least 200 molecules or 4000 interaction sites. The ﬁnite size dependency of the simulations was assessed by simulations with diﬀerent numbers of molecules and extrap- olating to the thermodynamic limit ( Nmol). Details are given in the Supplementary Material. The results show that no signiﬁcant ﬁnite size dependency applies to the results. The time step and integrator were chosen depending on the force ﬁelds and are summarized in Tab. 1. For the AA force ﬁeld (excluding the reactive force ﬁelds), the RESPA multiple Table 1: Time step and integrator used in the simulations. force ﬁeld type time step / fs integrator AA 0.25 (bonds, angles), 0.5 (dihedrals), 1 (nonbonded) RESPA54 AA (reactive) 0.25 UA 0.5 Velocity-Verlet55 CG 10 timestep method was used to reduce the computational costs. Thereby, the simulation time was reduced by a factor of at least 8 compared to the classical Velocity Verlet integrator with a time step of 0.25 fs. Details are given in the Supplementary Material. The cutoﬀ radius was set to 14Å throughout. An analytical long-range correction for the pressure was applied for the dispersive interactions.55Rigid bonds were realized using the SHAKE algorithm.56The electrostatic interactions were calculated by the particle-particle particle-mesh algorithm.57 6 All simulations were carried out with the LAMMPS package (Version 3 March 2020).58The initial conﬁguration was built up by randomly inserting molecules into the box. For the reactive force ﬁelds, the molecules were inserted on a lattice to eliminate overlapping, which could lead to false molecule geometries as the geometry of the molecule is only given by the initial position. This insertion was done at a low density to avoid unphysical overlapping. Subsequently, the energy of the system was minimized by moving the particles to avoid large forces at the start of the MD simulations. After the energy minimization, the box volume was reduced to the target density within a short MD run of 100 ps. After this initialization phase, the actual simulations started. For each substance and force ﬁeld, ﬁve state points were studied. The pressures were p"r0:1;10;100;200;400xMPaand the temperature was T373:15 Kin all cases. The density was determined by NpT simulations. The Nosé-Hoover thermostat and barostat59–61 were used to control the temperature and the pressure in the simulations. The damping times comprised 100 timesteps for the thermostat and 1000 timesteps for the barostat, as recommended in the LAMMPS documentation.58Each NpT simulation had an equilibration time of 1 ns and a production time of 2 ns. To ensure that the simulation is suﬃciently long enough, longer simulations have been run. More details are given in the Supplementary Ma- terial. To reduce the statistical uncertainty, each NpT simulation was replicated ﬁve times and the density of each state point was calculated as the average of these ﬁve replica simula- tions. This density was then taken as input for NVT simulations to determine the viscosity and the self-diﬀusion coeﬃcient. To calculate the transport properties, 20 independent NVT simulations were carried out for each state point. The equilibration time was 1 ns. The production time varied from 2 to 10 ns and was determined individually for each state point and substance, and was chosen to obtain convergence of the correlation functions. Details are summarized in Tab. 2. 7 Table 2: Simulation times in ns for the diﬀerent substances and state points. substance p©MPan-C10n-C20n-C30, TRI, SQU 0.1 23 4 10 23 4 100 24 6 200 26 8 400 46 10 The viscosity was calculated using the Green-Kubo relations62,63given by  1 VkBTE  0J pt~ J ptd~ ; (1) whereVis the box volume, kBthe Boltzmann constant, Tthe temperature, tthe reference time, thetimevariable,and J pthenon-diagonalentriesofthestresstensor.55Thepressure tensor was sampled every 4 fsduring the NVT simulations. Each of the sampled time steps was taken as reference time tto achieve best possible statistics. The time decomposition method (TDM)33was applied in this work, which enables the reliable calculation of high viscosities by EMD simulations. Therefore, the viscosity time evolution  was calculated for all 20 NVT simulations during the production run. For the calculation of the ACF, a fast Fourier transform (FFT) technique was used.64Contributions from all six independent entries of the pressure tensor p were averaged to calculate the viscosity as proposed by Ref.65The cutoﬀ time cutfor the ACFs (to circumvent the ACF long-term noise) were determined from the standard deviations  as cut  cut0:4as recommended in Ref.33 Diﬀerent numbers cut  cutfor the cutoﬀ time have been compared in the Supplementary Material. The bootstrapping method33with 200 iterations was used to calculate the ﬁnal value and its uncertainty for the viscosity. Therefore, 200 combinations of the 20 simulations were created randomly whereas multiple appearances of individual simulations were allowed. For each combination, an average viscosity  was calculated and the function given in 8 Eq. (2) was ﬁtted to  whereis the long term limit.   11e 11  21e 2 11  2(2) Fig. 3 illustrates the TDM method for two state points showing the viscosity sampled during the production phase  . The results for all 20 simulation replica are shown. The for η final results p = 200 MPa τ cut = 1719 psη / 10 -3 Pa s 0 10 20 30 τ / ps 0 500 1000 1500 p = 0.1 MPa τ cut = 384 psη / 10 -3 Pa s 0 2 4 0 100 200 300 400 Figure 3: Illustration of the TDM method. Results for the sampled viscosity from 20 inde- pendent runs  over time of TRI atp0:1 MPa(top) andp200 MPa (bottom) (T373:15 Kfor both). Histograms of the corresponding bootstrapping results and the ﬁnal result of the bootstrapping method. two simulations have diﬀerent cutoﬀ times cut, which compensates diﬀerences of long-term noise behavior. On the right-hand side of Fig. 3, also the ﬁnal values and the histograms consisting of 200 values as obtained from the bootstrapping method are shown. The statistical uncertainties were computed from the histogram as the 95 % conﬁdence interval. Theindividualresultsfor  showinpartsawidescatteringthatincreaseswithtime, which is a well-known phenomenon.66The time required to reach a stable plateau of the viscosity diﬀers strongly which can also be seen by the diﬀerent cutoﬀ times. Overall, the TDM method provides relatively low statistical uncertainties and enables a reliable calculation of 9 the viscosity. The self-diﬀusion coeﬃcient Dwas calculated by the Einstein relation55 Dlim d d 1 61 NN = i1¶rit rit¶2; (3) whereNis the number of particles in the box and riis the position vector of center of mass of the molecule or particle i. A calculation of the self-diﬀusion coeﬃcient by the Green-Kubo relation and the TDM method was not feasible as it would require to write out the velocity vector of all atoms with a high frequency. The Einstein relation requires the positional data of all atoms less frequently. In this work, the positional data were written out every 10 ps. The term ¶rit rit¶2is the mean square displacement (MSD). It is possible to calculate the MSD either from the atomic or the molecular center of mass.67In this work, the molecular center of mass was used. The MSD was calculated by making use of the FFT method.64The self-diﬀusion coeﬃcient Dwas calculated individually for each NVT simulation using Eq. (3). By averaging the results of the 20 independent simulations, a single value for Dwas obtained for every state point. The statistical uncertainty of Dwas calculated as the standard error as DÙ ÛÛÛÛÛÛÚ1 NrepNrep = i1DiD2; (4) withNrep20. The self-diﬀusion is subject to ﬁnite size eﬀects that are not negligible.68 There are two methods to obtain the self-diﬀusion coeﬃcient in the thermodynamic limit (i.e. inﬁnite box size). The classical brute force approach is conducting multiple simulations with diﬀerent box sizes and extrapolating to the thermodynamic limit. This is not feasible in this work since a large number of state points and force ﬁelds are studied. Instead, the 10 analytical correction term proposed by Yeh and Hummer69was applied as DLDLkBT 6L; (5) whereLis the box length of the simulation, the viscosity, and 2:837297is a dimen- sionless constant (adopted from Ref.69). Force ﬁelds Tab. 3 gives an overview of the force ﬁelds studied in this work and their assignment to force ﬁeld classes. All force ﬁeld parameters used in this study and validations of the implementations are given in the Supplementary Material. Table 3: Overview of all force ﬁelds used in this work. name type refs.published OPLS AA701996 COMPASS AA711998 L-OPLS AA722012 ReaxFF AA (reactive)732001 AIREBO-M AA (reactive)742015 TraPPE UA75,761998 Potoﬀ UA77,782009 TAMie UA79,802015 MARTINI CG812007 All-atom force ﬁelds (AA) Inthisstudy,theAAforceﬁeldsCOMPASS,71OPLS,70L-OPLS,72ReaxFF,82andAIREBO- M74,83were used. The latter two force ﬁelds are so-called reactive all-atom force ﬁelds that allow bond breaking. However, no bond breaking was observed in the simulations carried out in this work, hence, no special distinctions from the other all-atom force ﬁelds are necessary here. All AA force ﬁelds except AIREBO-M include Coulomb interactions. The COMPASS force ﬁeld models the dispersive interactions by a 9-6 Mie potential. The intramolecular 11 interactions include bonds, angles, dihedrals, impropers as well as cross terms, which include mixed bond, angle, and dihedral contributions. The parameters of the COMPASS force ﬁeld were determined from a ﬁt to quantum-mechanical (QM) data (charges and intramolecular potentials) and to liquid state properties (dispersive potential).71In the OPLS force ﬁeld, the classical Lennard-Jones (LJ) 12-6 potential is used for modeling the dispersive inter- actions. The bond and angle parameters of the OPLS force ﬁeld were adopted from the CHARMM/22 force ﬁeld.70,84The torsional and the dispersive parameters were ﬁtted to QM data and diﬀerent thermodynamic and structural properties, respectively. The L-OPLS force ﬁeld adopts the majority of the OPLS parameters but dihedral and dispersive param- eters as well as the charges were reﬁned to improve the modeling of longer hydrocarbons.72 Additionally, the bonds between carbon and hydrogen were constrained to a constant dis- tance in the L-OPLS force ﬁeld. It is known that the L-OPLS provides better predictions for bulk properties compared to the OPLS force ﬁeld.26As the OPLS force ﬁeld is still widely used in its original version, it was also included in this study. The ReaxFF82and AIREBO-M74,83force ﬁelds include explicit calculations of the bond order, which enables them to model chemical reactions.52TheCHON-2017_weak version of the ReaxFF force ﬁeld73was implemented in this work. It includes parameter adjustments thataimatimprovingthedensitypredictionofthecondensedphase.73TheAIREBO-Mforce ﬁeld74was derived from the AIREBO force ﬁeld83by replacing the LJ potential interactions by a Morse potential. Its parameters were adjusted to layer spacing of graphite up to pressures of 14 GPa and to QM data of small alkanes. Both reactive force ﬁelds were only applied to n-decane (n-C10) in the present work due to their extremely high computational costs. United-atom force ﬁelds (UA) The TraPPE,75Potoﬀ,77and TAMie79UA force ﬁelds were used in this study. The TraPPE force ﬁeld uses the classical 12-6 LJ potential. These LJ parameters were ﬁtted to critical 12 temperatures and saturated liquid densities. The Potoﬀ and the TAMie force ﬁelds use n-6 Mie potentials for the non-bonded interactions. In both force ﬁelds, the repulsive parameter nof the Mie potential was additionally adjusted in the ﬁtting procedure. For the TAMie and the Potoﬀ force ﬁelds, the non-bonded parameters were ﬁtted to saturated liquid den- sities and vapor pressures. The cutoﬀ distance of the Potoﬀ force ﬁeld was set to 10 Å as recommended.77In contrast to the AA force ﬁelds, there are no Coulomb interactions for UA models of hydrocarbons as all interaction sites are electrostatically neutral. In this work, ﬂexible bonds between the interaction sites were applied. Therefore, a harmonic potential was used with the given bond length as equilibrium bond length. The energy parameter of the bond potential was taken from the OPLS force ﬁeld. Despite the fact that the TraPPE, Potoﬀ, and TAMie force ﬁelds were originally developed using rigid bond lengths, it is com- mon practice to adapt them with ﬂexible bonds.85–87Moreover, it was shown that the ﬂexible bonds do not inﬂuence the results if the bond force constant is well chosen.49,88The angle and dihedral potentials are the same for all three UA force ﬁelds as they were mainly adopted from the OPLS-UA89and OPLS-AA70force ﬁelds. Coarse grained force ﬁeld (CG) In this study, the MARTINI CG force ﬁeld was used. There are two approaches to param- eterize CG force ﬁelds: bottom-up by derivation from a ﬁner resolved force ﬁeld (e.g. AA) and top-down by ﬁtting parameters directly to experimental data.90The MARTINI force ﬁeld81is based on the top-down approach. It was parameterized to the free energy of hydra- tion and vaporization as well as the partitioning free energies between water and multiple organic substances.81A four-to-one mapping is applied by the MARTINI force ﬁeld in gen- eral, which means four heavy atoms (all atoms except hydrogen atoms) are fused to a single interaction site. The force ﬁeld consists of the LJ potential for non-bonded parameters and harmonic potentials for the angle and bond interactions. There are no dihedral potentials. The mapping of molecules to CG beads is not strictly deﬁned. There are several methods to 13 create CG models from a given molecular topology.91–93In this study, the automatic map- ping algorithm from Potter et al.93was applied. As the mapping may in general have an inﬂuence on the force ﬁeld predictions, the MARTINI results should be considered as results obtained from "MARTINI + Potter mapping". For brevity, only "MARTINI" is used in the following. The bond lengths of the created models were adjusted as suggested in Ref.93 for the prediction of thermodynamic bulk properties. The MARTINI 281version was used in this work instead of the recently published MARTINI 3 force ﬁeld94as the MARTINI 2 version has been widely applied in recent years – and as diﬀerent auxiliary tools are available for MARTINI 2. Data analysis The simulation results for the density, the viscosity, and the self-diﬀusion coeﬃcient were analyzed in a consistent way as described in the following. As a basis, correlations to all available experimental literature data were developed. An overview of the available data for the ﬁve considered substances and the three considered thermophysical properties is given in Tab. 4. The functions, parameters, and accuracy of the correlations are reported in the Table 4: Overview of available experimental literature data (