To initialise the simulations, fluid particles were packed with random positions and orientations into a large box (periodic in x, y, and z, with x and y dimensions 10× the target values and z dimension equal to the target value). For systems to be simulated with the graphene substrate, the target x and y dimensions were 34.48σ and 34.125σ , respectively, to enable periodic packing of the substrate particles; for the silicon substrate, the final dimensions were x = y = 33.94σ for the same reason. The size in the z dimension was used to control the volume fraction of fluid particles and so depended on fluid particle size; the values for various systems can be found in the ESI Table . A soft potential U soft (r) = A 1 + cos πr r c for r < r c , where r c is is a cutoff distance and A an energy pre-factor that varies the "hardness" of the potential, was applied to remove particle overlaps (r c = 5σ , A increased linearly from 0 to 16.7ε over 20,000 timesteps). Interactions were changed from the soft potential to the GB potential, and the box dimensions shrunk linearly over 50,000 timesteps at a higher temperature of T = 0.88ε/k B to give the final dimensions outlined above. To initialise the system in an isotropic phase, the fluid particles were then allowed to equilibrate in the bulk system at T = 0.88ε/k B for another 50,000 time steps. The periodic boundary conditions in the z direction were then removed, and the box extended in the z direction to give a region of vacuum above the bulk fluid. The solid substrate was introduced below the fluid, and the harmonic wall at the vapour interface was positioned at the z coordinate required to give the desired overall volume fraction (ESI Table ). The energy of the system was then minimised to remove any overlaps between fluid and interfaces. The system was simulated for a further 510,000 timesteps at the specified temperature, with the equilibration and correlation time for the orientation of the top and bottom fluid layers determined using pymbar's timeseries module, which estimates the equilibration time as the time that maximises the number of uncorrelated samples (see ref. 71 for details). Only the data after the equilibration time was used for further analysis. Although a number of systems were simulated that were anisotropic in the bulk (typically for the most anisotropic shapes and interactions; see ESI Table for a list of systems that formed anisotropic bulk phases), the slow dynamics in these systems did not allow equilibrium properties to be measured reliably over the time scales that were simulated. Results are therefore only reported in the main paper for systems with an isotropic bulk fluid phase.