Theoretical Calculations of Acid Dissociation Constants: A Review Article

Impact of Intracellular Proteins on μ-Opioid Receptor Structure and Ligand Binding

Chronic pain is a prevalent problem affecting approximately one out of every five adults in the U.S. The most effective way to treat chronic pain is with opioids, but they cause dangerous side effects such as tolerance, addiction, and respiratory depression, which makes them quite deadly. Opioids, such as fentanyl, target the μ-opioid receptor (MOR), which can then bind to the intracellular Gi protein or the β-arrestin protein. The Gi pathway is primarily responsible for pain relief and potential side effects, but the β-arrestin pathway is chiefly responsible for the unwanted side effects. Ideally, an effective pain medication without side effects would bind to MOR, which would bias signaling solely through the Gi pathway. We used the Bio3D library to conduct principal component analysis to compare the cryo-electron microscopy MOR structure in complex with the Gi versus an X-ray crystallography MOR structure with a nanobody acting as a Gi mimic. Our results agree with a previous study by Munro, which concluded that nanobody-bound MOR is structurally different than Gi-bound MOR. Furthermore, we investigated the structural diversity of opioids that can bind to MOR. Quantum mechanical calculations show that the low energy solution structures of fentanyl differ from the one bound to MOR in the experimental structure, and pKa calculations reveal that fentanyl is protonated in aqueous solution. Glide docking studies show that higher energy structures of fentanyl in solution form favorable docking complexes with MOR. Our calculations show the relative abundance of each fentanyl conformation in solution as well as the energetic barriers that need to be overcome to bind to MOR. Docking studies confirm that multiple fentanyl conformations can bind to the receptor. Perhaps a variety of conformations of fentanyl can stabilize multiple conformations of the MOR, which can explain why fentanyl can induce different intracellular signaling and multiple physiological effects.

A Comprehensive Computational Study on the Thermodynamics and Kinetics of Tetrahydrobiopterin Regeneration Process

One of the most crucial enzymatic cofactors in the human body is tetrahydrobiopterin, which is acquired through biological synthesis and self-regeneration. During this regenerative process, it undergoes oxidation, deprotonation, further oxidation, and subsequent deprotonation, resulting in the formation of quinonoid-dihydrobiopterin, which then undergoes tautomerization to yield dihydrobiopterin. This study presents the thermodynamic and kinetic properties associated with each stage of the regeneration process using theoretical calculations. The redox potentials for oxidation steps and the pKa values for deprotonation steps are determined employing the Born-Haber cycle and the direct change of free energy in implicit solvent models. The redox metabolites are characterized and confirmed from their calculated absorption spectra using the time-dependent density functional theory method. For the tautomerization steps, an IRC calculation is executed, and rate constants are computed using Eyring's Transition State Theory (TST). The tunnelling probability of the H atom during the tautomerization process is incorporated using Wigner's tunnelling correction in the calculation of the rate constant. Notably, we identify the N3 atom as the most probable deprotonation site for H3B+ and predict its geometry based on our calculations. Furthermore, we elucidate the spectral properties of intermediates involved in the regeneration process, highlighting key electronic transitions responsible for their excitations. Our results indicate that each step of tautomerization occurs along vibrational bending modes. We have observed that these tautomerization processes have high activation energies by optimising transition states. Additionally, considering tunnelling correction can significantly affect the reaction rates associated with these processes. These results provide a comprehensive understanding of the thermodynamics and kinetics of the regeneration process of tetrahydrobiopterin, which will help in the modulation of its biological activity.

Automated preparation of nanoscopic structures: Graph-based sequence analysis, mismatch detection, and pH-consistent protonation with uncertainty estimates

Structure and function in nanoscale atomistic assemblies are tightly coupled, and every atom with its specific position and even every electron will have a decisive effect on the electronic structure, and hence, on the molecular properties. Molecular simulations of nanoscopic atomistic structures therefore require accurately resolved three-dimensional input structures. If extracted from experiment, these structures often suffer from severe uncertainties, of which the lack of information on hydrogen atoms is a prominent example. Hence, experimental structures require careful review and curation, which is a time-consuming and error-prone process. Here, we present a fast and robust protocol for the automated structure analysis and pH-consistent protonation, in short, ASAP. For biomolecules as a target, the ASAP protocol integrates sequence analysis and error assessment of a given input structure. ASAP allows for pK a  prediction from reference data through Gaussian process regression including uncertainty estimation and connects to system-focused atomistic modeling described in Brunken and Reiher (J. Chem. Theory Comput. 16, 2020, 1646). Although focused on biomolecules, ASAP can be extended to other nanoscopic objects, because most of its design elements rely on a general graph-based foundation guaranteeing transferability. The modular character of the underlying pipeline supports different degrees of automation, which allows for (i) efficient feedback loops for human-machine interaction with a low entrance barrier and for (ii) integration into autonomous procedures such as automated force field parametrizations. This facilitates fast switching of the pH-state through on-the-fly system-focused reparametrization during a molecular simulation at virtually no extra computational cost.

Identification of Fouling Occurring during Coupled Electrodialysis and Bipolar Membrane Electrodialysis Treatment for Tofu Whey Protein Recovery

