Many machine components, such as bearings and gears, operate under the so-called Thermo-Elasto-Hydrodynamic (TEHL) lubrication regime. In this specific regime, the surfaces in relative motion are entirely separated by a thin lubricant film ( < 1 μm) which can experience ultra-high pressures, high shear rates, and high temperatures. These state variables heavily influence the thermomechanical properties of the lubricants due to intensified interaction between molecules. As it is very difficult to perform experiments under such extreme conditions, Molecular Dynamics (MD) simulations offer a promising alternative for extracting thermomechanical properties of lubricants, especially the viscosity, which is one of the most crucial properties for lubricants under TEHL operating conditions. There are two main types of MD simulations that are used to describe the dependence of thermomechanical properties on the state variables (pressure, temperature, and shear rate), namely the Equilibrium MD (EMD) and Non-Equilibrium MD (NEMD). While the EMD methods are typically used to model the Newtonian behavior of fluids in the limit of vanishing shear rates, the NEMD methods are used to model the non-Newtonian behavior at higher shear rates. In MD simulations, the accuracy is mainly determined by the accuracy of force fields used for describing the inter-atomic interactions in the system. When EMD simulations are used, apart from the accuracy of force fields, the specific way of postprocessing the EMD trajectories has a serious impact on the accuracy of the calculated properties. In this study, the EMD simulations are used to determine the viscosity of lubricants under wide pressure and temperature range. A reliable postprocessing procedure was developed to analyze EMD trajectories for viscosity calculations with increased accuracy and robustness.