Ab Initio Quantum Mechanical Modeling Research Articles

Ab initio generalized Langevin equation

We introduce a machine learning-based approach called ab initio generalized Langevin equation (AIGLE) to model the dynamics of slow collective variables (CVs) in materials and molecules. In this scheme, the parameters are learned from atomistic simulations based on ab initio quantum mechanical models. Force field, memory kernel, and noise generator are constructed in the context of the Mori-Zwanzig formalism, under the constraint of the fluctuation-dissipation theorem. Combined with deep potential molecular dynamics and electronic density functional theory, this approach opens the way to multiscale modeling in a variety of situations. Here, we demonstrate this capability with a study of two mesoscale processes in crystalline lead titanate, namely the field-driven dynamics of a planar ferroelectric domain wall, and the dynamics of an extensive lattice of coarse-grained electric dipoles. In the first case, AIGLE extends the reach of ab initio simulations to a regime of noise-driven motions not accessible to molecular dynamics. In the second case, AIGLE deals with an extensive set of CVs by adopting a local approximation for the memory kernel and retaining only short-range noise correlations. The scheme is computationally more efficient than molecular dynamics by several orders of magnitude and mimics the microscopic dynamics at low frequencies where it reproduces accurately the dominant far-infrared absorption frequency.

Predicting potential SARS-CoV-2 mutations of concern via full quantum mechanical modelling.

Ab initio quantum mechanical models can characterize and predict intermolecular binding, but only recently have models including more than a few hundred atoms gained traction. Here, we simulate the electronic structure for approximately 13 000 atoms to predict and characterize binding of SARS-CoV-2 spike variants to the human ACE2 (hACE2) receptor using the quantum mechanics complexity reduction (QM-CR) approach. We compare four spike variants in our analysis: Wuhan, Omicron, and two Omicron-based variants. To assess binding, we mechanistically characterize the energetic contribution of each amino acid involved, and predict the effect of select single amino acid mutations. We validate our computational predictions experimentally by comparing the efficacy of spike variants binding to cells expressing hACE2. At the time we performed our simulations (December 2021), the mutation A484K which our model predicted to be highly beneficial to ACE2 binding had not been identified in epidemiological surveys; only recently (August 2023) has it appeared in variant BA.2.86. We argue that our computational model, QM-CR, can identify mutations critical for intermolecular interactions and inform the engineering of high-specificity interactors.

Electrochemical/computational analysis and the derivation of Pourbaix diagrams of redox-active structural and adsorbed species in uranyl peroxide clusters (U60)

Electrochemical analysis of U60 nanoclusters (Li40K20[UO2(O2)(OH)]60(H2O)214) and their natural analog, the mineral studtite (UO2)O2(H2O)2·(H2O)2, was carried out using cyclic voltammetry with a powder microelectrode setup. Voltammetric analysis was supplemented by ab initio quantum mechanical modeling to understand the mechanism, thermodynamics, and structural changes during redox transition (with a semi-quantitative kinetics evaluation). This research aims to better understand the redox behavior and transport of dissolved and structural (incorporated) uranyl units in the environment. Further analysis is performed to determine if redox switching of actinyl ions in clusters can be performed while leaving the cluster intact or at what redox transition it disintegrates. Quantum-mechanical calculations were applied to approximate electrochemical peak potential shifts and add them to standard reduction potentials to more accurately assign peaks to specific redox transitions between different species. These peak potentials are calculated from oxidation state-dependent binding energies between different uranyl species either within or adsorbed to U60. This theoretical approach incorporates extensive error cancellation and serves to predict electrochemical redox potentials and identify the reaction taking place. Voltammograms of U60 in electrolyte solutions exhibit kinetically-inhibited coupled redox peaks assigned to the U(VI)/U(V) transition of cluster structural uranyl units at −0.34 V (vs. standard hydrogen electrode), indicative of U60 reduction from UO22+ to UO2+. The U(V)/U(IV) transition at −0.71 V is assigned to the subsequent U60 reduction of UO2+ to UO20. A method of observing these transitions in situ using synchrotron x-ray absorption spectroscopy was explored. Voltammetry, modeling, and peak analysis indicate a kinetically-inhibited two-step, one-e- (per step) reduction of U60 from U(VI) to U(V) to U(IV). Voltammetry and x-ray absorption spectroscopy measurements performed before and after exhaustive cycling and the highly stable nature of U60 provide evidence that clusters typically do not break apart upon reduction. Voltammograms of U60 collected in varied uranyl solutions are indicative of UO22+ adsorption to the cluster surface and its subsequent reduction at +0.24 V (U(VI)/U(V)) and −0.05 V (U(V)/U(IV)). For uranyl solutions in electrolyte, the two-step, one-e- (per step) reduction of U60 is superseded by uranyl adsorption from solution to the electrode surface with subsequent reduction via U(V) disproportionation. This study represents the sole electrochemical investigation of U60 in the literature and only the second of any uranyl peroxide nanocluster. It is also one of the first to use quantum–mechanical approaches in combination with existing solution-based standard reduction potentials to calculate Pourbaix (EH-pH) diagrams for adsorbed and structural species.

