One of the key points to better understand the origins of friction is to know how two surfaces in contact adhere to one another. In this paper we present molecular dynamics (MD) simulations of two aluminium bodies in contact, exposed to a range of normal loads. The contact surfaces of both aluminium bodies have a self-affine fractal roughness, but the exact roughness varies from simulation to simulation. Both bodies are allowed to have an adhesive interaction and are fully deformable. Tracking important contact parameters (such as contact area, number of contact clusters, and contact pressure) during a simulation is challenging. We propose an algorithm (embedded within a parallel MD code) which is capable of accessing these contact statistics. As expected, our results show that contact area is increasing in proportion with applied load, and that a higher roughness reduces contact area. Contact pressure distributions are compared to theoretical models, and we show that they are shifted into the tensile regime due to the inclusion of adhesion in our model.