Ab Initio Quantum Mechanical Calculations Research Articles

Binding Free Energy Calculation Based on the Fragment Molecular Orbital Method and Its Application in Designing Novel SHP-2 Allosteric Inhibitors.

The accurate prediction of binding free energy is a major challenge in structure-based drug design. Quantum mechanics (QM)-based approaches show promising potential in predicting ligand-protein binding affinity by accurately describing the behavior and structure of electrons. However, traditional QM calculations face computational limitations, hindering their practical application in drug design. Nevertheless, the fragment molecular orbital (FMO) method has gained widespread application in drug design due to its ability to reduce computational costs and achieve efficient ab initio QM calculations. Although the FMO method has demonstrated its reliability in calculating the gas phase potential energy, the binding of proteins and ligands also involves other contributing energy terms, such as solvent effects, the 'deformation energy' of a ligand's bioactive conformations, and entropy. Particularly in cases involving ionized fragments, the calculation of solvation free energy becomes particularly crucial. We conducted an evaluation of some previously reported implicit solvent methods on the same data set to assess their potential for improving the performance of the FMO method. Herein, we develop a new QM-based binding free energy calculation method called FMOScore, which enhances the performance of the FMO method. The FMOScore method incorporates linear fitting of various terms, including gas-phase potential energy, deformation energy, and solvation free energy. Compared to other widely used traditional prediction methods such as FEP+, MM/PBSA, MM/GBSA, and Autodock vina, FMOScore showed good performance in prediction accuracies. By constructing a retrospective case study, it was observed that incorporating calculations for solvation free energy and deformation energy can further enhance the precision of FMO predictions for binding affinity. Furthermore, using FMOScore-guided lead optimization against Src homology-2-containing protein tyrosine phosphatase 2 (SHP-2), we discovered a novel and potent allosteric SHP-2 inhibitor (compound 8).

Investigation of the effects of solvent on oxygen evolution reactions on the surface of magnesium oxide

In this research, due to the importance of oxygen evolution reactions (OER) in the production of hydrogen gas, the effects of the solvent on the OER steps on the magnesium oxide surface have been studied. For this purpose, the difference in free energy value in OER steps on the Magnesium oxide surface for two systems (in the presence and the absence of solvent) have been investigated using molecular dynamics (MD) simulation and thermodynamic integration (TI). These results are in perfect agreement with the results obtained from ab initio molecular dynamics simulations and quantum mechanics calculations. Our findings show that the presence of water solvent around the surface of magnesium oxide has a crucial role on OER and leads to an increase more than twofold in the free energy of all steps. Also, the presence of the solvent has the most effect on the third step of OER, and its ΔG increases by about 1.9 eV. This fact can be attributed to the approach of the water molecules to the substrate spontaneously, which increases the probability of forming hydrogen bonds and the number of contacts, leading to more favorable thermodynamic reactions. Close inspection of the calculated binding energies between the substrate and the intermediates confirms that the binding energy of OER's second step is significantly higher than the other steps and is the velocity determining step.

Automated, Consistent, and Even-Handed Selection of Active Orbital Spaces for Quantum Embedding.

A widely used strategy to reduce the computational cost of quantum-chemical calculations is to partition the system into an active subsystem, which is the focus of the computational efforts, and an environment that is treated at a lower computational level. The system partitioning is mostly based on localized molecular orbitals. When reaction paths or energy differences are to be calculated, it is crucial to keep the orbital space consistent for all structures. Inconsistencies in orbital space can lead to unpredictable errors on the potential energy surface. While successful strategies to ensure this consistency have been established for organic and even metal-organic systems, these methods often fail for metal clusters or nanoparticles with a high density of near-degenerate and delocalized molecular orbitals. However, such systems are highly relevant for catalysis. Accurate yet feasible quantum-mechanical ab initio calculations are therefore highly desired. In this work, we present an approach based on the subsystem projected atomic orbital decomposition algorithm that allows us to ensure automated and consistent partitioning even for systems with delocalized and near-degenerate molecular orbitals and demonstrate the validity of this method for the binding energies of small molecules on transition-metal clusters.

Melting curves of ice polymorphs in the vicinity of the liquid-liquid critical point.