Absolute IR vibrational band intensities of hexafluoroacetone: Comparison of experiment and anharmonic ab initio calculation using the second-order operator canonical Van Vleck perturbation theory

Hexafluoroacetone (HFA, OC(CF3)2) is widely utilized in organic synthesis, chemistry of heat resistant polymers and biochemistry. In this work, infrared (IR) spectra of gaseous HFA were recorded in the range 2000–340 cm−1 with a resolution of 1 cm−1. For the first time, experimental absolute IR intensities of fourteen fundamental bands of HFA were measured and compared with intensities predicted by anharmonic ab initio calculations. Theoretical molecular structures, semi-diagonal quartic potential energy surfaces (PES), and cubic surfaces of dipole moment (DM) components of HFA were calculated using MP2/cc-pVTZ ab initio quantum-mechanical model. In addition, harmonic frequencies were also calculated at MP2/cc-pVQZ level of theory and used for refining the PES. The numerical-analytic implementation of the second-order operator canonical Van Vleck perturbation theory (CVPT2) was employed for predicting anharmonic IR spectra of HFA. The deviation from observations of predicted IR intensities for ten fundamentals (ν1, ν2, ν8, ν13–ν21) was less than 25% on average. The significant difference between ab initio harmonic intensities calculated with MP2/cc-pVTZ and MP2/cc-pVQZ models for some normal modes indicates that an inaccuracy of the quantum-mechanical model is probably the major source of errors in theoretical intensities. An interpretation of the IR experimental spectra and associated literature data for HFA is accomplished; the ν4, ν7, ν10 and ν15 fundamentals were reassigned. Three predicted Fermi resonances were identified in experimental spectra, ν3/ν5 + ν8, ν17/ν10 + ν19 and ν20/ν9 + ν23. An empirical fit of harmonic frequencies within anharmonic model assisted in a reliable identification of a number of overtones and binary combination tones. The fitted harmonic frequencies were compared with MP2/cc-pVQZ and a simplified complete basis set (CBS) extrapolation. It was found that MP2/cc-pVQZ frequencies were more reliable. The efficiency and convenience of the theoretical model employed and associated techniques for interpretation of complex spectra of HFA are demonstrated.

Ab initio modeling of Raman and infrared spectra of calcite

Calcite is the most thermodynamically stable polymorphic phase of the CaCO3 crystal. It is widely used in modern optical instruments operating in the infrared and visible wavelengths of electromagnetic radiation. In particular, due to its anisotropic properties, calcite is used in polarization optics devices. Ab initio quantum mechanical modeling of the Raman and infrared spectra of calcite makes it possible to better understand the structure and nature of the chemical bonds of the compound, and find the optimal conditions for the effective use of unique properties of calcite in photonics tools. The calculations are performed within the framework of the density functional theory in the CRYSTAL program, using the “hybrid” B3LYP functional and the all-electronic bases of the STO-6G, PO-TZVP and BSD atomic orbitals of the Gaussian type. The obtained results for the elastic constants, Raman and infrared spectra of calcite agree satisfactorily with the available experimental data for the basic sets of POB-TZVP and BSD.

Ab-initio quantum mechanical model for spin-strain coupling in ferroics

Ab-initio quantum mechanical model for spin-strain coupling in ferroics

Second Harmonic Generation Mediated by Aligned Water in Starch Granules.

