A high-throughput virtual screening pipeline has been extended from single energetically minimized structure Molecular Mechanics/Generalized Born Surface Area (MM/GBSA) rescoring to ensemble-average MM/GBSA rescoring. For validation, the binding affinities of a series of antithrombin ligands have been calculated by using the two MM/GBSA rescoring methods. The correlation coefficient (R2) of calculated and experimental binding free energies has been improved from 0.36 (for single-structure MM/GBSA rescoring) to 0.69 (for ensemble- average one). Decomposition of the calculated binding free energy reveals the electrostatic interactions in both solute and solvent play an important role in determining the binding free energy. The increasing negative charge of the compounds provides a more favorable electrostatic energy change but creates a higher penalty for the solvation free energy. Such a penalty is compensated by the electrostatic energy change, which results in a better binding affinity. The best binder has the highest ligand efficiency.