Initially, a 50.3 × 43.6 × 54.9 Å3 hematite surface was generated, which was then replicated in the z direction, parallel to the [100] crystallographic direction, to create additional systems of length 109.8, 164.8, 219.7 and 274.6 Å. The systems were energy minimised, followed by a 1 ns equilibration in the NPT ensemble at 300 K with an anisotropic pressure coupling of 1 atm in the z-axis, using a Nosé-Hoover thermostat and barostat with temperature and pressure damping coefficients of 100 fs and 1,000 fs respectively. Finally, the system underwent a 0.25 ns equilibration in the NVT ensemble at 300K using a Nosé-Hoover thermostat. Once the system had equilibrated, NEMD simulations were performed to calculate the thermal conductivity by applying a heat flux along the z-axis. 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 a Langevin thermostat. Previous comparisons suggested that the Langevin thermostat is preferrable to the Nosé-Hoover thermostat for this purpose. NEMD simulations were performed for 3 ns, with the initial 0.5 ns being discarded to allow the system to reach a non-equilibrium steady state. The temperature gradient, ∇T , and heat flux, J q , were calculated using method outlined previously, allowing the size dependant thermal conductivity to be calculated.