Reliable Viscosity Calculation from High-Pressure Equilibrium Molecular Dynamics: Case Study of 2,2,4-Trimethylhexane

Abstract

Viscosity is a fundamental property of liquid lubricants, yet it is challenging to determine accurately, especially at high pressures. Although equilibrium molecular dynamics (EMD) simulations are a promising alternative to resource-intensive experiments, practical challenges remain in assessing the sufficiency of simulation time and in controlling uncertainties in the Green-Kubo formalism due to the finite amount of trajectory data. In this work, we extend the STable AutoCorrelation Integral Estimator (STACIE), a recently developed algorithm for estimating transport properties. First, we introduce the Lorentz model to estimate the viscosity and the exponential correlation time from the low-frequency power spectrum of anisotropic pressure fluctuations. Second, we show how to supplement the three conventional off-diagonal elements of the pressure tensor ($P_{xy}$, $P_{yz}$ and $P_{zx}$) with two additional independent time series for shear viscosity calculations. Using these improvements, we apply STACIE to calculate the shear viscosity of 2,2,4-trimethylhexane from EMD simulations. We demonstrate STACIE's capability to reliably calculate viscosity under high-pressure conditions, offering a robust and automated solution with validated uncertainty quantification. Our results, when compared to the outcomes of the 10th International Fluid Properties Simulation Challenge (IFPSC), underscore the need for long EMD simulations. Large deviations from experimental viscosities in previous works were primarily due to insufficient simulation times and ad hoc post-processing choices, rather than the limitations of the force fields used. Unlike previous studies, our viscosity estimates agree well with experimental results (relative error < 6%) up to the highest pressure of 1 GPa, highlighting the improved reliability and accuracy of STACIE's systematic approach to viscosity predictions.

Figures (5)

• Figure 1: Preparation of five independent anisotropic pressure contributions per trajectory as input for STACIE's shear viscosity calculation. (a) Illustration of pressure tensor elements for a cubic simulation box containing 100 2,2,4-trimethylhexane molecules. (b) Block averages of pressure tensor elements printed using the "fix ave/time" command in LAMMPS. (c) Transformation of the pressure tensor into the five independent terms required as input for STACIE, resulting in $M=5 N_\text{traj}$ input time series in total.
• Figure 2: STACIE viscosity analysis of 2,2,4-trimethylhexane at ambient conditions (0.1MPa and 293K), with 50 trajectories, each 2ns long. (a) The final part of STACIE's screen output showing the shear viscosity estimate and recommended simulation settings. Note that the "time step" is not the MD integration time step, but the block size in units of time, as this is the input time series given to STACIE. (b) PSD averaged over the five independent pressure tensor elements and 50 MD trajectories, fitted with the Lorentz model. A detailed description of all elements of the spectrum plot can be found in Section S3 of the Supporting Information.
• Figure 3: Comparison of shear viscosity results with experimental data Bair2019Pressure and entries from the 10$^\text{th}$ IFPSC Bair2019PressureZheng2019Gong2019Kondratyuk2019Messerly2019Cunha2019. (a) Shear viscosity as a function of pressure up to 1GPa. (b) Relative error in viscosity versus pressure. The dashed green line indicates the experimental value as the zero-reference Bair2019Pressure. A "symlog" scale with a linear threshold of 100 is used on the y-axis to emphasize discrepancies. (c) EMD simulation times as a function of pressure, with the solid blue line representing the minimum required simulation time based on Eq. \ref{['eq:min-simtime']}. The legend in (a) applies to all subplots, with consistent gray shades used for each IFPSC entry.
• Figure 4: Shear viscosity with its standard error estimated from EMD trajectories as a function of simulation time. (a) $\bar{P}_\text{MD} = 494MPa$ with 60ns-long trajectories. (b) $\bar{P}_\text{MD} = 1017MPa$ with 500ns-long trajectories. The green line and dashed lines represent the experimental value and its uncertainty Bair2019Pressure, respectively, extrapolated to the average pressure of our production runs. The orange data points correspond to viscosity estimates of Kondratyuk et al.Kondratyuk2019.
• Figure 5: Simulation times as a function of viscosity. Blue points show the simulation times used in this study to compute viscosities up to 1GPa. STACIE's recommended simulation times based on the exponential correlation time, $\tau_\text{exp}$, are shown in orange. Red points represent the extrapolated viscosities and corresponding simulation times at 1.5GPa and 2GPa, assuming an exponential viscosity increase with pressure. The fitted model is shown as a dashed black line.