Protein–ligand free energies of binding from full-protein DFT calculations: convergence and choice of exchange–correlation functional

The accurate prediction of protein–ligand binding free energies with tractable computational methods has the potential to revolutionize drug discovery. Modeling the protein–ligand interaction at a quantum mechanical level, instead of relying on empirical classical-mechanics methods, is an important step toward this goal. In this study, we explore the QM-PBSA method to calculate the free energies of binding of seven ligands to the T4-lysozyme L99A/M102Q mutant using linear-scaling density functional theory on the whole protein–ligand complex. By leveraging modern high-performance computing we perform over 2900 full-protein (2600 atoms) DFT calculations providing new insights into the convergence, precision and reproducibility of the QM-PBSA method. We find that even at moderate sampling over 50 snapshots, the convergence of QM-PBSA is similar to traditional MM-PBSA and that the DFT-based energy evaluations are very reproducible. We show that in the QM-PBSA framework, the physically-motivated GGA exchange–correlation functional PBE outperforms the more modern, dispersion-including non-local and meta-GGA-nonlocal functionals VV10 and B97M-rV. Different empirical dispersion corrections perform similarly well but the three-body dispersion term, as included in Grimme's D3 dispersion, is significant and improves results slightly. Inclusion of an entropy correction term sampled over less than 25 snapshots is detrimental while an entropy correction sampled over the same 50 or 100 snapshots as the enthalpies improves the accuracy of the QM-PBSA method. As full-protein DFT calculations can now be performed on modest computational resources our study demonstrates that they can be a useful addition to the toolbox of free energy calculations.