The origin of second harmonic generation (SHG) in starch granules was investigated using ab initio quantum mechanical modeling and experimentally examined using polarization-in, polarization-out (PIPO) second harmonic generation microscopy. Ab initio calculations revealed that the largest contribution to the SHG signal from A- and B-type allomorphs of starch originates from the anisotropic organization of hydroxide and hydrogen bonds mediated by aligned water found in the polymers. The hypothesis was experimentally tested by imaging maize starch granules under various hydration and heat treatment conditions that alter the hydrogen bond network. The highest SHG intensity was found in fully hydrated starch granules, and heat treatment diminished the SHG intensity. The PIPO SHG imaging showed that dried starch granules have a much higher nonlinear optical susceptibility component ratio than fully hydrated granules. In contrast, deuterated starch granules showed a smaller susceptibility component ratio demonstrating that SHG is highly sensitive to the organization of the hydroxyl and hydrogen bond network. The polarization SHG imaging results of potato starch granules, representing starch allomorph B, were compared to those of maize starch granules representing allomorph A. The results showed that the amount of aligned water was higher in the maize granules. Nonlinear microscopy of starch granules provides evidence that varying hydration conditions leads to significant changes in the nonlinear susceptibility ratio as well as the SHG intensity, supporting the hypothesis from ab initio calculations that the dominant contribution to SHG is due to the ordered hydroxide and hydrogen bond network.

27Al NMR spectroscopy at multiple magnetic fields and ab initio quantum modeling for kaolinite

Diffraction methods indicate clearly that there are two crystallographic Al sites in kaolinite with different site symmetry. Many unsuccessful attempts have been made to resolve the two different Al sites in kaolinite by magic angle spinning nuclear magnetic resonance (MAS NMR) spectroscopy since the early 1980s. This study attempts to resolve these sites, derive accurate quadrupolar coupling parameters at the two sites, and explain the reasons for the lack of separation by combining NMR spectroscopy at various fields with ab initio quantum modeling. At a high field such as 21.06 T, large chemical shift anisotropy and dipolar effects can overwhelm the quadrupolar interaction effects for sites with relatively small C q values such as the two Al sites in kaolinite. The spectra will then display featureless lineshapes lacking the typical quadrupolar interaction characteristics. At lower fields, such as 4.68 T, the quadrupolar interaction effects become dominant, but due to the similarity of the two Al environments in kaolinite, as well as broadening by residual dipole-dipole and chemical shift anisotropy, the 27 Al MAS NMR peaks still could not be separated. To resolve the two different Al sites spectroscopically, 27 Al MAS NMR experiments at multiple magnetic fields, satellite transition NMR, and multiple-quantum MAS techniques have been utilized. Although these NMR experiments failed to resolve the two Al sites, the spectra have been interpreted with the aid of theoretical ab initio quantum mechanical modeling using full potential linearized augmented plane wave (FP-LAPW) model and the CASTEP software program. Average values of C q , η, and δ iso for the two Al sites were constrained to be 2.6 MHz, 0.75 and 6.25 ppm with high confidence.

Origins of the Different Metal Preferences of Escherichia coli Peptide Deformylase and Bacillus thermoproteolyticus Thermolysin: A Comparative Quantum Mechanical/Molecular Mechanical Study

The Escherichia coli peptide deformylase (PDF) and Bacillus thermoproteolyticus thermolysin (TLN) are two representative metal-requiring peptidases having remarkably similar active centers but distinctively different metal preferences. Zinc is a competent catalytic cofactor for TLN but not for PDF. Reaction pathways and the associated energetics for both enzymes were determined using combined semiempirical and ab initio quantum mechanical/molecular mechanical modeling, without presuming reaction coordinates. The results confirmed that both enzymes catalyze via the same chemical steps, and reproduced their different preferences for zinc or iron as competent cofactors. Further analyses indicated that different feasibility of the nucleophilic attack step leads to different metal preferences of the two enzymes. In TLN, the substrate is strongly activated and can serve as the fifth coordination ligand of zinc prior to the chemical steps. In PDF, the substrate carbonyl is activated by the chemical step itself, and becomes the fifth coordination partner of zinc only in a later stage of the nucleophilic attack. These leads to a much more difficult nucleophilic attack in PDF than in TLN. Different from some earlier suggestions, zinc has no difficulty in accepting an activated substrate as the fifth ligand to switch from tetra- to penta-coordination in either PDF or TLN. When iron replaces zinc, its stronger interaction with the hydroxide ligand may lead to higher activation barrier in TLN. In PDF, the stronger interactions of iron with ligands allow iron-substrate coordination to take place either before or at a very early stage of the chemical step, leading to effective catalysis. Our calculations also show combined semiempirical and ab initio quantum mechanical modeling can be efficient approaches to explore complicated reaction pathways in enzyme systems.

