One of the most common computational binding affinity prediction techniques is molecular mechanics Poisson-Boltzmann surface area (MMPBSA) . This is an approximate post-processing end-state method, which uses continuum solvent models to reduce the computational cost of obtaining results. The speed and ease of setup (compared to rigorous free energy calculations) makes MMPBSA an attractive candidate for use throughout the drug discovery pipeline. However, results are often seen to be system dependent and are widely perceived to be less accurate than those obtained from more expensive and theoretically rigorous approaches (such as free energy perturbation, FEP, and thermodynamic integration, TI) . Furthermore, the term MMPBSA as used in the literature permits a wide range of variants which incorporate different sampling strategies (for example, all ligand conformers can be drawn from simulation of the complex or from independent runs) and differing solvation and entropy terms. Our previous work has demonstrated that MMPBSA analysis of single simulations is highly unreliable with calculations initiated from the same structures varying by up to 12 kcal mol -1 for small molecules bound to HIV-1 protease and even more for flexible ligands binding to MHC . This served as the inspiration for our ESMACS (enhanced sampling of molecular dynamics with approximation of continuum solvent) protocols which use ensemble simulations that have been shown to produce results with reproducible uncertainties of less than 2 kcal mol -1 for a range of systems . In this work we seek to assess the performance of the approach in a challenging dataset containing a highly varied set of ligands which interact with water in the protein binding site. We assess the impact on protocol performance of multiple trajectory sampling, ligand parameterization, inclusion of explicit water molecules and a recently developed approach to calculating the entropic contribution to the binding free energy.