We present an adaptation of the MM-PB/SA method for the estimation of binding free energies in protein-ligand complexes simulated with QM/MM molecular dynamics. The method is applied to understand the binding of a set of tetrabromobenzimidazole inhibitors to the CK2 protein. We find that the QM/MM interaction energy alone cannot always be used as a predictor of the binding affinity, and the inclusion of solvation effects via the PB/SA method is essential in getting reliable results. In agreement with experimental observations, we show that the van der Waals interactions are the driving force for the binding, while the electrostatic interactions orient these inhibitors in the CK2 active site. Additionally a per-residue energy decomposition analysis was applied to determine the individual contributions to the protein-inhibitor interaction. Based on these results, we hypothesize that the inclusion of a sufficiently large polar group on the tetrabromobenzimidazole skeleton could increase the binding affinity. The results show that the QM/MM-PB/SA method can be successfully employed to understand complicated structure-activity relationships and to design new inhibitors.