The possible existence of a liquid-liquid critical point in deeply supercooled water has been a subject of debate due to the challenges associated with providing definitive experimental evidence. The pioneering work by Mishima and Stanley [Nature 392, 164-168 (1998)] sought to shed light on this problem by studying the melting curves of different ice polymorphs and their metastable continuation in the vicinity of the expected liquid-liquid transition and its associated critical point. Based on the continuous or discontinuous changes in the slope of the melting curves, Mishima [Phys. Rev. Lett. 85, 334 (2000)] suggested that the liquid-liquid critical point lies between the melting curves of ice III and ice V. We explore this conjecture using molecular dynamics simulations with a machine learning model based on abinitio quantum-mechanical calculations. We study the melting curves of ices III, IV, V, VI, and XIII and find that all of them are supercritical and do not intersect the liquid-liquid transition locus. We also find a pronounced, yet continuous, change in the slope of the melting lines upon crossing of the liquid locus of maximum compressibility. Finally, we analyze the literature in light of our findings and conclude that the scenario in which the melting curves are supercritical is favored by the most recent computational and experimental evidence. Although the preponderance of evidence is consistent with the existence of a second critical point in water, the behavior of ice polymorph melting lines does not provide strong evidence in support of this viewpoint, according to our calculations.

Femtosecond proton transfer in urea solutions probed by X-ray spectroscopy.

Proton transfer is one of the most fundamental events in aqueous-phase chemistry and an emblematic case of coupled ultrafast electronic and structural dynamics1,2. Disentangling electronic and nuclear dynamics on the femtosecond timescales remains a formidable challenge, especially in the liquid phase, the natural environment of biochemical processes. Here we exploit the unique features of table-top water-window X-ray absorption spectroscopy3-6 to reveal femtosecond proton-transfer dynamics in ionized urea dimers in aqueous solution. Harnessing the element specificity and the site selectivity of X-ray absorption spectroscopy with the aid of ab initio quantum-mechanical and molecular-mechanics calculations, we show how, in addition to the proton transfer, the subsequent rearrangement of the urea dimer and the associated change of the electronic structure can be identified with site selectivity. These results establish the considerable potential of flat-jet, table-top X-ray absorption spectroscopy7,8 in elucidating solution-phase ultrafast dynamics in biomolecular systems.

Development of AMOEBA Polarizable Force Field for Rare-Earth La3+ Interaction with Bioinspired Ligands.

Rare-earth metals (REMs) are crucial for many important industries, such as power generation and storage, in addition to cancer treatment and medical imaging. One promising new REM refinement approach involves mimicking the highly selective and efficient binding of REMs observed in relatively recently discovered proteins. However, realizing any such bioinspired approach requires an understanding of the biological recognition mechanisms. Here, we developed a new classical polarizable force field based on the AMOEBA framework for modeling a lanthanum ion (La3+) interacting with water, acetate, and acetamide, which have been found to coordinate the ion in proteins. The parameters were derived by comparing to high-level ab initio quantum mechanical (QM) calculations that include relativistic effects. The AMOEBA model, with advanced atomic multipoles and electronic polarization, is successful in capturing both the QM distance-dependent La3+-ligand interaction energies and experimental hydration free energy. A new scheme for pairwise polarization damping (POLPAIR) was developed to describe the polarization energy in La3+ interactions with both charged and neutral ligands. Simulations of La3+ in water showed water coordination numbers and ion-water distances consistent with previous experimental and theoretical findings. Water residence time analysis revealed both fast and slow kinetics in water exchange around the ion. This new model will allow investigation of fully solvated lanthanum ion-protein systems using GPU-accelerated dynamics simulations to gain insights on binding selectivity, which may be applied to the design of synthetic analogues.

A Density Functional Theory for the Average Electron Energy.

A formally exact density functional theory (DFT) determination of the average electron energy is presented. Our theory, which is based on a different accounting of energy functional terms, partially solves one well-known downside of conventional Kohn-Sham (KS) DFT: that electronic energies have but tenuous connections to physical quantities. Calculated average electron energies are close to experimental ionization potentials (IPs) in one-electron systems, demonstrating a surprisingly small effect of self-interaction and other exchange-correlation errors in established DFT methods. Remarkable agreement with ab initio quantum mechanical calculations of multielectron systems is demonstrated using several flavors of DFT, and we argue for the use of the average electron energy as a design criterion for density functional approximations.

