To calculate the bulk thermal conductivity, κ ∞ , the inverse of the thermal conductivity, 1 κ , was plotted against the inverse of the system length, 1 lz . The results were fitted linearly and κ ∞ was given by the inverse of the y intercept. The thermal conductivity was also calculated in the [001] and [010] directions to investigate the anisotropy of this property in hematite. Figure in the Supplementary Information illustrates a typical temperature profile along the z-axis for the Fe 2 O 3 systems during the NEMD simulations. The system was divided into 20n bins, where n was the number of times the initial 50 Å simulation was replicated in the z axis. The points on the graph represent the temperature of each bin, averaged over the final 5 ns of the simulation. As the box length, l z , increased, so did the size of the non-linear regions by the thermostats and as such, an exclusion zone of 5n Å was left on either side of the thermostats when calculating the temperature gradient to avoid fitting in the nonlinear region. The heat flux was obtained by calculating the rate of heat exchange, ⟨ Ė⟩, at the hot and cold thermostats and the thermal conductivity was calculated using Fourier's law. In all cases investigated here, there was a < 1 % difference in the total energy exchanged at the hot and cold thermostats during NEMD simulations, indicating good energy conservation. The thermal conductivity was calculated for five systems ranging from l z = 50 to 250 Å. For all four force fields, κ was found to increase as a function of l z . To obtain the bulk thermal conductivity, κ ∞ , 1 κ was plotted against 1 lz , as displayed for each force field in Figure in the Supplementary Information. For ClayFF, INTERFACE-FF and the Morse potential, all system sizes were found to fall within the linear regime; however, for the harmonic model only the systems between 100 and 250 Å fell in the linear regime. The linear regime was fitted and extrapolated to obtain the y intercept, which was equal to κ∞ .