Methyl parathion hydrolase (MPH) is an enzyme of the metallo-β-lactamase superfamily, which hydrolyses a wide range of organophosphates (OP). Recently, MPH has attracted attention as a promising enzymatic bioremediator. The crystal structure of MPH enzyme shows a dimeric form, with each subunit containing a binuclear metal ion center. MPH also demonstrates metal ion-dependent selectivity patterns. The origins of these patterns remain unclear but are linked to open questions about the more general role of metal ions in functional evolution and divergence within enzyme superfamilies. We aimed to investigate and compare the binding of different OP pesticides to MPH with cobalt(II) metal ions. In this study, MPH was modelled from Ochrobactrum sp. with two different classes of OP pesticides bound, including phosphomonoester (methyl paraoxon and dichlorvos) and S-substituted thiophsphotriester (profenofos). The docked structures for each substrate optimized by DFT calculation was selected and subjected to atomistic molecular dynamics simulations for 500 ns. It was found that alpha metal ions did not coordinate with all the pesticides. Rather, the pesticides coordinated with less buried beta metal ions. It was also observed that the coordination of beta metal ions was perturbed to accommodate the bulky pesticides. The binding free energy calculations and structure-based pharmacophore model revealed that all the three substrates could bind well at the active site. However, profenofos exhibits the stronger binding affinity to MPH in comparison to other two substrates. Therefore, the ability of the in silico analysis presented here could be informative for increasing enzyme stability and activity.