MD simulations were performed using the Large-scale Atomic/Molecular Massively Parallel Simulator (LAMMPS) software. The velocity-Verlet algorithm 69 was used to integrate the MD equations of motion. A time step of 1 fs was used for all of the simulations and all bonds involving H atoms were fixed using the SHAKE algorithm. Long-range electrostatic interactions were calculated using the slab implementation of the particle-particle, particle- mesh (PPPM) algorithm with a relative force accuracy of 10 -4 . The PMFs were extracted using the ABF method with the collective variables (colvars) module implemented in LAMMPS. Before starting the PMF production runs, the systems were energy minimized using the conjugate gradient algorithm, with all of the hematite atoms being fixed in place. Subsequently, the system was equilibrated at 300 K and 101 kPa for 2 ns. The outer layer of atoms in the bottom wall of hematite were fixed in all directions, whereas the outer layer of the top wall was fixed only in the x-and y-directions, but allowed to move in the z-direction during the equilibration phase. These outer layers consisted of all of the Fe and O atoms that were 1 Å or less from the bottom and the top of the simulation box respectively (440 atoms in both cases).