Development of a ReaxFF potential for Au–Pd

The bimetallic alloys often outperform their single-component counterparts due to synergistic effects. Being widely known, the Au–Pd alloy is a promising candidate for the novel heterogeneous nanocatalysts. Rational design of such systems requires theoretical simulations under ambient conditions. Ab initio quantum-mechanical calculations employ the density functional theory (DFT) and are limited to the systems with few tens of atoms and short timescales. The alternative solution implies development of reliable atomistic potentials. Among different approaches ReaxFF combines chemical accuracy and low computational costs. However, the development of a new potential is a problem without unique solution and thus requires accurate validation criteria. In this work we construct ReaxFF potential for the Au–Pd system based on ab initio DFT calculations for bulk structures, slabs and nanoparticles with different stoichiometry. The validation was performed with molecular dynamics and Monte-Carlo calculations. We present several optimal parametrizations that describe experimental bulk mechanical and thermal properties, atomic order-disorder phase transition temperatures and the resulting ordered crystal structures.

Protein-Ligand Binding Free-Energy Calculations with ARROW─A Purely First-Principles Parameterized Polarizable Force Field.

Protein-ligand binding free-energy calculations using molecular dynamics (MD) simulations have emerged as a powerful tool for in silico drug design. Here, we present results obtained with the ARROW force field (FF)─a multipolar polarizable and physics-based model with all parameters fitted entirely to high-level ab initio quantum mechanical (QM) calculations. ARROW has already proven its ability to determine solvation free energy of arbitrary neutral compounds with unprecedented accuracy. The ARROW FF parameterization is now extended to include coverage of all amino acids including charged groups, allowing molecular simulations of a series of protein-ligand systems and prediction of their relative binding free energies. We ensure adequate sampling by applying a novel technique that is based on coupling the Hamiltonian Replica exchange (HREX) with a conformation reservoir generated via potential softening and nonequilibrium MD. ARROW provides predictions with near chemical accuracy (mean absolute error of ∼0.5 kcal/mol) for two of the three protein systems studied here (MCL1 and Thrombin). The third protein system (CDK2) reveals the difficulty in accurately describing dimer interaction energies involving polar and charged species. Overall, for all of the three protein systems studied here, ARROW FF predicts relative binding free energies of ligands with a similar accuracy level as leading nonpolarizable force fields.

Investigations of an Unexpected [2+2] Photocycloaddition in the Synthesis of (-)-Scabrolide A from Quantum Mechanics Calculations.

We utilize ab initio quantum mechanics calculations to evaluate a range of plausible mechanistic pathways for the unexpected formation of a [6-4-4] ring system from an enone-olefin photocycloaddition in the synthesis of (-)-scabrolide A, previously reported by our group. We present a mechanistic analysis that is consistent with all current experimental observations, including the photoexcitation, the C-C bond formation, and the associated chemo- and diastereoselectivity.

Theoretical Investigation of 1H and 13C NMR Spectra of Diethanol Amine Dithiocarbamate RAFT Agent

In this study, quantum chemical studies of diethanolamine dithiocarbamate (DADC) RAFT compound were carried out. For this purpose, the composite was optimized by using the DFT / B3LYP method 6-311G and the B3PW91 method SDD basis set. By using the optimized structure obtained according to the GIAO method, 1H and 13C NMR chemical shift values in the gas phase were calculated. According to the results obtained, it was seen that the theoretical data were coherent with the experimental data. In addition, in the theoretical part of the study, the FT-IR frequency values of the compound synthesized by using the same methods and basic set were compared experimentally and theoretically. In addition, the structure details and analysis of molecules, electronic properties such as HOMO and LUMO energies, molecular electrostatic potential (MEP) and thermodynamic properties have been performed. The electric dipole moment (μ) and the initial hyperpolarizability (β) values of the studied molecule were calculated by using ab initio quantum mechanical calculations. In addition, the calculated results show that the (DADC) molecule can have nonlinear optical (NLO) behavior with nonzero values. Neutral Band Orbital (NBO) analysis has been calculated with DFT / B3PW91/SDD basis set.

