Molecular Mechanics Poisson-Boltzmann Surface Area (MMPBSA) methods have become widely adopted in estimating protein-ligand binding affinities due to their efficiency and high correlation with experiment. Here different computational alternatives were investigated to assess their impact to the agreement of MMPBSA calculations with experiment. Seven receptor families with both high-quality crystal structures and binding affinities were selected. First the performance of nonpolar solvation models was studied and it was found that the modern approach that separately models hydrophobic and dispersion interactions dramatically reduces RMSD's of computed relative binding affinities. The numerical setup of the Poisson-Boltzmann methods was analyzed next. The data shows that the impact of grid spacing to the quality of MMPBSA calculations is small: the numerical error at the grid spacing of 0.5 Å is already small enough to be negligible. The impact of different atomic radius sets and different molecular surface definitions was further analyzed and weak influences were found on the agreement with experiment. The influence of solute dielectric constant was also analyzed: a higher dielectric constant generally improves the overall agreement with experiment, especially for highly charged binding pockets. The data also showed that the converged simulations caused slight reduction in the agreement with experiment. Finally the direction of estimating absolute binding free energies was briefly explored. Upon correction of the binding-induced rearrangement free energy and the binding entropy lost, the errors in absolute binding affinities were also reduced dramatically when the modern nonpolar solvent model was used, although further developments were apparently necessary to further improve the MMPBSA methods. © 2016 Wiley Periodicals, Inc.