{'fraction_non_alphanumeric': 0.03903464667702715, 'fraction_numerical': 0.02651064996800439, 'mean_word_length': 4.650826446280992, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Comparing different force fields by viscosity prediction for branched alkane at 0.1 and 400 MPa

Shear viscosity is one of the key subjects of molecular modeling studies since this quality is used in the development of lubricants. In this paper, molecular dynamics methods are used to compute viscosity at 298 K and 0.1 MPa, 400 MPa for 2,2,4-trimethylpentane for which experimental data are available. Three different force fields are compared by the ability to predict the shear viscosity coefficient at these two pressures.


Introduction
Modern industry is strongly interested in studying the properties of liquid hydrocarbons, since they are components of the lubricants, insulation and fuel liquids. The properties that play a major role in the process of elasto-hydrodynamic lubrication (EHL) are the shear viscosity dependence on pressure and the pressure viscosity coefficient [1]. The classical film thickness formulas require the latter property but it is not always clear how to define it [2]. The viscosity dependence on pressure determines the shape of friction curve (friction coefficient vs. sliding speed) of oil films [3]. However, such dependencies are often not measured experimentally in the EHL analysis and the demand for such kind of data is growing in recent years.
The molecular dynamics (MD) techniques proved useful for the estimation of such properties [4][5][6]. Allen in 1997 [7] made an attempt to predict the 2,2,4-trimethylpentane shear viscosity at 0.1 MPa from MD simulations with OPLS-AA [8] and OPLS-UA [9] force fields using NEMD. The OPLS-AA potential value lies very close to the experimental data at 0.1 MPa. Lahtela and coathours have studied how branching influences the viscous properties of eicosane isomers using NEMD [10], they show that a branch at carbon 2 in chain significantly increases viscosity. The calculations of kinematic viscosity of three C 30 H 62 isomers which are lubricant basestocks are done in [11]. The detailed comparison of the classical force fields is made in [12] in terms of prediction of lubricants viscosity. The pressure dependence of viscosity for squalane is predicted in the work by Jadhao and Robbins [13] which references the extensive experimental work [14]. However, they have no simulations in the region of the super-Arrhenius pressure dependence of viscosity which is important to friction [3]. Also, the predictive power of the different force fields is not studied enough at the high pressures.
The choice of atomic interaction model determines the accuracy of properties prediction and computational cost. For the theoretical studies, authors usually use simple models such as point particle models [15][16][17][18], corse-grained [19] and united atom models [20] that give qualitative results to verify or expand theory. If ones need accurate estimation of properties, more complex models should be used [21][22][23][24][25][26]. The 2,2,4-trimethylpentane liquid is a good choice to verify different force fields since a lot of experimental data are available. The dependence of 2,2,4trimethylpentane viscosity on pressure is studied experimentally by Dymond in 1985 [27] (for pressures up to 500 MPa) and Padua in 1996 [28] (for pressures up to 100 MPa).
Here, the molecular dynamics methods are used to compute the shear viscosity coefficient for 2,2,4-trimethylpentane at 298 K and two pressures (0.1 MPa and 400 MPa). Three force fields are compared by prediction of the viscosity. The rest of the paper is arranged as follows. In Section 2, details about force fields, molecular dynamics details, system equilibration and Green-Kubo method are presented. The calculated densities and viscosities at 0.1 MPa and 400 MPa for 2,2,4-trimethylpentane and their comparison with the experimental data are discussed in Section 3.

Simulation details 2.1. Force fields for branched alkane
The potential energy of the 2,2,4-trimethylpentane liquid can be modeled as follows: where the first three terms are valence bond, valence angle and dihedral interactions correspondingly. Non-bonded forces are described by Lennard-Jones and Coulomb terms. Two all-atom class I force fields are used: OPLS-AA [8] and its modified version for long hydrocarbons L-OPLS-AA [29]. The bond and angle interactions are described by harmonic oscillators, the dihedral terms are presented via Fourier function. The Lennard-Jones interactions are treated via classical 6-12 term.
These two force fields are compared in terms of prediction of the density and viscosity with the COMPASS force field [30]. In the COMPASS potential, the bond and angle interactions are treated as non-harmonic oscillators. The dihedral interactions are described by the similar functions as in the OPLS-AA and L-OPLS-AA force fields. The main difference is that all the cross terms between bonds, angles, dihedrals are included. The Lennard-Jones interactions are presented via 6-9 form which, as expected, give more trustworthy results at the high pressures. All the interaction parameters are taken from the LAMMPS user forum [31].

Molecular dynamics details
All the calculations are performed for the equilibrated systems of the 1000 2,2,4-trimethylpentane molecules. The non-bonded 1-4 interactions are scaled with 0.5 factor in OPLS-AA and L-OPLS-AA. The cutoff radius for non-bonded forces is 12Å, long-range Coulumb interactions are calculated using the particle-particle particle-mesh method [32] with the desired relative error in forces 10 −5 . The tail corrections for the pressure are also applied which are analitically calculated long-range Van der Waals contributions to the energy and pressure [30]. The equations of motion are integrated using rRESPA algorithm [33] with a timestep of 1 fs for Lennard-Jones and Coulomb interactions. The bonds and angles are treated with 0.125 fs, the dihedrals and impropers with 0.5 fs timestep. The simulations are performed in periodic boundary conditions (PBC) using the LAMMPS package [34] with the use of GPU package [35][36][37].

System equilibration
The molecular typologies for all the force fields are generated using msi2lmp tool in LAMMPS. Each equilibrated system is obtained in the following way. The initial configuration consists of randomly packed 1000 molecules at density equal to 0.7 g/cm 3 . At the first stage, the 3 minimization of the system energy is carried out and system relaxes in the Langevin thermostat for 50 ps at T = 298 K. At the second stage, the NPT simulation for 1 ns is performed for the calculation of the equilibrium density at the given pressure and T = 298 K. After that, the system is equilibrated in the NVT ensemble for 1 ns at the corresponding equilibrium density and T = 298 K.

Green-Kubo method for viscosity
The Green-Kubo formula for the shear viscosity η αβ in αβ-plane is [38]: where C σ (t) is an autocorrelation function, σ αβ are off-diagonal components of the stress tensor, V and T are system volume and temperature, k B is Boltzmann's constant. In practice, the integral in (2) is typically calculated up to some time t when C σ decays to zero to the accuracy of numerical simulation. . . . in (3) is an average over the canonical ensemble. The shear viscosity η is found as an average of the η xy , η xz and η yz . The stress tensor σ αβ is calculated from the following equation: where N is a number of atoms, r iα and v iα are α-components of coordinate and velocity of the i-th atom, and f iα is α-component of the force that acts on the i-th atom. The Green-Kubo integral converges without any problems in the cases of atomic or simple molecular liquids [39][40][41][42]. In liquids with high viscosities, it becomes a challenge to achieve convergence of the Green-Kubo integral (2) because of high correlation times of C σ (t). Zhang and Maginn proposed the time decomposition method (TDM) [43] which allows to take Green-Kubo integral accurately in the case of ionic liquids and successfully applied it for hydrocarbons [44]. Here, we also use TDM for the viscosity calculation of 2,2,4-trimethylpentane.
In the TDM method, the Green-Kubo integral is fitted by the double exponential function where A, α, τ 1 and τ 2 are the fitting parameters. The C σ (t) long time numerical errors are taken into account with 1/σ 0.5 weight in the fitting procedure.

Density prediction
The equilibrium densities of 2,2,4-trimethylpentane obtained from the NPT simulations at 298 K for two pressures (0.1 MPa and 400 MPa) are presented in table 1. At the normal pressure, OPLS-AA overestimates density. The modified force field L-OPLS-AA gives the value that lies closer to the experimental data. It matches the expectations because L-OPLS-AA is fitted to reproduce the equation of state of linear alkanes. COMPASS prediction (obtained in [45]) is the best in this case. At the higher pressure, the L-OPLS-AA potential represents more accurate value then OPLS-AA and COMPASS.

Viscosity prediction
The viscosity coefficients are calculated using the Green-Kubo relation. Since the viscosity does not depend on the size effects [44,46] that take place in the diffusivity case [44,[47][48][49], only 1000 molecules in the simulation box are considered. The convergence of the Green-Kubo integral is obtained by the time decomposition method [43]. For the normal pressure, ten trajectories of 1 ns length are enough for TDM to achieve the accurate dependence of the Green-Kubo integral on upper time limit. For 400 MPa, twenty trajectories are required due to the higher relaxation times of the Green-Kubo integral at this pressure. The TDM fitting procedure is made in Gnuplot. The calculated values of A, α, τ 1 and τ 2 are shown in table 2.   Since COMPASS requires almost the same number of the computational resources as OPLS-AA and L-OPLS-AA, it has much better predictive power in terms of the shear viscosity calculations.

Conclusions
The molecular dynamics calculations of the shear viscosity coefficient are carried out for 2,2,4trimethylpentane at 293 K and two pressures (0.1 MPa and 400 MPa). The Green-Kubo method is used for the shear viscosity computation. The results obtained in the OPLS-AA, L-OPLS-AA and COMPASS force fields are compared. The simulations are performed in LAMMPS package.
Both OPLS-AA and L-OPLS-AA overestimate the viscosity at the pressures under study. As pressure increases, the OPLS-AA overestimation becomes higher. Thus, the viscosity is 25.5 ± 2 cP in OPLS-AA and the experimental value is 9.5 ± 0.7 cP. The L-OPLS-AA value (17.2 ± 1.4 cP) becomes closer to the experimental viscosity at 400 MPa.
COMPASS gives accurate predictions of the shear viscosity at both pressures. At 0.1 MPa, the calculated viscosity matches the experimental data within the standard error. Since it requires almost the same number of computational resources as OPLS-AA and L-OPLS-AA, it has much better predictive power in terms of the shear viscosity calculations.