Tofu whey, a by-product of tofu production, is rich in nutrients such as proteins, minerals, fats, sugars and polyphenols. In a previous work, protein recovery from tofu whey was studied by using a coupled environmental process of ED + EDBM to valorize this by-product. This process allowed protein recovery by reducing the ionic strength of tofu whey during the ED process and acidifying the proteins to their isoelectric point during EDBM. However, membrane fouling was not investigated. The current study focuses on the fouling of membranes at each step of this ED and EDBM process. Despite a reduction in the membrane conductivities and some changes in the mineral composition of the membranes, no scaling was evident after three runs of the process with the same membranes. However, it appeared that the main fouling was due to the presence of isoflavones, the main polyphenols in tofu whey. Indeed, a higher concentration was observed on the AEMs, giving them a yellow coloration, while small amounts were found in the CEMs, and there were no traces on the BPMs. The glycosylated forms of isoflavones were present in higher concentrations than the aglycone forms, probably due to their high amounts of hydroxyl groups, which can interact with the membrane matrices. In addition, the higher concentration of isoflavones on the AEMs seems to be due to a combination of electrostatic interactions, hydrogen bonding, and π-π stacking, whereas only π-π stacking and hydrogen bonds were possible with the CEMs. To the best of our knowledge, this is the first study to investigate the potential fouling of BPMs by polyphenols, report the fouling of IEMs by isoflavones and propose potential interactions.

Dissociation constants of relevant secondary organic aerosol components in the atmosphere

The presented studies focus on measuring the determination of the acidity constant (pKa) of relevant secondary organic aerosol components. For our research, we selected important oxidation products (mainly carboxylic acids) of the most abundant terpene compounds, such as α-pinene, β-pinene, β-caryophyllene, and δ-3-carene. The research covered the synthesis and determination of the acidity constant of selected compounds. We used three methods to measure the acidity constant, i.e., 1H NMR titration, pH-metric titration, Bates-Schwarzenbach spectrophotometric method. Moreover, the pKa values were calculated with Marvin 21.17.0 software to compare the experimentally derived values with those calculated from the chemical structure. pKa values measured with 1H NMR titration ranged from 3.51 ± 0.01 for terebic acid to 5.18 ± 0.06 for β-norcaryophyllonic acid. Moreover, the data determined by the 1H NMR method revealed a good correlation with the data obtained with the commonly used potentiometric and UV-spectroscopic methods (R2 = 0.92). In contrast, the comparison with in silico results exhibits a relatively low correlation (R2Marvin = 0.66). We found that most of the values calculated with the Marvin Program are lower than experimental values obtained with pH-metric titration with an average difference of 0.44 pKa units. For di- and tricarboxylic acids, we obtained two and three pKa values, respectively. A good correlation with the literature values was observed, for example, Howell and Fisher (1958) used pH-metric titration and measured pKa1 and pKa2 to be 4.48 and 5.48, while our results are 4.24 ± 0.10 and 5.40 ± 0.02, respectively.

Determination of Amino Acids' pKa: Importance of Cavity Scaling within Implicit Solvation Models and Choice of DFT Functionals

Protonation states of molecules significantly influence the thermodynamics and kinetics of chemical reactions. This is especially important in biochemical processes, where appropriate protonation states of amino acids control the exo/endoergicity of practically all biochemical cycles. This paper is focused on appraisal of the impact of DFT functionals and PCM solvation models on the accuracy of pKa evaluations for all proteinogenic amino acids. Eight functionals (B3LYP, PBE0, revPBE0, M06-2X, M11, M11-L, TPSSh, and ωB97X-D) and four basis sets are considered, together with four kinds of implicit solvation models when additional attention is paid to a cavity construction. An influence of nonelectrostatic contributions and Wertz's corrections on Gibbs free energy is investigated together with accuracy of provided proton solvation energy. The best model is based on the M06-2X/6-311++G**/D-PCM/UAKS computational level. The fitting procedure is utilized to improve the accuracy of the evaluated models. All of these results are also compared with values obtained from the COSMOtherm program and CCSD(T) calculations. Results for cysteine and histidine are discussed individually, as they can be found in different protonation states at neutral pH.

Leveraging DFT and Molecular Fragmentation for Chemically Accurate pKa Prediction Using Machine Learning

We present a quantum mechanical/machine learning (ML) framework based on random forest to accurately predict the pKas of complex organic molecules using inexpensive density functional theory (DFT) calculations. By including physics-based features from low-level DFT calculations and structural features from our connectivity-based hierarchy (CBH) fragmentation protocol, we can correct the systematic error associated with DFT. The generalizability and performance of our model are evaluated on two benchmark sets (SAMPL6 and Novartis). We believe the carefully curated input of physics-based features lessens the model's data dependence and need for complex deep learning architectures, without compromising the accuracy of the test sets. As a point of novelty, our work extends the applicability of CBH, employing it for the generation of viable molecular descriptors for ML.

Accurate Computation of Aqueous pKas of Biologically Relevant Organic Acids: Overcoming the Challenges Posed by Multiple Conformers, Tautomeric Equilibria, and Disparate Functional Groups with the Fully Black-Box pK-Yay Method