First-Principles Prediction of Equilibrium Potentials for Water Activation by a Series of Metals

In this paper, we describe and apply an ab initio quantum-mechanical model to examine the electrochemical dissociation of water over a number of close-packed and near-close-packed metal surfaces. We have examined a series of metals important for electrochemical energy generation and corrosion-resistant alloy development including Ni, Cu, Ru, Au, Mo, Pd, and Pt. Advances in electrochemical theory depend on a detailed understanding of the interplay between surface chemistry and electrochemical phenomena. Herein, we have used ab initio quantum-mechanical methods to examine the double-layer regions for over these metals and compare the equilibrium potentials for the initial steps of water reduction and oxidation at the surface with known experimental quantities. We find that in its current form, the model developed herein is semiquantitative: it allows for the correct prediction of trends and the size of the double-layer regions, but in some cases it results in significant deviation with known absolute equilibrium potentials. Additional discussion is provided that outlines steps that may be taken to improve the quantitative accuracy of this model.

Mechanism for the dissolution of olivine series minerals in acidic solutions

The complexity of silicate minerals makes prediction of their dissolution rates a challenging problem. A combination of large cluster ab-initio quantum mechanical models and chemical probe dissolution experiments are used to understand the dissolution process for olivine-group minerals. Rapid release of M2+ cations by precursor reactions involving H+ attack at µ3-O surface sites produces a silica-enriched surface. Slower rates of silica release via a ligand exchange reaction involving a proton in the activated complex controls the overall rate of olivine dissolution. Our results provide a physical explanation for the correlation among olivine dissolution rates and water exchange rates for the corresponding aqueous cation.

Structural Rationale for Low-Nanomolar Binding of Transition State Mimics to a Family GH3 β-d-Glucan Glucohydrolase from Barley,

The interactions of a transition state mimic anilinomethyl glucoimidazole (AmGlcIm), with a K(i) constant of 0.6 x 10(-)(9) M and a Gibbs free energy value of -53.5 kJ/mol, with a family GH3 beta-d-glucan glucohydrolase from barley have been analyzed crystallographically and by ab initio quantum mechanical modeling. AmGlcIm binds 3 times more tightly to the beta-d-glucan glucohydrolase than a previously investigated phenyl glucoimidazole. In the enzyme-AmGlcIm complex, an additional residue, Tyr253, and a water molecule positioned between subsites -1 and +1 are recruited for binding. Analyses of the two binary complexes reveal the following. (i) An intricate network exists in which hydrogen bonds between the enzyme's catalytic pocket residues Lys206, His207, Tyr253, Asp285, and Glu491 and the glucoimidazoles are shorter by 0.15-0.53 A, compared with distances of hydrogen bonds in the Michaelis complex. (ii) The "glucose" moiety of the glucoimidazoles adopts a (4)E conformation that is vital for the low-nanomolar binding. (iii) The N1 atoms of the glucoimidazoles are positioned nearly optimally for in-line protonation by the Oepsilon1 atom of the catalytic acid/base Glu491. (iv) The enzyme derives binding energies from both glycone and aglycone components of the glucoimidazoles. (iv) The prevalent libration motion of the two domains of the enzyme could play a significant role during induced fit closure in the active site. (v) Modeling based on the structural data predicts that protons could be positioned on the N1 atoms of the glucoimidazoles, and the catalytic acid/base Glu491 could carry an overall negative charge. (vi) The enzyme-AmGlcIm complex reveals the likely structure of an early transition state during hydrolysis. Finally, the high-resolution structures enabled us to define minimal structures of oligosaccharides attached to Asn221, Asn498, and Asn600 N-glycosylation sites.

Ab initio quantum mechanical modeling of infrared vibrational frequencies of the OH group in dioctahedral phyllosilicates. Part II: Main physical factors governing the OH vibrations

