The fluid regions of our systems were generated using the Materials and Process Simulations platform (MAPS) from Scienomics SARL (Paris, France). In order to identify a valid model for the fluid, a series of benchmark simulations were performed using different force fields, with each being assessed on its ability to predict the density, thermal conductivity, and viscosity of the fluid. Four force fields were selected for this investigation, three atom-istic (LOPLS-AA, CHARMM36, AMBER 54 ) and one united-atom (TraPPE-UA). MD simulations were performed in LAMMPS 56 using a velocity-Verlet integrator. For the atomistic force fields, a time step of 0.5 fs was used, while for the united-atom force field this was increased to 5 fs. Note that satisfactory energy conservation was not observed for the atomistic force fields with a timestep of 1 fs. We did not constrain the C-H bonds with the SHAKE algorithm due to the difficulties in assigning the local temperature with rigid constraints. We used a Lennard-Jones (LJ) cut-off of 12 Å for all of the force fields. The unlike LJ interactions were calculated using arithmetic mean mixing rules for CHARMM36 and AMBER, while the geometric mean (Lorentz-Berthelot) mixing rules were used for LOPLS-AA and TraPPE-UA. The long-range Coulomb interactions between the partial charges on the carbon and hydrogen atoms in the atomistic force fields were calculated using the particle-particle particle-mesh (PPPM) algorithm with a relative error in the forces of 1 To calculate the thermal conductivity, a box was created with periodic boundary conditions in all three directions, containing 400 PAO4 molecules at an initial density of 0.8 g cm -1 . The systems were energy minimised, followed by a 1 ns equilibration in the isobaricisothermal (NPT) ensemble at 300 K and a pressure of 1 atm. We used the Nosé-Hoover thermostat and barostat with temperature and pressure damping coefficients of 0.1 ps and 1.0 ps, respectively. Next, the systems underwent a 0.25 ns equilibration in the canonical (NVT) ensemble at 300 K using a Nosé-Hoover thermostat. NEMD simulations were then performed to calculate the thermal conductivity by applying a heat flux along the z axis using a boundary driven approach. Hot and cold thermostat regions were defined in the xy plane, with a thickness of 5 Å, at the centre and edges of the simulation box. Atoms falling within these regions were thermostatted to 325 K and 275 K for the hot and cold regions respectively, using the Langevin thermostat. A temperature damping coefficient of 1.0 ps was used, which ensured good energy conservation during then NEMD simulations. The NEMD simulations were performed for 5.5 ns, with the initial 0.5 ns discarded to allow the system to reach a non-equilibrium steady state.