Various approaches have been proposed to include the effect of pH in Molecular Dynamics (MD) simulations. Among these, the lambda-dynamics approach proposed by Brooks and co-workers can be performed with little computational overhead and hence be used to routinely perform MD simulations at microsecond timescales, as shown in the accompanying paper. At such timescales, however, the accuracy of the molecular mechanics force field and the parameterization becomes critical. Here, we address these issues and provide the community with guidelines on how to set up and perform long timescale constant pH MD simulations. We found that barriers associated with the torsions of side chains in the CHARMM36m force field are too high for reaching convergence in constant pH MD simulations on microsecond timescales. To avoid the high computational cost of extending the sampling, we propose small modifications to the force field to selectively reduce the torsional barriers. We demonstrate that with such modifications we obtain converged distributions of both protonation and torsional degrees of freedom and hence consistent pKa estimates, while the sampling of the overall configurational space accessible to proteins is unaffected as compared to normal MD simulations. We also show that the results of constant pH MD depend on the accuracy of the correction potentials. While these potentials are typically obtained by fitting a low-order polynomial to calculated free energy profiles, we find that higher-order fits are essential to provide accurate and consistent results. By resolving problems in accuracy and sampling, the work described in this and the accompanying paper paves the way to the widespread application of constant pH MD.