In this study, interactions of the catalytically active binuclear form of glycerophosphodiesterase (GpdQ) with four chemically diverse substrates, i.e. *** (a phosphomonoester), BNPP and GPE (both phosphodiesters), and paraoxon (a phosphotriester) have been investigated using all-atom molecular dynamics (MD) simulations. The roles of metal ions and key amino acid residues, coordination flexibility, and dynamic transformations in all enzyme–substrate complexes have been elucidated. The roles of important first and second coordination shell residues in substrate binding and coordination flexibility of the enzyme suggested by simulations are supported by experimental data. The chemical nature of the substrate is found to influence the mode of binding, electrostatic surface potential, metal–metal distance, and reorganization of the active site. The experimentally proposed association between the substrate binding and coordination flexibility is analyzed using principal component analysis (PCA), movements of loops, and root-mean-square-fluctuations (RMSF) as parameters. The PCA of these substrates provides different energy basins, i.e. one, three, two and five for ***, BNPP, GPE, and paraoxon, respectively. Additionally, the area of an irregular hexagon (268.3, 288.9, 350.8, and 362.5 Ų) formed by the residues on these loops illustrates their distinct motions. The substrate binding free energies of ***, BNPP, and GPE are quite close (22.4–24.3 kcal mol⁻¹), but paraoxon interacts with the smallest binding free energy (14.1 kcal mol⁻¹). The metal binding energies in the presence of these substrates are substantially different, i.e. the lowest for *** and the highest for paraoxon. These results thus provide deeper insight into the chemical promiscuity and coordination flexibility of this important enzyme.