Exploring Deep Learning of Quantum Chemical Properties for Absorption, Distribution, Metabolism, and Excretion Predictions.

Quantum mechanical (QM) descriptors of small molecules have wide applicability in understanding organic reactivity and molecular properties, but the substantial compute cost required for ab initio QM calculations limits their broad usage. Here, we investigate the use of deep learning for predicting QM descriptors, with the goal of enabling usage of near-QM accuracy electronic properties on large molecular data sets such as those seen in drug discovery. Several deep learning approaches have previously been benchmarked on a published data set called QM9, where 12 ground-state properties have been calculated for molecules with up to nine heavy atoms, limited to C, H, N, O, and F elements. To advance the work beyond the QM9 chemical space and enable application to molecules encountered in drug discovery, we extend the QM9 data set by creating a QM9-extended data set covering an additional ∼20,000 molecules containing S and Cl atoms. Using this extended set, we generate new deep learning models as well as leverage ANI-2x models to provide predictions on larger, more diverse molecules common in drug discovery, and we find the models estimate 11 of 12 ground-state properties reasonably. We use the predicted QM descriptors to augment graph convolutional neural network (GCNN) models for selected ADME end points (rat microsomal clearance, hepatic clearance, total clearance, and P-glycoprotein efflux) and found varying degrees of performance improvement compared to nonaugmented GCNN models, including pronounced improvement in P-glycoprotein efflux prediction.

Calculations of pKa Values for a Series of Naturally Occurring Modified Nucleobases.

Modified nucleobases are found in functionally important regions of RNA and are often responsible for essential structural roles. Many of these nucleobase modifications are dynamically regulated in nature, with each modification having a different biological role in RNA. Despite the high abundance of modifications, many of their characteristics are still poorly understood. One important property of a nucleobase is its pKa value, which has been widely studied for unmodified nucleobases, but not for the modified versions. In this study, the pKa values of modified nucleobases were determined by performing ab initio quantum mechanical calculations using a B3LYP density functional with the 6-31+G(d,p) basis set and a combination of implicit-explicit solvation systems. This method, which was previously employed to determine the pKa values of unmodified nucleobases, is applicable to a variety of modified nucleobases. Comparisons of the pKa values of modified nucleobases give insight into their structural and energetic impacts within nucleic acids.

Seeding the multi-dimensional nonequilibrium pulling for Hamiltonian variation: indirect nonequilibrium free energy simulations at QM levels.

The combination of free energy simulations in the alchemical and configurational spaces provides a feasible route to access the thermodynamic profiles under a computationally demanding target Hamiltonian. Normally, due to the significant differences between the computational cost of ab initio quantum mechanics (QM) calculations and those of semi-empirical quantum mechanics (SQM) and molecular mechanics (MM), this indirect method could be used to obtain the QM thermodynamics by combining the SQM or MM results and the SQM-to-QM or MM-to-QM corrections. In our previous work, a multi-dimensional nonequilibrium pulling framework for Hamiltonian variations was introduced based on bidirectional pulling and bidirectional reweighting. The method performs nonequilibrium free energy simulations in the configurational space to obtain the thermodynamic profile along the conformational change pathway under a selected computationally efficient Hamiltonian, and uses the nonequilibrium alchemical method to correct or perturb the thermodynamic profile to that under the target Hamiltonian. The BAR-based method is designed to achieve the best generality and transferability and thus leads to modest (∼20 fold) speedup. In this work, we explore the possibility of further accelerating the nonequilibrium free energy simulation by employing unidirectional pulling and using the selection criterion to obtain the initial configurations used to initiate nonequilibrium trajectories following the idea of adaptive steered molecular dynamics (ASMD). A single initial condition is used to seed the whole multi-dimensional nonequilibrium free energy simulation and the sampling is performed fully in the nonequilibrium ensemble. Introducing very short ps-length equilibrium sampling to grab more initial seeds could also be helpful. The ASMD scheme estimates the free energy difference with the unidirectional exponential average (EXP), but it does not follow exactly the requirements of the EXP estimator. Another deficiency of the seeding simulation is the inherently sequential or serial pulling due to the inter-segment dependency, which triggers some problems in the parallelizability of the simulation. Numerical tests are performed to grasp some insights and guidelines for using this selection-criterion-based ASMD scheme. The presented selection-criterion-based multi-dimensional ASMD scheme follows the same perturbation network of the BAR-based method, and thus could be used in various Hamiltonian-variation cases.

