As anticipated, structures using covalent metal models remained stable over the simulation time (Figure ; see Figure -S16 for RMSD, coordination sphere analysis and snapshots of individual trajectories), with a median RMSD = 1.9 Å relative to the crystal structure. In contrast, most of the simulations utilizing non-bonded and dummy models resulted in disassembly. Only ZIF-8 and ZIF-67 MOFs produced stable simulations with these models, possibly due to the absence of competing interactions with counter ions and solvent. Removing counterions in the other systems resulted in a marginal improvement in stability. While the stability of the covalent model is expected as dissociation is not allowed, the poor performance of other models for almost all systems was surprising. This is likely the result of parameters being overfitted to reproduce aqueous complexes rather than interactions with other heteroatoms containing molecules, leading to imbalanced metal-ligand and metal-counterions interactions.