AbstractAntibiotic resistance is a major threat to global public health. β‐lactamases, which catalyze breakdown of β‐lactam antibiotics, are a principal cause. Metallo β‐lactamases (MBLs) represent a particular challenge because they hydrolyze almost all β‐lactams and to date no MBL inhibitor has been approved for clinical use. Molecular simulations can aid drug discovery, for example, predicting inhibitor complexes, but empirical molecular mechanics (MM) methods often perform poorly for metalloproteins. Here we present a multiscale approach to model thiol inhibitor binding to IMP‐1, a clinically important MBL containing two catalytic zinc ions, and predict the binding mode of a 2‐mercaptomethyl thiazolidine (MMTZ) inhibitor. Inhibitors were first docked into the IMP‐1 active site, testing different docking programs and scoring functions on multiple crystal structures. Complexes were then subjected to molecular dynamics (MD) simulations and subsequently refined through QM/MM optimization with a density functional theory (DFT) method, B3LYP/6‐31G(d), increasing the accuracy of the method with successive steps. This workflow was tested on two IMP‐1:MMTZ complexes, for which it reproduced crystallographically observed binding, and applied to predict the binding mode of a third MMTZ inhibitor for which a complex structure was crystallographically intractable. We also tested a 12‐6‐4 nonbonded interaction model in MD simulations and optimization with a SCC‐DFTB QM/MM approach. The results show the limitations of empirical models for treating these systems and indicate the need for higher level calculations, for example, DFT/MM, for reliable structural predictions. This study demonstrates a reliable computational pipeline that can be applied to inhibitor design for MBLs and other zinc‐metalloenzyme systems.