Modelling metalloproteins using the classical force fields is challenging. Several methods have been devised to model metalloproteins in force fields. Of these methods, the bonded model, combined with Restrained Electrostatic Potential (RESP) charge fitting, proved its superiority. The latter method was facilitated by the development of the python-based Metal Centre Parameter Builder (MCPB.py) AmberTool. However, the standard bonded model method offered by the MCPB.py tool may not be appropriate for validating and refining the binding modes predicted by docking when crystal structures are lacking. That is because the representation of coordination interactions between any bound ligand and metal ions by covalent bonds can hinder the flexibility of the ligand. Therefore, a new modification to the standard bonded model approach is proposed here. Molecular dynamics (MD) simulations based on the new modified bonded model (MBM) approach avoid the bias caused by coordination bonds and, unlike hybrid QM/MM MD, allow for sufficient sampling of the binding mode given the currently available computational power. The MBM MD approach reproduced the studied crystal structure conformations of New Delhi Metallo-β-lactamase 1 (NDM-1). Furthermore, the MBM approach described the binding interactions of intact β-lactams with NDM-1 reasonably, and predicted a non--productive binding mode for the poor NDM-1 substrate aztreonam whilst predicting productive binding modes for known good substrates. This study presents a useful MD method for metallo-β-lactamases and provides better understanding of β-lactam substrates recognition by NDM-1. The proposed MBM approach might also be useful in the investigation of other metal-containing protein targets.