Skip to content
Open access

Benchmark Study for Calculations of pK a Values of Metal Ligands in Proteins

Jul 2026 · Journal of Chemical Theory and Computation · Vol 22, pp. 7420 - 7435 · 0 citations · 156 references
Medicine

Abstract

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.

Read PDF

Similar papers

Open access Jul 2026

Determination of the Formation Constants of Copper(II) Complexes Using Quantum Chemical Calculations

Gibbs free energy values obtained from quantum-chemical calculations were used to determine the formation constants of homo- and heteroligand copper(II) complexes with various ligands (amino acids, diimines, and phosphorylated dithiocarbamates) in an aqueous medium. A computationally robust, yet moderately expensive, level of theory—B3LYP/def2-TZVPPD—was employed. Solvent effects were accounted for using two models: C-PCM and SMD. A key prerequisite for obtaining reliable results is the use of a reference complex with a known formation constant that is structurally and solvation-wise similar to the compound under study. In this context, “similarity” implies identical stoichiometry, the same coordination number, the same number of water molecules in the inner coordination sphere; matching charges, however, is considerably less critical. The importance of considering conformers and isomers to obtain the most accurate values is demonstrated. The more rigid the structure and the greater the degree of similarity between the studied and reference compounds, the better the agreement between calculated and experimental data; the discrepancy can be as low as 0.1–0.2 logarithmic units. When an appropriate reference is selected, the average deviation of the calculated stability constants is less than 1 logarithmic unit.

N. S. Aksenin, Mikhail S. Bukharov, V. Shtyrlin et al. · 0 citations
Open access Jul 2026

Benchmarking semiempirical quantum chemical methods on liquid water.

Stimulated by the renewed interest and recent developments in semiempirical quantum chemical (SQC) methods for noncovalent interactions, we examine the properties of liquid water under ambient conditions by means of molecular dynamics (MD) simulations, both with the conventional neglect of diatomic differential overlap-type methods, e.g., AM1 and PM6, and with DFTB-type (density-functional tight-binding) methods, e.g., DFTB2 and GFN-xTB (Geometry-Frequency-Noncovalent eXtended Tight-Binding). Besides the original parameter sets, some specifically reparameterized SQC methods (denoted as AM1-W, PM6-fm, and DFTB2-iBi) targeting various smaller water systems ranging from molecular clusters to bulk are considered as well. The quality of these different SQC methods for describing liquid water properties under ambient conditions is assessed by comparison with well-established experimental data and also with BLYP-D3 density functional theory-based ab initio MD simulations. Our analyses reveal that static and dynamic properties of bulk water are poorly described by all considered SQC methods with the original parameters, regardless of the underlying theoretical models, with most of the methods suffering from too weak hydrogen bonds and hence predicting a far too fluid water with highly distorted hydrogen bond kinetics. Meanwhile, the reparameterized force-matched PM6-fm method is shown to be able to quantitatively reproduce the static and dynamic features of liquid water and thus can be used as a computationally efficient alternative to electronic structure-based MD simulations for liquid water that requires extended length and time scales. DFTB2-iBi predicts a slightly overstructured water with reduced fluidity, whereas AM1-W gives an amorphous ice-like structure for water under ambient conditions.

Xin Wu, Hossam Elgabarty, Vahideh Alizadeh et al. · 0 citations
Jul 2026

Efficient Force Evaluation via Fragmentation: Toward Geometry Optimization of Protein-Ligand Systems with Quantum Mechanical Accuracy.

This work presents a hybrid fragmentation-based approach, FragQMMM, that enables efficient and accurate force evaluation for large molecular systems, and successfully reproduces the hydrogen-bond network and preserving the underlying structure-activity relationship.

R. Han, Zonghua Bo, Jiawei Yan · 0 citations
Open access Jul 2026

Stability of Some Ternary 13-Atom Icosahedral Clusters Assessed with Geometric, Electronic, and Thermodynamic Criteria