The use of static electronic structure calculations to compute solution-phase pKas offers a great advantage in that a macroscopic bulk property could be computed via microscopic computations involving very few molecules. There are various sources of errors in the quantum chemical calculations though. Overcoming these errors to accurately compute pKas of a plethora of acids is an active area of research in physical chemistry pursued by both computational as well as experimental chemists. We recently developed the pK-Yay method in our attempt to accurately compute aqueous pKas of strong and weak acids. The method is fully black-box, computationally inexpensive, and is very easy for even a nonexpert to use. However, the method was thus far tested on very few molecules (only 16 in all). Herein, in order to assess the future applicability of pK-Yay, we study the effect of multiple conformers, the presence of tautomers under equilibrium, and the impact of a wide variety of functional groups (derivatives of acetic acid with substituents at various positions, dicarboxylic acids, aromatic carboxylic acids, amines and amides, phenols and thiols, and fluorine bearing organic acids). Starting with more than 1000 conformers and tautomers, this study establishes that overall errors of ∼ 1.0 pKa units are routinely obtained for a majority of the molecules. Larger errors are noted in cases where multiple charges, intramolecular hydrogen bonding, and several ionizable functional groups are simultaneously present. An important conclusion to emerge from this work is that, the computed pKas are insensitive (difference <0.5) to whether we consider multiple conformers/tautomers or only choose the most stable conformer/tautomer. Further, pK-Yay captures the stereoelectronic effects arising due to differing axial vs equatorial pattern, and is useful to predict the dominant acid-base equilibrium in a system featuring several equilibria. Overall, pK-Yay may be employed in several chemical applications featuring organic molecules and biomonomers.

Theoretical Study on the Alkylimino-Substituted Sulfonamides with Potential Biological Activity

Antibiotics play a key role in the fight against bacterial diseases. However, bacteria quickly learn how to minimize the effects of antibiotics and strengthen their resistance. Thus, the fight against them becomes more and more difficult and there is a constant search for new bactericidal compounds. It is important in this type of search to determine the basic properties of compounds such as pKa, hydrogen bond formation, or hydrophobicity. Here, we present the results of our in silico study of five sulfonamide derivatives differing in alkylamine substituent length. Based on our results, we propose a model of three possible pKa values for each of the studied compounds. Interestingly, the use of Muckerman's approach for pKa determination exhibits that theoretical and experimental results are in very good agreement. Intramolecular hydrogen bond formation affects pKa. The strength of the H-bond interaction increases from ethyl to butylamine and then decreases with the elongation of the alkylamine chain. The obtained partition coefficients (expressed here in the value of log P) increase with the number of carbon atoms in the alkylamine chain following Lipinski's rule of five. The presented results provide important structural, physicochemical, and thermodynamic information that allows for the understanding of the influence of some sulfonamides and their possible activity.

Approaching Coupled Cluster Accuracy with Density Functional Theory Using the Generalized Connectivity-Based Hierarchy

This Perspective reviews connectivity-based hierarchy (CBH), a systematic hierarchy of error-cancellation schemes developed in our group with the goal of achieving chemical accuracy using inexpensive computational techniques ("coupled cluster accuracy with DFT"). The hierarchy is a generalization of Pople's isodesmic bond separation scheme that is based only on the structure and connectivity and is applicable to any organic and biomolecule consisting of covalent bonds. It is formulated as a series of rungs involving increasing levels of error cancellation on progressively larger fragments of the parent molecule. The method and our implementation are discussed briefly. Examples are given for the applications of CBH involving (1) energies of complex organic rearrangement reactions, (2) bond energies of biofuel molecules, (3) redox potentials in solution, (4) pK_a predictions in the aqueous medium, and (5) theoretical thermochemistry combining CBH with machine learning. They clearly show that near-chemical accuracy (1-2 kcal/mol) is achieved for a variety of applications with DFT methods irrespective of the underlying density functional used. They demonstrate conclusively that seemingly disparate results, often seen with different density functionals in many chemical applications, are due to an accumulation of systematic errors in the smaller local molecular fragments that can be easily corrected with higher-level calculations on those small units. This enables the method to achieve the accuracy of the high level of theory (e.g., coupled cluster) while the cost remains that of DFT. The advantages and limitations of the method are discussed along with areas of ongoing developments.

Electronic Substitution Effect on the Ground and Excited State Properties of Indole Chromophore: A Computational Study

Indole, being the main chromophore of amino acid tryptophan and several other biologically relevant molecules like serotonin, melatonin, has prompted considerable theoretical and experimental interest. The current work focuses on the investigation of substitution effect on the ground and excited electronic states of indole using computational quantum chemistry. Having three close-lying excited electronic states, the vibronic coupling effect becomes extremely important yet challenging for the photophysics and photochemistry of indole. Here, we have evaluated the performance of time-dependent density functional theory against available experimental and ab initio results from the literature. The electronic effects on the excited states of indole and indole derivatives e. g. tryptophan, serotonin and melatonin are reported. A bathochromic shift has been observed in the absorption spectrum for the L_a state. The absorption wavelength increases in the order of indole<tryptophan <serotonin <melatonin. While the contribution of the in-plane small adjacent groups increases the electron density of the indole ring, the out-of-plane long substituent groups have minor effect. The absorption spectra calculated including the vibronic coupling are in good agreement with experiments. These results can be used to estimate the error in photophysical observables of indole derivatives calculated considering indole as a prototypical system.

Universal Trends between Acid Dissociation Constants in Protic and Aprotic Solvents

pKa values in non-aqueous solvents are of critical importance in many areas of chemistry. Our knowledge is, despite their relevance, still limited to the most fundamental properties and few pKa values in the most common solvents. Taking advantage of a recently introduced computationally efficient procedure we computed the pKa values of 182 compounds in 21 solvents. This data set is used to establish for the first time universal trends between all solvents. Our computations indicate, that the total charge of the molecule and the charge of the acidic group combined with the Kamlet-Taft solvatochromic parameters are sufficient to predict pKa values with at least semi- quantitative accuracy. We find, that neutral acids such as alcohols are strongly affected by the solvent properties. This is contrasted by cationic acids like ammonium ions whose pKa is often almost completely independent from the choice of solvent.

