With a higher accuracy of the underlying deprotonation free energy profiles, we found that the correct sampling of protonation states also depends critically on the order of the polynomial fit used to obtain an analytical form for these correction potentials. We show that the commonly accepted first-order fit, 12 although firmly based on linear response theory, is not sufficiently accurate and can lead to erroneous protonation dynamics in constant pH MD simulations.