Hexagonal boron nitride sp2layers grown and supported on the Rh(111) metal surface attractedquite some interest thanks to the structura l and electronic peculiarities of this quasi-2D system . Thehoneycomb regular corrugation is the key feature at the origin of several properties and applications innanotechnology, e.g., the selective adsorption and functionalisation related to the modulation of the elec-tronic structure. Atomistic simulation s play an important role, since they can shed light on the natu re ofsuch a complex interface, providing resolutio n of details that cannot be achieved experimentally. However,the studies by electronic structure calculations have been mostly limited to static models of the optimi zedsystem. The sampling of conﬁgurations at ﬁnite temperature by ab-initio molecular dynamics requires sig-niﬁcantly larger computational eﬀort, and can become unfea sible for large scale and metallic models, as itis the case of h-BN/Rh(111). In this work, we employ a recently developed Car-Parrinello-like approach toovercome the performance limitations of the standard Born-Oppenheimer molecular dynamics scheme, thusobtaining a speed-up of 17×. We report on the set-up and the application of this approach to simulate the h-BN/Rh(111) interfac e at diﬀerent temperatures and disc uss the thermal stability of the corru gated pattern.