Predicting octanol/water partition coefficients and pKa for the SAMPL7 challenge using the SM12, SM8 and SMD solvation models

Blind predictions of octanol/water partition coefficients and pKa at 298.15 K for 22 drug-like compounds were made for the SAMPL7 challenge. Octanol/water partition coefficients were predicted from solvation free energies computed using electronic structure calculations with the SM12, SM8 and SMD solvation models. Within these calculations we compared the use of gas- and solution-phase optimized geometries of the solute. Based on these calculations we found that in general the use of solution phase-optimized geometries increases the affinity of the solutes for water as compared to octanol, with the use of gas-phase optimized geometries resulting in the better agreement with experiment. The pKa is computed using the direct approach, scaled solvent-accessible surface model, and the inclusion of an explicit water molecule, where the latter two methods have previously been shown to offer improved predictions as compared to the direct approach. We find that the use of an explicit water molecule provides superior predictions, and that the predicted macroscopic pKa is sensitive to the employed microstates.

Theoretical Study on the Photoacidity of Hydroxypyrene Derivatives in DMSO Using ADC(2) and CC2

This work applies the thermodynamic Förster cycle to theoretically investigate the pK_a^*, i.e., excited-state pK_a values of pyranine-derived superphotoacids developed by Jung and co-workers. The latter photoacids are strong enough to transfer a proton to the aprotic solvent dimethyl sulfoxide (DMSO). The Förster cycle provides access to pK_a^* via the ground-state pK_a and the electronic excitation energies. We use the conductor-like screening model for real solvents (COSMO-RS) to compute the ground-state pK_a and the correlated wavefunction-based methods ADC(2) and CC2 with the continuum solvation model COSMO to calculate the pK_a change upon excitation. A comparison of the calculated UV/Vis absorption and fluorescence emission energies to the experimental results leads us to infer that this approach allows for a proper description of the electronic excitations. In particular, implicit solvation by means of the COSMO model appears to be sufficient for the treatment of these photoacids in DMSO. The calculations confirm the presumption that a charge redistribution from the hydroxy group to the aromatic ring and the electron-withdrawing substituents is the origin of photoacidity for these photoacids. Moreover, the calculations with the continuum solvation model predict that the pK_a jump upon excitation decreases with increasing solvent polarity, as rationalized based on the Förster cycle.

Prediction of acid pKa values in the solvent acetone based on COSMO-RS

In this contribution we extent the use of the conductor-like screening model for realistic solvation (COSMO-RS) to the prediction of pK_a values in acetone, a commonly used dipolar aprotic solvent. For this, we calculated the Gibbs free energy of dissociation of 120 organic acids (nine acrylic acids, 87 benzoic acids, nine phenols, and 15 benzenesulfonamides) using COSMO-RS at the two levels BP-TZVP and BP-TZVPD-FINE and determined the parameters for a linear free energy relation for the p K a prediction by performing linear fits to experimental values. Our results suggest that the data set dissects into two groups, with the phenols being different from the other three subsets. The acrylic and benzoic acids and the sulfonamides can be treated together and yield an excellent linear correlation ( r 2 > 0.95 ) with an RMSD of only ~0.3. The slope is found to be significantly smaller than the theoretical value ( 1 / RT ln 10 ), only 45% of it, which is in accordance with findings in the literature. The phenols, however, while similarly well correlated in their own subset with an RMSD of 1.7-1.9, exhibit a slope larger than one. We discuss both a higher uncertainty in the reference values as well as physical origins as possible reasons.

How to Predict the pK a of Any Compound in Any Solvent

Acid-base properties of molecules in nonaqueous solvents are of critical importance for almost all areas of chemistry. Despite this very high relevance, our knowledge is still mostly limited to the pK a of rather few compounds in the most common solvents, and a simple yet truly general computational procedure to predict pK a's of any compound in any solvent is still missing. In this contribution, we describe such a procedure. Our method requires only the experimental pK a of a reference compound in water and a few standard quantum-chemical calculations. This method is tested through computing the proton solvation energy in 39 solvents and by comparing the pK a of 142 simple compounds in 12 solvents. Our computations indicate that the method to compute the proton solvation energy is robust with respect to the detailed computational setup and the construction of the solvation model. The unscaled pK a's computed using an implicit solvation model on the other hand differ significantly from the experimental data. These differences are partly associated with the poor quality of the experimental data and the well-known shortcomings of implicit solvation models. General linear scaling relationships to correct this error are suggested for protic and aprotic media. Using these relationships, the deviations between experiment and computations drop to a level comparable to that observed in water, which highlights the efficiency of our method.

Brønsted acidity in zeolites measured by deprotonation energy

Acid forms of zeolites have been used in industry for several decades but scaling the strength of their acid centers is still an unresolved and intensely debated issue. In this paper, the Brønsted acidity strength in aluminosilicates measured by their deprotonation energy (DPE) was investigated for FAU, CHA, IFR, MOR, FER, MFI, and TON zeolites by means of periodic and cluster calculations at the density functional theory (DFT) level. The main drawback of the periodic DFT is that it does not provide reliable absolute values due to spurious errors associated with the background charge introduced in anion energy calculations. To alleviate this problem, we employed a novel approach to cluster generation to obtain accurate values of DPE. The cluster models up to 150 T atoms for the most stable Brønsted acid sites were constructed on spheres of increasing diameter as an extension of Harrison's approach to calculating Madelung constants. The averaging of DPE for clusters generated this way provides a robust estimate of DPE for investigated zeolites despite slow convergence with the cluster size. The accuracy of the cluster approach was further improved by a scaled electrostatic embedding scheme proposed in this work. The electrostatic embedding model yields the most reliable values with the average deprotonation energy of about 1245 ± 9 kJ·mol^-1 for investigated acidic zeolites. The cluster calculations strongly indicate a correlation between the deprotonation energy and the zeolite framework density. The DPE results obtained with our electrostatic embedding model are highly consistent with the previously reported QM/MM and periodic calculations.