Superatoms are clusters with physical or chemical properties that bear certain similarities to atoms. We calculated equilibrium geometries and properties using spin-polarized density functional theory (DFT) within the generalized gradient approximation (PBE). We also performed Ab Initio Molecular Dynamics (AIMD) simulations in the NVT ensemble with a Nosé-Hoover thermostat. In addition, Monte Carlo simulations with a potential fitted to DFT were carried out to assess a number of clusters as possible superatoms. The clusters we studied have a chemical composition AB2C10, where A, B, and C are metals, and a structure derived from that of a 13-atom icosahedron. We adopt an operational definition of superatom to assess clusters. According to this definition, superatoms should have (i) a quasispherical equilibrium geometry and (ii) a set of valence orbitals delocalized over the cluster; (iii) they should maintain these two attributes at elevated temperatures (mechanical and thermal stability); and (iv) they should be thermodynamically favored relative to similar clusters. By these criteria, Sn3Y10, CrSn2Zr10, and possibly CrSn2Ti10 can be classified as “superatoms”. We show that cluster stability arises from (i) the intrinsic stability of the C 10 fragment; (ii) its interaction energy with the three “dopant atoms”; and (iii) compatible atomic sizes of the three elements. Clusters such as CrSn2Ti10 and CrSn2Zr10 exhibit an optimal balance among these contributions. They have a high effective coordination, an absence of very low vibrational frequencies, a favorable ratio of their elements’ cohesive energy, and exceptional thermal resilience, with melting temperatures that are well above the weighted average of their elements’ bulk melting points. Electronic analysis shows that CrSn2Ti10 and CrSn2Zr10 display discrete density-of-states features and partial electron delocalization, consistent with superatomic behavior. These electronic properties vanish as the clusters approach the melting point and undergo structural distortion. We show the importance of looking at several, often interrelated, properties in assessing possible superatoms.

Anirudh Krishnadas, João Marcos Tomaz Palheta, Jonathan Bekele Mekonnen et al. · 0 citations
2026

Quantum Methods in Drug Discovery: An Introduction to the Theoceptor (Theoretical Receptor) Approach.

Theoceptor method can be applied to X-ray crystallography- or cryo-electron microscopy-derived structures, can be applied equally to covalent and noncovalent binders, and can provide excellent predictions for the changes in activity between matched molecular pairs.

Nikeel Cull, Yixin He, Sylvia C. Kaempf et al. · 0 citations
Open access Aug 2026

Evaluating molecular docking for binding affinity predictions: a systematic analysis of key parameters and the utility of AlphaFold2 structures for the Schrödinger dataset

Molecular docking is one of the most established methods in computational drug discovery, due to its balance of speed and accuracy. However, the accuracy of docking results depends on a number of different parameters, and systematic reference data for comparisons to more advanced methods for binding affinity prediction are still scarce. This study assesses the impact of key parameters on the accuracy of binding free energy estimates from docking, using nine benchmark systems with 278 high-affinity ligands. Using the Molecular Operating Environment (MOE), we evaluated combinations of three receptor structures (two crystal structures, one AlphaFold2 model), two force fields, two scoring functions, two receptor flexibility settings, and two statistical evaluation schemes. The performance of the docking approaches is measured based on the squared Pearson’s correlation coefficient (R²), the root mean square error (RMSE) with respect to the experimental binding affinities, as well as the mean signed error (MSE) and Kendall’s tau for individual targets and the full dataset. The results show that the scoring function and the protein structure are the most important factors for binding affinity accuracy in rigid docking with the MOE software. Amber10:EHT and MMFF94x force fields had the same average Rmean2 value, but Amber10:EHT had a lower average RMSEmean. AlphaFold2 protein models yielded lower binding affinity accuracy and higher errors compared to experimental crystal structures, although induced fit docking improved results. Using the original benchmark, we also compared several docking programs. DOCK6 and MOE performed best, with mean R² values of about 0.49 and 0.40, respectively. The remaining docking programs did not outperform a molecular weight regression baseline. For a subset of four targets (CDK2, JNK1, P38, TYK2) evaluated in previous work, the performance of the optimized DOCK6 and MOE protocols produced correlation coefficients similar to those reported for certain MM/PBSA, FMO, and Boltz2 implementations evaluated on the same target subset. This raises questions about potential dataset biases, the structural preparation, or the implementation of those methods. Docking therefore should be considered as an important and computationally inexpensive reference baseline for binding affinity prediction.

Konstantinos Tornesakis, J. Essex, Paul A. Cox et al. · 0 citations