First-Principles Simulation of Dielectric Function in Biomolecules.

The dielectric spectra of complex biomolecules reflect the molecular heterogeneity of the proteins and are particularly important for the calculations of electrostatic (Coulomb) and electrodynamic (van der Waals) interactions in protein physics. The dielectric response of the proteins can be decomposed into different components depending on the size, structure, composition, locality, and environment of the protein in general. We present a new robust simulation method anchored in rigorous ab initio quantum mechanical calculations of explicit atomistic models, without any indeterminate parameters to compute and gain insight into the dielectric spectra of small proteins under different conditions. We implement this methodology to a polypeptide RGD-4C (1FUV) in different environments, and the SD1 domain in the spike protein of SARS-COV-2. Two peaks at 5.2–5.7 eV and 14.4–15.2 eV in the dielectric absorption spectra are observed for 1FUV and SD1 in vacuum as well as in their solvated and salted models.

Magnetic Anisotropy in a Cubane-like Ni4 Complex: An Ab Initio Perspective.

Magnetic anisotropy, in the absence of an external magnetic field, relates to the degeneracy lift of energy levels. In the standard case of transition metal complexes, this property is usually modeled by an anisotropic spin Hamiltonian and one speaks of "zero-field splitting" (ZFS) of spin states. While the case of mononuclear complexes has been extensively described by means of ab initio quantum mechanical calculations, the literature on polynuclear complexes studied with these methodologies is rather scarce. In this work, advanced multiconfigurational wave function theory methods are applied to compute the ZFS of the ground S = 4 state of an actual tetranickel(II) complex, displaying a magnet behavior below 0.5 K. First, the isotropic couplings are computed in the absence of the spin-orbit coupling operator, in the full complex and also in clusters with only two active nickel(II) centers, confirming the occurrence of weak ferromagnetic couplings in this system. Second, the single-site magnetic anisotropies are computed on a cluster bearing only one active nickel(II) site, showing that the single-site anisotropy axes are not oriented in an optimal fashion for generating a large uniaxial molecular anisotropy. Furthermore, the possibility for involving only a few local orbital excited states in the calculation is assessed, actually opening the way for a consistent and manageable treatment of the ZFS of the ground S = 4 state. Third, multiconfigurational calculations are performed on the full complex, confirming the weak uniaxial anisotropy occurring for this state and also, interestingly, revealing a significant contribution of the lowest-lying orbitally excited S = 3 states. Overall, by comparison with the experiment, the reported results question the common habit of using only one structure, in particular derived from a crystallography experiment, to compute magnetic anisotropy parameters.

Solid-liquid coexistence of neon, argon, krypton, and xenon studied by simulations.

The noble elements constitute the simplest group of atoms. At low temperatures or high pressures, they freeze into the face-centered cubic (fcc) crystal structure (except helium). This paper investigates neon, argon, krypton, and xenon by molecular dynamics using the simplified atomic potentials recently proposed by Deiters and Sadus [J. Chem. Phys. 150, 134504 (2019)], which are parameterized using data from accurate ab initio quantum-mechanical calculations by the coupled-cluster approach at the single-double-triple level. We compute the fcc freezing lines and find good agreement with the empirical values. At low pressures, predictions are improved by including many-body corrections. Hidden scale invariance of the potential-energy function is established by showing that mean-squared displacement and the static structure factor are invariant along the lines of constant excess entropy (isomorphs). The isomorph theory of melting [Pedersen et al., Nat. Commun. 7, 12386 (2016)] is used to predict from simulations at a single state point the freezing line's shape, the entropy of melting, and the Lindemann parameter of the crystal at melting. Finally, our results suggest that the body-centered cubic crystal is the thermodynamically stable phase at high pressures.

Determining the atomic charge of calcium ion requires the information of its coordination geometry in an EF-hand motif.

