Throughout this section, position restraints of 100 000 kJ/mol/nm 2 were applied on all protein heavy atoms to bias results to the reference 1XLW geometry. Systems were first subjected to up to 1000 steps of conjugate gradient minimization, followed by 1 ps QM/MM equilibration with 1 fs time step. Finally, production runs were performed in 3 replicas of 14 parallel walkers each. A time step of 0.5 fs was used, and each walker sampled for 50 ps. Velocity-rescale thermostat was used with no pressure control. Multiple-walkers well-tempered metadynamics was used to sample proton location possibilities and reconstruct free energy profiles. Collective variables used are illustrated in Figures and, and more detail is given in the Supplementary Information section. Metadynamics potential 1.046 kJ/mol high was applied each 100th step. Potential width was 0.25 radians for water molecule rotation and 0.2Å for all distance-based variables. Bias factor was set to 20 except for runs only simulating water molecule rotation, in which it was 12. Free energy profiles were reconstructed with Tiwary reweighting, with the first 5 ps of each walker discarded as initial transient. 15 systems corresponding to deepest local minima were selected for subsequent analysis. For that, positional restraints were lifted from the QM region, and each system was sampled for 10 ps with 0.5 fs step in 50 replicas. dRMSD values were assessed by aggregating the last 1 ps of each replica.