Surface Coverage Simulation and 3D Plotting of Main Process Parameters for Molybdenum and Vanadium Adsorption onto Ferrihydrite

Ferrihydrite, FHY, was synthesized and characterized for morphology, mineralogy, surface area, hydrodynamic diameter and surface charge properties before molybdenum (Mo) and vanadium (V) adsorption. The potentiometric titration results showed first direct evidence that CO₂ affects FHY surface sites at pH 6–9. Beside CO₂, particles concentration may affect surface properties with an impact on adsorption performance. Additional new adsorption simulation results on theoretical surface coverage vs. experimental results obtained at varying particles concentration help theoreticians and experimentalists to better estimate and apply anion adsorption processes to real environments and suggest that simulation may not always be entirely reliable. Uptake capacities obtained experimentally, varying pH, particles and metals concentrations, were plotted to assess their synergetic effect and derive trends for future process optimization. Adsorption kinetics and isotherms were also considered. Experimentally derived values for maximum uptake capacities (0.43 and 1.20 mmol g⁻¹, for Mo and V, respectively) and partitioning coefficients have applications, such as in making decisions for anions removal from wastewaters to achieve depollution efficiency or concentration required for effluents discharge and also implications in elements cycling from a geochemical perspective. In this work, the 3D plotting of the main adsorption process parameters obtained experimentally showed inter-correlations between significant process parameters that influence the adsorption process, and provides guidelines for its optimization and indicates that laboratory data can be transposed to real systems.

Computational Estimation of the Acidities of Pyrimidines and Related Compounds

Pyrimidines are key components in the genetic code of living organisms and the pyrimidine scaffold is also found in many bioactive and medicinal compounds. The acidities of these compounds, as represented by their pKas, are of special interest since they determine the species that will prevail under different pH conditions. Here, a quantum chemical quantitative structure-activity relationship (QSAR) approach was employed to estimate these acidities. Density-functional theory calculations at the B3LYP/6-31+G(d,p) level and the SM8 aqueous solvent model were employed, and the energy difference ∆EH2O between the parent compound and its dissociation product was used as a variation parameter. Excellent estimates for both the cation → neutral (pKa1, R2 = 0.965) and neutral → anion (pKa2, R2 = 0.962) dissociations were obtained. A commercial package from Advanced Chemical Design also yielded excellent results for these acidities.

Predicting pKa Values of Quinols and Related Aromatic Compounds with Multiple OH Groups

Quinonoid compounds play central roles as redox-active agents in photosynthesis and respiration and are also promising replacements for inorganic materials currently used in batteries. To design new quinonoid compounds and predict their state of protonation and redox behavior under various conditions, their pK_a values must be known. Methods that can predict the pK_a values of simple phenols cannot reliably handle complex analogues in which multiple OH groups are present and may form intramolecular hydrogen bonds. We have therefore developed a straightforward method based on a linear relationship between experimental pK_a values and calculated differences in energy between quinols and their deprotonated forms. Simple adjustments allow reliable predictions of pK_a values when intramolecular hydrogen bonds are present. Our approach has been validated by showing that predicted and experimental values for over 100 quinols and related compounds differ by an average of only 0.3 units. This accuracy makes it possible to select proper pK_a values when experimental data vary, predict the acidity of quinols and related compounds before they are made, and determine the sites and orders of deprotonation in complex structures with multiple OH groups.

DFT-Machine Learning Approach for Accurate Prediction of pKa

In this study, we propose a novel method of pK_a prediction in a diverse set of acids, which combines density functional theory (DFT) method with machine learning (ML) methods. First, the DFT method with B3LYP/6-31++G**/SM8 is used to predict pK_a, yielding a mean absolute error of 1.85 pK_a units. Subsequently, such pK_a values predicted from the DFT method are employed as one of 10 molecular descriptors for developing ML models trained on experimental data. Kernel Ridge Regression (KRR), Gaussian Process Regression, and Artificial Neural Network are optimized using three Pipelines: Pipeline 1 involving only hyperparameter optimization (HPO), Pipeline 2 involving HPO followed by a relative contribution analysis (RCA) and recursive feature elimination (RFE), and Pipeline 3 involving HPO followed by RCA and RFE on an expanded set of composite features. Finally, it is demonstrated that KRR with Pipeline 3 yields optimal pK_a prediction at an MAE of 0.60 log units. This algorithm was then utilized to predict the pK_a of 37 novel acids. The two most important features were determined to be the number of hydrogen atoms in the molecule and the degree of oxidation of the acid. The predicted pK_a values were documented for future reference.

Prediction of n-octanol/water partition coefficients and acidity constants (pKa) in the SAMPL7 blind challenge with the IEFPCM-MST model