It is challenging to parameterize the force field for calcium ions (Ca2+) in calcium-binding proteins because of their unique coordination chemistry that involves the surrounding atoms required for stability. In this work, we observed a wide variation in Ca2+ binding loop conformations of the Ca2+-binding protein calmodulin, which adopts the most populated ternary structures determined from the molecular dynamics simulations, followed by ab initio quantum mechanical (QM) calculations on all 12 amino acids in the loop that coordinate Ca2+ in aqueous solution. Ca2+ charges were derived by fitting to the electrostatic potential in the context of a classical or polarizable force field (PFF). We discovered that the atomic radius of Ca2+ in conventional force fields is too large for the QM calculation to capture the variation in the coordination geometry of Ca2+ in its ionic form, leading to unphysical charges. Specifically, we found that the fitted atomic charges of Ca2+ in the context of PFF depend on the coordinating geometry of electronegative atoms from the amino acids in the loop. Although nearby water molecules do not influence the atomic charge of Ca2+, they are crucial for compensating for the coordination of Ca2+ due to the conformational flexibility in the EF-hand loop. Our method advances the development of force fields for metal ions and protein binding sites in dynamic environments.

In vivo uptake of antisense oligonucleotide drugs predicted by ab initio quantum mechanical calculations

Liver and kidney uptake and antisense activity is studied for a series of Locked Nucleic Acid (LNA) oligonucleotides with fully stereo-defined, internucleoside linkages. These stereo-specific phosphorothioates are made with a newly developed synthesis method and are being analyzed both theoretically and experimentally. Their structures are obtained theoretically by using many-body Schrödinger equations applied to a group of 11 stereo-defined LNA antisense oligonucleotides selected for biological experiments. The fully converged electronic structures were obtained from ab initio quantum calculations providing the specific electronic structures. One important result was the observation that the calculated electronic structure, represented by the iso-surface area of the electron density in Å2, correlated linearly with LNA oligonucleotide uptake in the liver and kidney. This study also shows that more complex biological phenomena, such as drug activity, will require more molecular and cellular identifiers than used here before a correlation can be found. Establishing biological correlations between quantum mechanical (QM) calculated structures and antisense oligonucleotides is novel, and this method may constitute new tools in drug discovery.

Computation of host-guest binding free energies with a new quantum mechanics based mining minima algorithm.

A new method called QM-VM2 is presented that efficiently combines statistical mechanics with quantum mechanical (QM) energy potentials in order to calculate noncovalent binding free energies of host-guest systems. QM-VM2 efficiently couples the use of semi-empirical QM (SEQM) energies and geometry optimizations with an underlying molecular mechanics (MM) based conformational search, to find low SEQM energy minima, and allows for processing of these minima at higher levels of ab initio QM theory. A progressive geometry optimization scheme is introduced as a means to increase conformational sampling efficiency. The newly implemented QM-VM2 is used to compute the binding free energies of the host molecule cucurbit[7]uril and a set of 15 guest molecules. The results are presented along with comparisons to experimentally determined binding affinities. For the full set of 15 host-guest complexes, which have a range of formal charges from +1 to +3, SEQM-VM2 based binding free energies show poor correlation with experiment, whereas for the ten +1 complexes only, a significant correlation (R2 = 0.8) is achieved. SEQM-VM2 generation of conformers followed by single-point ab initio QM calculations at the dispersion corrected restricted Hartree-Fock-D3(BJ) and TPSS-D3(BJ) levels of theory, as post-processing corrections, yields a reasonable correlation with experiment for the full set of host-guest complexes (R2 = 0.6 and R2 = 0.7, respectively) and an excellent correlation for the +1 formal charge set (R2 = 1.0 and R2 = 0.9, respectively), as long as a sufficiently large basis set (triple-zeta quality) is employed. The importance of the inclusion of configurational entropy, even at the MM level, for the achievement of good correlation with experiment was demonstrated by comparing the calculated ΔE values with experiment and finding a considerably poorer correlation with experiment than for the calculated free energy ΔE - TΔS. For the complete set of host-guest systems with the range of formal charges, it was observed that the deviation of the predicted binding free energy from experiment correlates somewhat with the net charge of the systems. This observation leads to a simple empirical interpolation scheme to improve the linear regression of the full set.