The system was equilibrated in the NPT ensemble for a total of 100 ns, using the Parrinello-Rahman barostat and a stochastic velocity rescaling thermostat 41 at 300 K with a 1 ps coupling time. Hydrogen-containing bonds were constrained using the LINCS algorithm and the SETTLE algorithm for water molecules, which allowed the use of a 2 fs propagation time step. Particle mesh Ewald 44 with a cutoff distance of 12 Å was used to treat long-range electrostatic interactions, and the same cutoff was used for Van der Waals parameters. Position restraints on heavy atoms were gradually decreasing until complete removal (Table for the detailed procedure), while the restraints on the Mg 2+ coordination sphere were kept during the entire equilibration. From the end of this equilibration protocol, in order to test the performance of different force fields to describe the active site structure, an additional 20 ns NPT equilibration was performed with each of the tested force fields, followed by a 500 ns equilibration keeping only restraints (k = 20000 kcal/mol/nm 2 ) on the distance between Mg 2+ ions and the non-water oxygen atoms in their hydration sphere. These restraints were then removed and an additional 500 ns NPT production run was performed with each force field in 4 replicates, with different initial velocities.