The calculation of protein-ligand binding free energy (ΔG) is of great importance for virtual screening and drug design. Molecular dynamics (MD) simulation has been an attractive tool to investigate this scientific problem. However, the reliability of such approach is affected by many factors including electrostatic interaction calculation. Here, we present a practical protocol using quantum mechanics/molecular mechanics (QM/MM) calculations to generate polarizable QM protein charge (QMPC). The calculated QMPC of some atoms in binding pockets was obviously different from that calculated by AMBER ff03, which might significantly affect the calculated ΔG. To evaluate the effect, the MD simulations and MM/GBSA calculation with QMPC for 10 protein-ligand complexes, and the simulation results were then compared to those with the AMBER ff03 force field and experimental results. The correlation coefficient between the calculated ΔΔG using MM/GBSA under QMPC and the experimental data is .92, while that with AMBER ff03 force field is .47 for the complexes formed by streptavidin or its mutants and biotin. Moreover, the calculated ΔΔG with QMPC for the complexes formed by ERβ and five ligands is positively related to experimental result with correlation coefficient of .61, while that with AMBER ff03 charge is negatively related to experimental data with correlation coefficient of .42. The detailed analysis shows that the electrostatic polarization introduced by QMPC affects the electrostatic contribution to the binding affinity and thus, leads to better correlation with experimental data. Therefore, this approach should be useful to virtual screening and drug design.