Benchmark Study for Calculations of pK a Values of Metal Ligands in Proteins
We have compared the performance of 64 different computational methods, based on combined quantum mechanical (QM) and molecular mechanical (QM/MM) or QM-cluster calculations in a continuum solvent, to estimate the acid constant (pK a) of metal-bound ligands in proteins. As a calibration set, we use 12 experimental pK a values from six different proteins that involve Zn2+, Fe3+, or Fe4+. We employ two different density functional theory (DFT) methods (TPSS and B3LYP), two basis sets (def2-SV(P) and def2-TZVPD), QM regions of three different sizes (∼40, ∼100, and ∼350 atoms), relaxed or fixed surroundings, and three different values of the dielectric constant of the continuum-solvation model (ε = 4, 20, or 80). The results clearly show that QM-cluster+continuum-solvation is much better than QM/MM. In general, the most accurate results are obtained with ε = 80 and the minimal QM region. The two DFT methods, the two basis sets, and relaxing or fixing the surroundings give similar results. The best-performing method is TPSS with the minimal QM region, def2-TZVPD, relaxed surroundings, and ε = 80, yielding a mean absolute deviation (after removal of a systematic error of 11.6 pK a units, pu) of 2.0 pu and a maximum deviation of 5.0 pu. The coefficient of determination (R 2) and Kendall’s τ with respect to the experimental pK a values are both 0.64, while Spearman’s rank correlation coefficient is 0.78. This level of accuracy should be sufficient to reliably determine the protonation states of metal-bound ligands in QM-based studies of enzymatic reaction mechanisms.