The use of conformational ensembles provided by nuclear magnetic resonance (NMR) experiments or generated by molecular dynamics (MD) simulations has been regarded as a useful approach to account for protein motions in the context of pK_a calculations, yet the idea has been tested occasionally. This is the first report of systematic comparison of pK_a estimates computed from long multiple MD simulations and NMR ensembles. As model systems, a synthetic leucine zipper, the naturally occurring coiled coil GCN4, and barnase were used. A variety of conformational averaging and titration curve-averaging techniques, or combination thereof, was adopted and/or modified to investigate the effect of extensive global conformational sampling on the accuracy of pK_a calculations. Clustering of coordinates is proposed as an approach to reduce the vast diversity of MD ensembles to a few structures representative of the average electrostatic properties of the system in solution. Remarkable improvement of the accuracy of pK_a predictions was achieved by the use of multiple MD simulations. By using multiple trajectories the absolute error in pK_a predictions for the model leucine zipper was reduced to as low as approximately 0.25 pK_a units. The validity, advantages, and limitations of explicit conformational sampling by MD, compared with the use of an average structure and a high internal protein dielectric value as means to improve the accuracy of pK_a calculations, are discussed.