Periodic boundary conditions were applied in all systems. Electrostatic interactions were modeled with the particle mesh Ewald method, while van der Waals interactions were modeled with Lennard-Jones potentials which were smoothly switched to zero in the range The leap-frog integrator with an integration step of 2 fs was used. Bond lengths to hydrogens in the solute were constrained with the LINCS algorithm, while the internal degrees of the CHARMM TIP3P water molecules were constrained with the SETTLE algorithm. Prior to the constant pH MD simulations, the energy of all systems was minimized using the steepest descent method, followed by a 1 ns equilibration.