The physical factors responsible for the variability observed in OH infrared (IR) fundamentals in dioctahedral phyllosilicates, due to octahedral substitution of Al3+ by Mg2+, Fe2+, and Fe3+, are discussed here. The data analyzed consist of experimental frequencies as well as frequencies modeled using Density Functional Theory (DFT) calculations. The charge of the octahedral cations surrounding the OH is one of the main factors affecting both the OH stretch and the in-plane bend; cationic electronegativity and ionic radius play important roles in the stretch and bend modes, respectively. The mass of the octahedral cations does not affect the OH fundamental vibrations. The nature of the octahedral cations alone can explain most of the variability observed in the OH in-plane bend, making this fundamental vibration the most suitable for assessing octahedral composition. Discrepancies between modeled and experimental OH stretch frequencies indicate the existence of other factors governing this fundamental vibration. Further DFT calculations indicate that apical O atoms of the tetrahedral sheet with unsatisfied charges due to octahedral and/or tetrahedral substitutions can explain these discrepancies. The modeling results are utilized to predict the frequency of the OH stretch and in-plane-bend combination band that occurs near 4545 cm−1 (2.2 μm) in phyllosilicates. This band can be observed in imaging spectrometer data, allowing for the detection and analysis of phyllosilicates and other minerals in large natural systems. The modeling results confirm that the variability observed in the combination band of dioctahedral phyllosilicates reflects octahedral and, to a certain degree, tetrahedral composition, but not interlayer composition.

Ab initio quantum mechanical models of peptide helices and their vibrational spectra.

Structural parameters for standard peptide helices (alpha, 3(10), 3(1) left-handed) were fully ab initio optimized for Ac-(L-Ala)(9)-NHMe and for Ac-(L-Pro)(9)-NHMe (poly-L-proline-PLP I and PLP II-forms), in order to better understand the relative stability and minimum energy geometries of these conformers and the dependence of the ir absorption and vibrational CD (VCD) spectra on detailed variation in these conformations. Only the 3(10)-helical Ala-based conformation was stable in vacuum for this decaamide structure, but both Pro-based conformers minimized successfully. Inclusion of solvent effects, by use of the conductor-like screening solvent model (COSMO), enabled ab initio optimizations [at the DFT/B3LYP/SV(P) level] without any constraints for the alpha- and 3(10)-helical Ala-based peptides as well as the two Pro-based peptides. The geometries obtained compare well with peptide chain torsion angles and hydrogen-bond distances found for these secondary structure types in x-ray structures of peptides and proteins. For the simulation of VCD spectra, force field and intensity response tensors were obtained ab initio for the complete Ala-based peptides in vacuum, but constrained to the COSMO optimized torsional angles, due to limitations of the solvent model. Resultant spectral patterns reproduce well many aspects of the experimental spectra and capture the differences observed for these various helical types.

Ab initio quantum mechanical modeling of infrared vibrational frequencies of the OH group in dioctahedral phyllosilicates. Part I: Methods, results and comparison to experimental data

The infrared (IR) spectra of small clusters of atoms ([MM′(OH)2] and [MM′ (OH)2 (H2O)6], where M, M′ = Al3+, Mg2+, Fe2+, Fe3+) mimicking the environment of the OH group in dioctahedral phyllosilicates have been modeled using ab initio quantum mechanical calculations. These modeling results are relevant to establishing the connections between IR spectra of phyllosilicates and their composition, and to investigate the utility of quantum mechanical models for calculating IR frequencies of minerals. This study focused on the OH stretch and in-plane bend fundamentals, because they give rise to a combination band near 4545 cm−1 (2.2 μm) that can be observed in imaging spectrometer or hyperspectral remote sensing data. A comparison among results obtained using both ab initio [Hartree-Fock (HF) and Density Functional Theory (DFT)], and semi-empirical [PM3(tm)] methods, showed that the DFT model approaches IR frequency experimental values most closely. IR spectra of phyllosilicates were modeled using the DFT method. The modeled frequencies were scaled using a mode-dependent linear transformation, and experimental frequencies were reproduced satisfactorily. The modeling results show that most of the variability observed in the OH in-plane bend fundamental of dioctahedral phyllosilicates can be explained by the effects of neighboring octahedral cations alone. Discrepancies between modeling and experimental results in the case of the OH stretch point to the existence of factors other than the nature of the neighboring octahedral cations, such as tetrahedral substitution, affecting this fundamental mode.

Determination of effective pair potentials from ab initio simulations: application to liquid water

We present a novel approach to construct pairwise site–site intermolecular interaction potentials from the results of ab initio quantum-mechanical modeling to be used in classical molecular simulations. In the present application, as a first step, we have carried out ab initio Car–Parrinello molecular dynamics (MD) simulations of liquid water to calculate the radial distribution functions (RDF) between all the pairs of nuclei. A second step involves use of the inverse Monte Carlo method, recently developed by us [Phys. Rev. E 52 (1995) 3730], to construct a full set of site–site interaction potentials from the ab initio RDFs. The resulting three-site SPC-like water model, which reproduces the liquid structure obtained in the Car–Parrinello MD simulations, can be used as a potential function in classical molecular simulations. We have also found that there exist no pair potential for a rigid three-site water model which exactly satisfy the experimental water RDFs.