Within the scope of SAMPL7 challenge for predicting physical properties, the Integral Equation Formalism of the Miertus-Scrocco-Tomasi (IEFPCM/MST) continuum solvation model has been used for the blind prediction of n-octanol/water partition coefficients and acidity constants of a set of 22 and 20 sulfonamide-containing compounds, respectively. The log P and pK_a were computed using the B3LPYP/6-31G(d) parametrized version of the IEFPCM/MST model. The performance of our method for partition coefficients yielded a root-mean square error of 1.03 (log P units), placing this method among the most accurate theoretical approaches in the comparison with both globally (rank 8th) and physical (rank 2nd) methods. On the other hand, the deviation between predicted and experimental pK_a values was 1.32 log units, obtaining the second best-ranked submission. Though this highlights the reliability of the IEFPCM/MST model for predicting the partitioning and the acid dissociation constant of drug-like compounds compound, the results are discussed to identify potential weaknesses and improve the performance of the method.

Evaluation of log P, pKa, and log D predictions from the SAMPL7 blind challenge

The Statistical Assessment of Modeling of Proteins and Ligands (SAMPL) challenges focuses the computational modeling community on areas in need of improvement for rational drug design. The SAMPL7 physical property challenge dealt with prediction of octanol-water partition coefficients and pK_a for 22 compounds. The dataset was composed of a series of N-acylsulfonamides and related bioisosteres. 17 research groups participated in the log P challenge, submitting 33 blind submissions total. For the pK_a challenge, 7 different groups participated, submitting 9 blind submissions in total. Overall, the accuracy of octanol-water log P predictions in the SAMPL7 challenge was lower than octanol-water log P predictions in SAMPL6, likely due to a more diverse dataset. Compared to the SAMPL6 pK_a challenge, accuracy remains unchanged in SAMPL7. Interestingly, here, though macroscopic pK_a values were often predicted with reasonable accuracy, there was dramatically more disagreement among participants as to which microscopic transitions produced these values (with methods often disagreeing even as to the sign of the free energy change associated with certain transitions), indicating far more work needs to be done on pK_a prediction methods.

Quantum Chemistry Calculations for Metabolomics

A primary goal of metabolomics studies is to fully characterize the small-molecule composition of complex biological and environmental samples. However, despite advances in analytical technologies over the past two decades, the majority of small molecules in complex samples are not readily identifiable due to the immense structural and chemical diversity present within the metabolome. Current gold-standard identification methods rely on reference libraries built using authentic chemical materials ("standards"), which are not available for most molecules. Computational quantum chemistry methods, which can be used to calculate chemical properties that are then measured by analytical platforms, offer an alternative route for building reference libraries, i.e., in silico libraries for "standards-free" identification. In this review, we cover the major roadblocks currently facing metabolomics and discuss applications where quantum chemistry calculations offer a solution. Several successful examples for nuclear magnetic resonance spectroscopy, ion mobility spectrometry, infrared spectroscopy, and mass spectrometry methods are reviewed. Finally, we consider current best practices, sources of error, and provide an outlook for quantum chemistry calculations in metabolomics studies. We expect this review will inspire researchers in the field of small-molecule identification to accelerate adoption of in silico methods for generation of reference libraries and to add quantum chemistry calculations as another tool at their disposal to characterize complex samples.

Thermodynamic prediction of proton and hydrogen atom abstraction in dehydroascorbic acid and its bicyclic form

Conformation analysis of dehydroascorbic acid and its bicyclic form was performed using the density functional theory. For the energetically preferred conformations, ionization potentials (IP) and bond dissociation enthalpies (BDE) were calculated using the B3LYP functional and 6-311++G** basis set. The effects of aqueous solution were estimated using the solvation model based on density (SMD) and the polarizable continuum model (IEF-PCM). The obtained results were compared with available experimental data for reference L-ascorbic acid (vitamin C). Our calculations indicate that the investigated bicyclic metabolic product of vitamin C can also exhibit limited radical scavenging ability due to the thermodynamically preferred dissociation of tertiary —CH bonds.

Machine learning-accelerated quantum mechanics-based atomistic simulations for industrial applications

Atomistic simulations have become an invaluable tool for industrial applications ranging from the optimization of protein-ligand interactions for drug discovery to the design of new materials for energy applications. Here we review recent advances in the use of machine learning (ML) methods for accelerated simulations based on a quantum mechanical (QM) description of the system. We show how recent progress in ML methods has dramatically extended the applicability range of conventional QM-based simulations, allowing to calculate industrially relevant properties with enhanced accuracy, at reduced computational cost, and for length and time scales that would have otherwise not been accessible. We illustrate the benefits of ML-accelerated atomistic simulations for industrial R&D processes by showcasing relevant applications from two very different areas, drug discovery (pharmaceuticals) and energy materials. Writing from the perspective of both a molecular and a materials modeling scientist, this review aims to provide a unified picture of the impact of ML-accelerated atomistic simulations on the pharmaceutical, chemical, and materials industries and gives an outlook on the exciting opportunities that could emerge in the future.

Mechanistic Insights into Co and Fe Quaterpyridine-Based CO2 Reduction Catalysts: Metal-Ligand Orbital Interaction as the Key Driving Force for Distinct Pathways

Both [Co^II(qpy)(H₂O)₂]²⁺ and [Fe^II(qpy)(H₂O)₂]²⁺ (with qpy = 2,2':6',2″:6'',2‴-quaterpyridine) are efficient homogeneous electrocatalysts and photoelectrocatalysts for the reduction of CO₂ to CO. The Co catalyst is more efficient in the electrochemical reduction, while the Fe catalyst is an excellent photoelectrocatalyst ( ACS Catal. 2018, 8, 3411-3417). This work uses density functional theory to shed light on the contrasting catalytic pathways. While both catalysts experience primarily ligand-based reductions, the second reduction in the Co catalyst is delocalized onto the metal via a metal-ligand bonding interaction, causing a spin transition and a distorted ligand framework. This orbital interaction explains the experimentally observed mild reduction potential and slow kinetics of the second reduction. The decreased hardness and doubly occupied d_z²-orbital facilitate a σ-bond with the CO₂-π* in an η¹-κC binding mode. CO₂ binding is only possible after two reductions resulting in an EEC mechanism (E = electron transfer, C = chemical reaction), and the second protonation is rate-limiting. In contrast, the Fe catalyst maintains a Lewis acidic metal center throughout the reduction process because the metal orbitals do not strongly mix with the qpy-π* orbitals. This allows binding of the activated CO₂ in an η²-binding mode. This interaction stabilizes the activated CO₂ via a π-type interaction of a Fe-t_2g orbital and the CO₂-π* and a dative bond of the oxygen lone pair. This facilitates CO₂ binding to a singly reduced catalyst resulting in an ECE mechanism. The barrier for CO₂ addition and the second protonation are higher than those for the Co catalyst and rate-limiting.

Method for the accurate prediction of electron transfer potentials using an effective absolute potential

A protocol for the accurate computation of electron transfer (ET) potentials from ab initio and density functional theory (DFT) calculations is described. The method relies on experimental pKa values, which can be measured accurately, to compute a computational setup dependent effective absolute potential. The effective absolute potentials calculated using this protocol display strong variations between the different computational setups and deviate in several cases significantly from the "generally accepted" value of 4.28 V. The most accurate estimate, obtained from CCSD(T)/aug-ccpvqz, indicates an absolute potential of 4.14 V for the normal hydrogen electrode (nhe) in water. Using the effective absolute potential in combination with CCSD(T) and a moderately sized basis, we are able to predict ET potentials accurately for a test set of small organic molecules (σ = 0.13 V). Similarly we find the effective absolute potential method to perform equally good or better for all considered DFT functionals compared to using one of the literature values for the absolute potential. For, M06-2X, which comprises the most accurate DFT method, standard deviation of 0.18 V is obtained. This improved performance is a result of using the most appropriate effective absolute potential for a given method.

pK-Yay: A Black-Box Method Using Density Functional Theory and Implicit Solvation Models to Compute Aqueous pKa Values of Weak and Strong Acids

The concept of pK_a is very important in chemistry and biology. Over the past two decades, electronic structure calculations have made huge strides toward becoming reliable counterparts to experiments in determining solution-phase pK_a values. However, invariably, the computation of pK_a values involves the use of methods whose error bars are intrinsically larger (definitely >2 kcal/mol with density functionals) than the accuracy desired for the estimation of pK_a values (<1 or 1.5 pK_a units). This scenario presents an ample scope for innovation in developing systematic error cancellation methods even today. In this work, we have developed the pK-Yay method. It is a user-friendly black-box method used to compute aqueous pK_a values of strong and weak acids employing routinely used and computationally inexpensive density functionals and implicit solvation models. It does not require the use of any explicit solvent molecule or modifying any other parameter in an electronic structure program. As part of evaluating the method, a comprehensive test set of 26 weak and strong organic and mineral acids covering 35 pK_a units (20 to -15) was assembled. The detailed description of the method, its performance with different functionals (ωB97X-D performed best for organic acids with a mean absolute error (MAE) of 0.8 pK_a units, and B2-PLYP performed best for strong acids with an MAE of 1.6 pK_a units), the strengths and limitations in the present version, and a future scope to improve the accuracy and reduce its empirical nature are presented herein.

Nucleation mechanisms and speciation of metal oxide clusters

The self-assembly mechanisms of polyoxometalates (POMs) are still a matter of discussion owing to the difficult task of identifying all the chemical species and reactions involved. We present a new computational methodology that identifies the reaction mechanism for the formation of metal-oxide clusters and provides a speciation model from first-principles and in an automated manner. As a first example, we apply our method to the formation of octamolybdate. In our model, we include variables such as pH, temperature and ionic force because they have a determining effect on driving the reaction to a specific product. Making use of graphs, we set up and solved 2.8 × 105 multi-species chemical equilibrium (MSCE) non-linear equations and found which set of reactions fitted best with the experimental data available. The agreement between computed and experimental speciation diagrams is excellent. Furthermore, we discovered a strong linear dependence between DFT and empirical formation constants, which opens the door for a systematic rescaling.

Chemical Equilibrium of Zinc Acetate Complexes in Ethanol Solution. A Theoretical Description through Thermodynamic Cycles

The Gibbs free energy of complexation between the Zn(II) species and acetate ligands, forming the [Zn(OAc)_n]^2-n complexes with n = 1, 2 in an ethanol solution, was assessed by two different theoretical protocols based on thermodynamic cycles. In both approaches, the solution phase Gibbs free energy of each reaction is computed by summing up contributions from gas phase thermochemistry calculations to solvation Gibbs free energies obtained in a hybrid fashion, i.e., each (neutral or electrically charged) solute was first solvated by explicit solvent molecules in order to capture relevant (micro) solute-solvent and/or solvent-solvent interactions and then, a continuum model calculation is performed in order to get the corresponding bulky solute-solvent contributions. For our first thermodynamic protocol, here denominated as variant 1, a set of x independent solvent molecules are used to screen each of the involved solutes, while the variant 2 strategy uses the fact that a set of solvent molecules may exist as aggregates (or molecular clusters) in the solvent macroscopic media, before the solvation process of solutes. Our selected quantum theoretical protocol was the M05-2X/6-31+G(d)/SMD level. We made a systematic exploration about the influence of several sources of errors, such as the solvent conformation, the number of solvent molecules used to screen each of the involved solutes, the coordination geometry of the metallic center before and after the complexation process, and the pertinence of using molecular geometries optimized in gas phase and in ethanol solution, for the computation of the Gibbs free energy variation regarding the two chemical reactions under study. We set an accuracy threshold equal or less than 4.0 kcal·mol^-1, with respect to the corresponding experimental records. The robustness of our thermodynamic strategies was then tested by computing the gas phase free energy contributions to the (solution phase) reaction free energies here assessed, using different density functional approximations, namely the M05-2X, BH&HLYP, PBE0, ωb97X-D and M06-2X functionals in conjunction with the larger 6-311+G(d,p) basis set.