An ab initio study of the comparative stabilities of the simple phosphine oxide tautomers

High level ab initio quantum mechanical model chemistries have been used to study the relative stabilities of the simple phosphine oxides tautomers OPH3/PH2OH and OPH(OH)2/P(OH)3 and related molecules. It is found that the two PH2OH isomers are virtually equi-energetic and only slightly more stable than OPH3 itself, and the oxides with one or more hydroxyl groups are significantly more stable than their hydroxide tautomers. The difference in stability arises mainly from the variation of the P=O bond energy among the various species, the P=O bond energy increasing with the number of hydroxyl groups present. The effects of solvation in aqueous solution are estimated by reaction field polarized continuum models and by explicit calculations of OPH3 and PH2OH hydrogen bound to a single water molecule. Both approaches indicate that the oxide form is significantly stabilized by its aqueous environment. © 2000 John Wiley & Sons, Inc. Heteroatom Chem 11:73–80, 2000

Torsional Control of Epoxidation Stereoselectivity in 1,2-Dihydronaphthalenes. Transition State Modeling with Semiempirical Quantum Mechanics

Recently, highly diastereoselective epoxidations of chiral 1,2-dihydronaphthalenes were reported. Crystallographic characterization of the major product indicated that the selectivity cannot be attributed to long-range steric effects. We examined this system through ab initio and semiempirical quantum mechanical modeling. Corroborating experimental data, our investigations explain the apparent irrelevance of direct interactions with substituents and address the insensitivity of the stereoselectivity to oxidant structure. Furthermore, extension of this model to other related epoxidation reactions indicates that pi-facial discrimination may be interpreted as a result of a preference for a staggered asynchronous transition state, obtainable only through axial attack.

Ab Initio Quantum Mechanical Model Calculations on the Catalytic Mechanism of Aspartylglucosaminidase (AGA): A Serine Protease‐Like Mechanism with an N‐terminal Threonine and Substrate‐Assisted Catalysis

AbstractAb initio quantum mechanical model calculations were used in studying the acylation and deacylation steps of the catalytic mechanism of aspartylglucosaminidase (AGA). AGA catalyses the hydrolysis of an amide linkage between oligosaccharide and asparagine by utilising an N‐terminal threonine as a catalytic amino acid. Results are reported for the model enzyme reaction at the MP2/6–31+G*//HF/6–31 + G* + ΔZPE level. Contribution of aqueous solvation to the relative energies was estimated by using the continuum solvation model of Tomasi. The serine protease‐like catalytic mechanism was found to be feasible for AGA. The protonated α‐amino group of the substrate of AGA was suggested to enhance the catalysis by stabilising the anionic oxygen of the substrate, which is formed in the reaction, and by lowering the pKa of the nucleophilic oxygen of the N‐terminal threonine. Finally, the similarities in the catalytic mechanisms of AGA and other amidohydrolases were discussed.

Model assembly study of the ligand binding by p-hydroxybenzoate hydroxylase: correlation between the calculated binding energies and the experimental dissociation constants.

The energies of binding of seven ligands by p-hydroxybenzoate hydroxylase (PHBH) were calculated theoretically. Direct enzyme-ligand interaction energies were calculated using the ab initio quantum mechanical model assembly of the active site at the 3-21G level. Solvation energies of the ligands needed in the evaluation of the binding energies were calculated with the semiempirical AM1-SM2 method and the long-range electrostatic interaction energies between the ligands and the protein matrix classically using the static charge distributions of the ligands and the protein. Energies for proton-transfer between the ligands' OH or SH substituent at position 4 and the active-site tyrosine within the ab initio model assemblies were calculated and compared to the corresponding pKas in aqueous solution. Excluding 3,4-dihydroxybenzoate, the natural product of PHBH, a linear relationship between the calculated binding energies and the experimental binding free energies was found with a correlation coefficient of 0.90. Contributions of the direct enzyme-ligand interaction energies, solvation energies and the long-range electrostatic interaction energies to the calculated binding energies were analyzed. The proton-transfer energies of the ligands with substituents ortho to the ionized OH were found to be perturbed less in the model calculations than the energies of their meta isomers as deduced from the corresponding pKas.