Relative Order of Acidity among Hydroxyl Groups of Oxyluciferin and Emission Light Colors in Aqueous Solution

The magnitude of the acidity of the oxyluciferin in water in the ground and excited state is investigated, and it is found for the first time using computational approach that the enol group of the phenol-enol species is the most acidic in the ground state, but the deprotonation of the phenol of the phenol-keto form is the most favored in the excited state. The relative order of the acidity among the hydroxyl groups in the oxyluciferin is attributed to the sequence of the O-H bond lengths in the enol and phenol group of the phenol-enol form, and the phenol group of the phenol-keto species. The mechanism of determining the dominant emissive species in the excited state is proposed, and the dependence of emission light colors on the photoexcitation energy is elucidated by the high relative concentration of six chemical forms in the ground state and the absorption efficiency.

An experimental and theoretical study on the degradation of clonidine by hydroxyl and sulfate radicals

Emerging contaminants such as pharmaceuticals that cannot be completely removed by traditional biological treatments are ubiquitously present in water bodies with detected concentrations ranging from ng L^-1 to mg L^-1. Advanced oxidation technologies (AOTs) are promising, efficient, and environmentally friendly for the removal of these pharmaceuticals. In this study, we investigated the degradation kinetics of a model pharmaceutical, clonidine (CLD), via hydroxyl radical (OH) in UV/H₂O₂ and sulfate radical (SO₄^•-) in UV/peroxydisulfate (PS) systems for the first time. The second-order rate constants (k) of protonated cationic CLD with OH and SO₄^•- were measured to be (2.15 ± 0.07) × 10⁹ M^-1 s^-1 and (1.12 ± 0.03) × 10⁹ M^-1 s^-1, respectively. We also calculated the pK_a value of CLD and thermodynamic behaviors for reactions of CLD/HCLD⁺ with OH and SO₄^•- at M05-2X/6-311++G**//M05-2X/6-31+G** level with SMD solvation model. The pK_a value was calculated to be 8.14, confirming the literature value. H atom abstraction pathway was the most favorable pathway for both OH and SO₄^•-, while single electron transfer pathway was thermodynamically feasible only for SO₄^•- for CLD but not for HCLD⁺. In addition, the reactivities of both tautomeric forms of CLD (i.e., amino and imino CLD) with both radicals were also investigated. This study contributed to a better understanding on the degradation mechanisms of CLD and proposed the possibilities of the elimination of pharmaceuticals by applying AOTs during wastewater treatment processes.

A theoretical study on the pKa values of selenium compounds in aqueous solution

The pK_a values of different kinds of selenium compounds (R-SeH) were investigated by using the ωB97XD method with a SMD model.

Synthesis and characterization of bis(4-amino-2-bromo-6-methoxy)azobenzene derivatives

Aminoazobenzene derivatives with four ortho substituents with respect to the N–N double bond are a relatively unexplored class of azo compounds that show promise for use as photoswitches in biology. Tetra-ortho-methoxy-substituted aminoazobenzene compounds in particular can form azonium ions under physiological conditions and exhibit red-light photoswitching. Here, we report the synthesis and characterization of two bis(4-amino-2-bromo-6-methoxy)azobenzene derivatives. These compounds form red-light-absorbing azonium ions, but only under very acidic conditions (pH < 1). While the low pK_a makes the azonium form unsuitable, the neutral versions of these compounds undergo trans-to-cis photoisomerization with blue-green light and exhibit slow (τ_1/2 ≈ 10 min) thermal reversion and so may find applications under physiological conditions.

Combining proton and silaborane-based superhalogen anions - an effective route to new superacids as verified via systematic DFT calculations

Based on systematic DFT calculations, silaborane-based superhalogen anions, which obey the Wade-Mingos rule, are shown to be capable of giving rise to superacids via their combination with protons. Compared to previous carborane-based systems, the acidities of the composites here are stronger in both the gas phase and solution phase. Thus, the potential of candidates based on silaborane could be greater than those based on carborane in the search for ultra-strong acidic systems. Within a given group, a higher superhalogen anion vertical electron detachment energy (VDE) generally leads to stronger acidity. This consistency arises from the dominant role of the VDE, as established through the decomposition of the gas-phase acidity into different contributions. Thus, constructing superacids from superhalogens is a rational route whose future should be positive. Besides the VDE, other effects, i.e., the deformation energy (DE) and bond dissociation energy (BDE), could also be crucial, especially in terms of the differences between the acidities of composites belonging to different groups. A comparison between the results in the gas phase and solution phase indicates that complete calculations of both gas-phase ΔG_acid and solution-phase pK_a values are necessary to obtain an unbiased description of the acidity. The solvation free energies of the participants in the deprotonation process, especially the conjugate acid, are responsible for the discrepancies between gas phase and solution phase behavior.