Search

Accurate elucidation of gas-phase chemical structures using collision cross section (CCS) values obtained from ion-mobility mass spectrometry benefits from a synergism between experimental and in silico results. We have shown in recent work that for a molecule of modest size with a proscribed conformational space we can successfully capture a conformation(s) that can match experimental CCS values. However, for flexible systems such as fatty acids that have many rotatable bonds and multiple intramolecular London dispersion interactions, it becomes necessary to sample a much greater conformational space. Sampling more conformers, however, accrues significant computational cost downstream in optimization steps involving quantum mechanics. To reduce this computational expense for lipids, we have developed a novel machine learning (ML) model to facilitate conformer filtering according to the estimated gas-phase CCS values. Herein we report that the implementation of our CCS knowledge-based approach for conformational sampling resulted in improved structure prediction agreement with experiment by achieving favorable average CCS prediction errors of ∼2% for lipid systems in both the validation set and the test set. Moreover, most of the gas-phase candidate conformations obtained by using CCS focusing achieved lower energy-minimum geometries than the candidate conformations without focusing. Altogether, the implementation of this ML model into our modeling workflow has proven to be beneficial for both the quality of the results and the turnaround time. Finally, while our approach is limited to lipids, it can be readily extended to other molecules of interest.

The Zrt/Irt-like protein (ZIP) metal transporters are key players in maintaining the homeostasis of a panel of essential microelements. The prototypical ZIP from Bordetella bronchiseptica (BbZIP) is an elevator transporter, but how the metal substrate moves along the transport pathway and how the transporter changes conformation to allow alternating access remain to be elucidated. Here, we combined structural, biochemical, and computational approaches to investigate the process of metal substrate translocation along with the global structural rearrangement. Our study revealed an upward hinge motion of the transport domain in a high-resolution crystal structure of a cross-linked variant, elucidated the mechanisms of metal release from the transport site into the cytoplasm and activity regulation by a cytoplasmic metal-binding loop, and unraveled an unusual elevator mode in enhanced sampling simulations that distinguishes BbZIP from other elevator transporters. This work provides important insights into the metal transport mechanism of the ZIP family., Competing Interests: Conflicting Interests The authors declare that they have no competing interests.

The design of transition-metal complexes (TMCs) has drawn much attention over the years because of their important applications as metallodrugs and functional materials. In this work, we present an extension of our recently reported approach, LigandDiff [Jin et al. J. Chem. Theory Comput. 20 , 4377(2024)]. The new model, which we call multi-LigandDiff, is more flexible and greatly outperforms its predecessor. This scaffold-based diffusion model allows de novo ligand design with either existing ligands or without any ligand. Moreover, it allows users to predefine the denticity of the generated ligand. Our results indicate that multi-LigandDiff can generate well-defined ligands and is transferable to multiple transition metals and coordination geometries. In terms of its application, multi-LigandDiff successfully designed 338 Fe(II) spin-crossover (SCO) complexes from only 47 experimentally validated SCO complexes. And these generated complexes are configurationally diverse and structurally reasonable. Overall, the results show that multi-LigandDiff is an ideal tool to design novel TMCs from scratch.

One commonly observed binding motif in metalloproteins involves the interaction between a metal ion and histidine's imidazole side chains. Although previous imidazole-M(II) parameters established the flexibility and reliability of the 12-6-4 Lennard-Jones (LJ)-type nonbonded model by simply tuning the ligating atom's polarizability, they have not been applied to multiple-imidazole complexes. To fill this gap, we systematically simulate multiple-imidazole complexes (ranging from one to six) for five metal ions (Co(II), Cu(II), Mn(II), Ni(II), and Zn(II)) which commonly appear in metalloproteins. Using extensive (40 ns per PMF window) sampling to assemble free energy association profiles (using OPC water and standard HID imidazole charge models from AMBER) and comparing the equilibrium distances to DFT calculations, a new set of parameters was developed to focus on energetic and geometric features of multiple-imidazole complexes. The obtained free energy profiles agree with the experimental binding free energy and DFT calculated distances. To validate our model, we show that we can close the thermodynamic cycle for metal-imidazole complexes with up to six imidazole molecules in the first solvation shell. Given the success in closing the thermodynamic cycles, we then used the same extended sampling method for six other metal ions (Ag(I), Ca(II), Cd(II), Cu(I), Fe(II), and Mg(II)) to obtain new parameters. Since these new parameters can reproduce the one-imidazole geometry and energy accurately, we hypothesize that they will reasonably predict the binding free energy of higher-level coordination numbers. Hence, we did not extend the analysis of these ions up to six imidazole complexes. Overall, the results shed light on metal-protein interactions by emphasizing the importance of ligand-ligand interaction and metal-π-stacking within metalloproteins.

We report the development and testing of new integrated cyberinfrastructure for performing free energy simulations with generalized hybrid quantum mechanical/molecular mechanical (QM/MM) and machine learning potentials (MLPs) in Amber. The Sander molecular dynamics program has been extended to leverage fast, density-functional tight-binding models implemented in the DFTB+ and xTB packages, and an interface to the DeePMD-kit software enables the use of MLPs. The software is integrated through application program interfaces that circumvent the need to perform "system calls" and enable the incorporation of long-range Ewald electrostatics into the external software's self-consistent field procedure. The infrastructure provides access to QM/MM models that may serve as the foundation for QM/MM-ΔMLP potentials, which supplement the semiempirical QM/MM model with a MLP correction trained to reproduce ab initio QM/MM energies and forces. Efficient optimization of minimum free energy pathways is enabled through a new surface-accelerated finite-temperature string method implemented in the FE-ToolKit package. Furthermore, we interfaced Sander with the i-PI software by implementing the socket communication protocol used in the i-PI client-server model. The new interface with i-PI allows for the treatment of nuclear quantum effects with semiempirical QM/MM-ΔMLP models. The modular interoperable software is demonstrated on proton transfer reactions in guanine-thymine mispairs in a B-form deoxyribonucleic acid helix. The current work represents a considerable advance in the development of modular software for performing free energy simulations of chemical reactions that are important in a wide range of applications., (© 2024 Author(s). Published under an exclusive license by AIP Publishing.)

Transition metal complexes are a class of compounds with varied and versatile properties, making them of great technological importance. Their applications cover a wide range of fields, either as metallodrugs in medicine or as materials, catalysts, batteries, solar cells, etc. The demand for the novel design of transition metal complexes with new properties remains of great interest. However, the traditional high-throughput screening approach is inherently expensive and laborious since it depends on human expertise. Here, we present LigandDiff, a generative model for the de novo design of novel transition metal complexes. Unlike the existing methods that simply extract and combine ligands with the metal to get new complexes, LigandDiff aims at designing configurationally novel ligands from scratch, which opens new pathways for the discovery of organometallic complexes. Moreover, it overcomes the limitations of current methods, where the diversity of new complexes highly relies on the diversity of available ligands, while LigandDiff can design numerous novel ligands without human intervention. Our results indicate that LigandDiff designs unique and novel ligands under different contexts, and these generated ligands are synthetically accessible. Moreover, LigandDiff shows good transferability by generating successful ligands for any transition metal complex.

Phosphate derivatives and their interaction with metal cations are involved in many important biological phenomena, so an accurate characterization of the phosphate-metal interaction is necessary to properly understand the role of phosphate-metal contacts in mediating biological function. Herein, we improved the standard 12-6 Lennard-Jones (LJ) potential via the usage of the 12-6-4 LJ model, which incorporates ion-induced dipole interactions. Via parameter scanning, we fine-tuned the 12-6-4 LJ polarizability values to obtain accurate absolute binding free energies for the phosphate anions H 2 PO 4 - , HPO 4 2- , PO 4 3- coordinating with Ca 2+ and Mg 2+ . First, we modified the phosphate 12-6-4 LJ parameters to reproduce the solvation free energies of the series of phosphate anions using the thermodynamic integration (TI) method. Then, using the potential mean force (PMF) method, the polarizability of the metal-phosphate interaction was obtained. We show that the free energy profiles of phosphate ions coordinated to Ca 2+ and Mg 2+ generally show similar trends at longer metal-phosphate distances, while the absolute binding energy values increased with deprotonation. The resulting parameters demonstrate the flexibility of the 12-6-4 LJ-type nonbonded model and its usefulness in accurately describing cation-anion interactions.

Understanding the energetic landscapes of large molecules is necessary for the study of chemical and biological systems. Recently, deep learning has greatly accelerated the development of models based on quantum chemistry, making it possible to build potential energy surfaces and explore chemical space. However, most of this work has focused on organic molecules due to the simplicity of their electronic structures as well as the availability of data sets. In this work, we build a deep learning architecture to model the energetics of zinc organometallic complexes. To achieve this, we have compiled a configurationally and conformationally diverse data set of zinc complexes using metadynamics to overcome the limitations of traditional sampling methods. In terms of the neural network potentials, our results indicate that for zinc complexes, partial charges play an important role in modeling the long-range interactions with a neural network. Our developed model outperforms semiempirical methods in predicting the relative energy of zinc conformers, yielding a mean absolute error (MAE) of 1.32 kcal/mol with reference to the double-hybrid PWPB95 method.

We report a Fe(II) data set of more than 23000 conformers in both low-spin (LS) and high-spin (HS) states. This data set was generated to develop a neural network model that is capable of predicting the energy and the energy splitting as a function of the conformation of a Fe(II) organometallic complex. In order to achieve this, we propose a type of scaled electronic embedding to cover the long-range interactions implicitly in our neural network describing the Fe(II) organometallic complexes. For the total energy prediction, the lowest MAE is 0.037 eV, while the lowest MAE of the splitting energy is 0.030 eV. Compared to baseline models, which only incorporate short-range interactions, our scaled electronic embeddings improve the accuracy by over 70% for the prediction of the total energy and the splitting energy. With regard to semiempirical methods, our proposed models reduce the MAE, with respect to these methods, by 2 orders of magnitude.

Accurately determining the global minima of a molecular structure is important in diverse scientific fields, including drug design, materials science, and chemical synthesis. Conformational search engines serve as valuable tools for exploring the extensive conformational space of molecules and for identifying energetically favorable conformations. In this study, we present a comparison of Auto3D, CREST, Balloon, and ETKDG (from RDKit), which are freely available conformational search engines, to evaluate their effectiveness in locating global minima. These engines employ distinct methodologies, including machine learning (ML) potential-based, semiempirical, and force field-based approaches. To validate these methods, we propose the use of collisional cross-section (CCS) values obtained from ion mobility-mass spectrometry studies. We hypothesize that experimental gas-phase CCS values can provide experimental evidence that we likely have the global minimum for a given molecule. To facilitate this effort, we used our gas-phase conformation library (GPCL) which currently consists of the full ensembles of 20 small molecules and can be used by the community to validate any conformational search engine. Further members of the GPCL can be readily created for any molecule of interest using our standard workflow used to compute CCS values, expanding the ability of the GPCL in validation exercises. These innovative validation techniques enhance our understanding of the conformational landscape and provide valuable insights into the performance of conformational generation engines. Our findings shed light on the strengths and limitations of each search engine, enabling informed decisions for their utilization in various scientific fields, where accurate molecular structure determination is crucial for understanding biological activity and designing targeted interventions. By facilitating the identification of reliable conformations, this study significantly contributes to enhancing the efficiency and accuracy of molecular structure determination, with particular focus on metabolite structure elucidation. The findings of this research also provide valuable insights for developing effective workflows for predicting the structures of unknown compounds with high precision.

Metal ions play crucial roles in protein- and ligand-mediated interactions. They not only act as catalysts to facilitate biological processes but are also important as protein structural elements. Accurately predicting metal ion interactions in computational studies has always been a challenge, and various methods have been suggested to improve these interactions. One such method is the 12-6-4 Lennard-Jones (LJ)-type nonbonded model. Using this model, it has been possible to successfully reproduce the experimental properties of metal ions in aqueous solution. The model includes induced dipole interactions typically ignored in the standard 12-6 LJ nonbonded model. In this we expand the applicability of this model to metal ion-carboxylate interactions. Using 12-6-4 parameters that reproduce the solvation free energies of the metal ions leads to an overestimation of metal ion-acetate interactions, thus, prompting us to fine-tune the model to specifically handle the latter. We also show that the standard 12-6 LJ model significantly falls short in reproducing the experimental binding free energy between acetate and 11 metal ions (Ni(II), Mg(II), Cu(II), Zn(II), Co(II), Cu(I), Fe(II), Mn(II), Cd(II), Ca(II), and Ag(I)). In this study, we describe optimized C 4 parameters for the 12-6-4 LJ nonbonded model to be used with three widely employed water models (Transferable Intermolecular Potential with 3 Points (TIP3P), Simple Point Charge Extended (SPC/E), and Optimal Point Charge (OPC) water models). These parameters can accurately match the experimental binding free energy between 11 metal ions and acetate. These parameters can be applied to the study of metalloproteins and transition metal ion channels and transporters, as acetate serves as a representative of the negatively charged amino acid side chains from aspartate and glutamate.

Based on a series of energy minimizations with starting structures obtained from the Baker test set of 30 organic molecules, a comparison is made between various open-source geometry optimization codes that are interfaced with the open-source QUantum Interaction Computational Kernel (QUICK) program for gradient and energy calculations. The findings demonstrate how the choice of the coordinate system influences the optimization process to reach an equilibrium structure. With fewer steps, internal coordinates outperform Cartesian coordinates, while the choice of the initial Hessian and Hessian update method in quasi-Newton approaches made by different optimization algorithms also contributes to the rate of convergence. Furthermore, an available open-source machine learning method based on Gaussian process regression (GPR) was evaluated for energy minimizations over surrogate potential energy surfaces with both Cartesian and internal coordinates with internal coordinates outperforming Cartesian. Overall, geomeTRIC and DL-FIND with their default optimization method as well as with the GPR-based model using Hartree-Fock theory with the 6-31G** basis set needed a comparable number of geometry optimization steps to the approach of Baker using a unit matrix as the initial Hessian to reach the optimized geometry. On the other hand, the Berny and Sella offerings in ASE outperformed the other algorithms. Based on this, we recommend using the file-based approaches, ASE/Berny and ASE/Sella, for large-scale optimization efforts, while if using a single executable is preferable, we now distribute QUICK integrated with DL-FIND.

AmberTools is a free and open-source collection of programs used to set up, run, and analyze molecular simulations. The newer features contained within AmberTools23 are briefly described in this Application note.

We implemented an ab initio CCS prediction workflow which incrementally refines generated structures using molecular mechanics, a deep learning potential, conformational clustering, and quantum mechanics (QM). Automating intermediate steps for a high performance computing (HPC) environment allows users to input the SMILES structure of small organic molecules and obtain a Boltzmann averaged collisional cross section (CCS) value as output. The CCS of a molecular species is a metric measured by ion mobility spectrometry (IMS) which can improve annotation of untargeted metabolomics experiments. We report only a minor drop in accuracy when we expedite the CCS calculation by replacing the QM geometry refinement step with a single-point energy calculation. Even though the workflow involves stochastic steps (i.e., conformation generation and clustering), the final CCS value was highly reproducible for multiple iterations on L-carnosine. Finally, we illustrate that the gas phase ensembles modeled for the workflow are intermediate files which can be used for the prediction of other properties such as aqueous phase nuclear magnetic resonance chemical shift prediction. The software is available at the following link: https://github.com/DasSusanta/snakemake_ccs.

Atomic radii play important roles in scientific research. The covalent radii of atoms, ionic radii of ions, and van der Waals (VDW) radii of neutral atoms can all be derived from crystal structures. However, the VDW radii of ions are a challenge to determine because the atomic distances in crystal structures were determined by a combination of VDW interactions and electrostatic interactions, making it unclear how to define the VDW sphere of ions in such an environment. In the present study, we found that VDW radii, which were determined based on the 0.0015 au electron density contour through a wavefunction analysis on atoms, have excellent agreement with the VDW radii of noble-gas atoms determined experimentally. Based on this criterion, we calculated the VDW radii for various atomic ions across the periodic table, providing a systematic set of VDW radii of ions. Previously we have shown that the 12-6 Lennard-Jones nonbonded model could not simultaneously reproduce the hydration free energy (HFE) and ion-oxygen distance (IOD) for an atomic ion when its charge is +2 or higher. Because of this, we developed the 12-6-4 model to reproduce both properties at the same time by explicitly considering the ion-induced dipole interactions. However, recent studies showed it was possible to use the 12-6 model to simulate both properties simultaneously when an ion has the R min /2 parameter ( i.e ., the VDW radius) close to the Shannon ionic radius. In the present study, we show that such a "success" is due to an unphysical overfitting, as the VDW radius of an ion should be significantly larger than its ionic radius. Through molecular dynamics simulations, we show that such overfitting causes significant issues when transferring the parameters from ion-water systems to ion-ligand and metalloprotein systems. In comparison, the 12-6-4 model shows significant improvement in comparison to the overfitted 12-6 model, showing excellent transferability across different systems. In summary, although both the 12-6-4 and 12-6 models could reproduce HFE and IOD for an ion, the 12-6-4 model accomplishes such a task based on the consideration of the physics involved, while the 12-6 model accomplishes this through overfitting, which brings significant transferability issues when simulating other systems. Hence, we strongly recommend the use of the 12-6-4 model (or even more sophisticated models) instead of overfitted 12-6 models when simulating complex systems such as metalloproteins.

We have ported and optimized the graphics processing unit (GPU)-accelerated QUICK and AMBER-based ab initio quantum mechanics/molecular mechanics (QM/MM) implementation on AMD GPUs. This encompasses the entire Fock matrix build and force calculation in QUICK including one-electron integrals, two-electron repulsion integrals, exchange-correlation quadrature, and linear algebra operations. General performance improvements to the QUICK GPU code are also presented. Benchmarks carried out on NVIDIA V100 and AMD MI100 cards display similar performance on both hardware for standalone HF/DFT calculations with QUICK and QM/MM molecular dynamics simulations with QUICK/AMBER. Furthermore, with respect to the QUICK/AMBER release version 21, significant speedups are observed for QM/MM molecular dynamics simulations. This significantly increases the range of scientific problems that can be addressed with open-source QM/MM software on state-of-the-art computer hardware.

We develop a framework for the design of optimized alchemical transformation pathways in free energy simulations using nonlinear mixing and a new functional form for so-called "softcore" potentials. We describe the implementation and testing of this framework in the GPU-accelerated AMBER software suite. The new optimized alchemical transformation pathways integrate a number of important features, including (1) the use of smoothstep functions to stabilize behavior near the transformation end points, (2) consistent power scaling of Coulomb and Lennard-Jones (LJ) interactions with unitless control parameters to maintain balance of electrostatic attractions and exchange repulsions, (3) pairwise form based on the LJ contact radius for the effective interaction distance with separation-shifted scaling, and (4) rigorous smoothing of the potential at the nonbonded cutoff boundary. The new softcore potential form is combined with smoothly transforming nonlinear λ weights for mixing specific potential energy terms, along with flexible λ-scheduling features, to enable robust and stable alchemical transformation pathways. The resulting pathways are demonstrated and tested, and shown to be superior to the traditional methods in terms of numerical stability and minimal variance of the free energy estimates for all cases considered. The framework presented here can be used to design new alchemical enhanced sampling methods, and leveraged in robust free energy workflows for large ligand data sets.

Ion channels are specialized proteins located on the plasma membrane and control the movement of ions across the membrane. Zn ion plays an indispensable role as a structural constituent of various proteins, moreover, it plays an important dynamic role in cell signaling. In this chapter, we discuss computational insights into zinc efflux and influx mechanism through YiiP (from Escherichia coli and Shewanella oneidensis) and BbZIP (Bordetella bronchiseptica) transporters, respectively. Gaining knowledge about the mechanism of zinc transport at the molecular level can aid in developing treatments for conditions such as diabetes and cancer by manipulating extracellular and intracellular levels of zinc ions., (Copyright © 2023. Published by Elsevier Inc.)

Modeling the interaction between a metal ion and small molecules can provide pivotal information to bridge and close the gap between two types of simulations: metal ions in water and metal ions in metalloproteins. As previously established, the 12-6-4 Lennard-Jones (LJ)-type nonbonded model, because of its ability to account for the induced dipole effect, has been highly successful in simulating metal ion systems. Using the potential of mean force (PMF) method, the polarizability of the metal-chelating nitrogen from two types of imidazole molecules, delta nitrogen protonated (HID) and epsilon nitrogen protonated (HIE), has been parametrized against experiment for 11 metals (Ag(I), Ca(II), Cd(II), Co(II), Cu(I), Cu(II), Fe(II), Mg(II), Mn(II), Ni(II), and Zn(II)) in conjunction with three commonly used water models (TIP3P, SPC/E, and OPC). We show that the standard 12-6 and unmodified 12-6-4 models are not able to accurately model these interactions and, indeed, predict that the complex should be unstable. The resultant parameters further establish the flexibility and the reliability of the 12-6-4 LJ-type nonbonded model, which can correctly describe three-component interactions between a metal, ligand, and solvent by simply tuning the polarizability of the chelating atom. Also, the transferability of this model was tested, showing the capability of describing metal-ligand interactions in various environments.

The reactive force field (ReaxFF) model bridges the gap between traditional classical models and quantum mechanical (QM) models by incorporating dynamic bonding and polarizability. To achieve realistic simulations using ReaxFF, model parameters must be optimized against high fidelity training data which typically come from QM calculations. Existing parameter optimization methods for ReaxFF consist of black box techniques using genetic algorithms or Monte Carlo methods. Due to the stochastic behavior of these methods, the optimization process oftentimes requires millions of error evaluations for complex parameter fitting tasks, thereby significantly hampering the rapid development of high quality parameter sets. Rapid optimization of the parameters is essential for developing and refining Reax force fields because producing a force field which exhibits empirical accuracy in terms of dynamics typically requires multiple refinements to the training data as well as to the parameters under optimization. In this work, we present JAX-ReaxFF, a novel software tool that leverages modern machine learning infrastructure to enable fast optimization of ReaxFF parameters. By calculating gradients of the loss function using the JAX library, JAX-ReaxFF utilizes highly effective local optimization methods that are initiated from multiple guesses in the high dimensional optimization space to obtain high quality results. Leveraging the architectural portability of the JAX framework, JAX-ReaxFF can execute efficiently on multicore CPUs, graphics processing units (GPUs), or even tensor processing units (TPUs). As a result of using the gradient information and modern hardware accelerators, we are able to decrease ReaxFF parameter optimization time from days to mere minutes. Furthermore, the JAX-ReaxFF framework can also serve as a sandbox environment for domain scientists to explore customizing the ReaxFF functional form for more accurate modeling.

Quantum chemistry enables to study systems with chemical accuracy (<1 kcal/mol from experiment) but is restricted to a handful of atoms due to its computational expense. This has led to ongoing interest to optimize and simplify these methods while retaining accuracy. Implementing quantum mechanical (QM) methods on modern hardware such as multiple-GPUs is one example of how the field is optimizing performance. Multiscale approaches like the so-called QM/molecular mechanical method are gaining popularity in drug discovery because they focus the application of QM methods on the region of choice (e.g., the binding site), while using efficient MM models to represent less relevant areas. The creation of simplified QM methods is another example, including the use of machine learning to create ultra-fast and accurate QM models. Herein, we summarize recent advancements in the development of optimized QM methods that enhance our ability to use these methods in computer aided drug discovery., (Copyright © 2022 Elsevier Ltd. All rights reserved.)

The COVID-19 outbreak has been devastating, with hundreds of millions of infections and millions of deaths reported worldwide. In response, the application of structure-activity relationships (SAR) upon experimentally validated inhibitors of SARS-CoV-2 main protease (M pro ) may provide an avenue for the identification of new lead compounds active against COVID-19. Upon the basis of information gleaned from a combination of reported crystal structures and the docking of experimentally validated inhibitors, four "rules" for designing potent M pro inhibitors have been proposed. The aim here is to guide medicinal chemists toward the most probable hits and to provide guidance on repurposing available structures as M pro inhibitors. Experimental examination of our own previously reported inhibitors using the four "rules" identified a potential lead compound, the cathepsin inhibitor GB111-NH 2 , that was 2.3 times more potent than SARS-CoV-2 M pro inhibitor N3 .

A novel locally polarizable multisite model based on the original cation dummy atom (CDA) model is described for molecular dynamics simulations of ions in condensed phases. Polarization effects are introduced by the electronegativity equalization model (EEM) method where charges on the metal ion and its dummy atoms can fluctuate to respond to the environment. This model includes explicit polarization and ion-induced interactions and can be coupled with nonpolarizable or polarizable water models, making it more transferable to simpler force fields. This approach allows us to enhance the original fixed charge CDA model where the charge distribution cannot adapt to the local solvent structure. To illustrate the new CDA pol model, we examined properties of the Zn 2+ , Al 3+ , and Zr 4+ ions in aqueous solution. The polarizable model and Lennard-Jones parameters were refined for octahedrally coordinated Zn 2+ , Al 3+ , and Zr 4+ CDAs to reproduce thermodynamic and geometrical properties. Using this locally polarizable model, we were able to obtain the experimental hydration free energy, ion-oxygen distance, and coordination number coupled with the standard 12-6 Lennard-Jones model. This model can be used in myriad additional applications where local polarization and charge transfer is important.

As a fundamental property of all fluids, diffusion plays myriad roles in both science and our daily lives. Diffusive properties of many liquids including water have been extensively studied both experimentally and theoretically, while for transition metal ions, there exist significant experimental data that have not been extensively studied theoretically. Hence, high-confidence predictions for challenging systems like radioactive ions that are biohazardous cannot be reliably generated. In this work, a workflow named ISAIAH (Ion Simulation using AMBER for dIffusion Action when Hydrated) was designed to accurately simulate the diffusion coefficients of 15 monoatomic ions with charges varying from -1 to +3 in four water models. As the results indicate, good agreement with experimental values was achieved, leading us to select 239 Pu 4+ (for which no experimental data are available) as a candidate ion to make a theoretical prediction of its diffusion coefficient in water. Among all the force field parameter sets, the ones parametrized using an augmented 12-6-4 Lennard-Jones (LJ) potential showed lower average unsigned errors (AUE) for ions of various radii and electron configurations relative to some 12-6 LJ parameters. This observation agrees well with the fact that diffusion is affected by both the hydration free energy (HFE) and the ion-oxygen distance (IOD) between solute and solvent molecules, both of which are handled well by the 12-6-4 model.

The interpretation of ion mobility coupled to mass spectrometry (IM-MS) data to predict unknown structures is challenging and depends on accurate theoretical estimates of the molecular ion collision cross section (CCS) against a buffer gas in a low or atmospheric pressure drift chamber. The sensitivity and reliability of computational prediction of CCS values depend on accurately modeling the molecular state over accessible conformations. In this work, we developed an efficient CCS computational workflow using a machine learning model in conjunction with standard DFT methods and CCS calculations. Furthermore, we have performed Traveling Wave IM-MS (TWIMS) experiments to validate the extant experimental values and assess uncertainties in experimentally measured CCS values. The developed workflow yielded accurate structural predictions and provides unique insights into the likely preferred conformation analyzed using IM-MS experiments. The complete workflow makes the computation of CCS values tractable for a large number of conformationally flexible metabolites with complex molecular structures.

Scoring functions are the essential component in molecular docking methods. An accurate scoring function is expected to distinguish the native ligand pose from decoy poses. Our previous experience (Pei et al. J. Chem. Inf. Model. 2019, 59 (7), 3305-3315) proved that combining the random forest (RF) algorithm with knowledge-based potential functions can emphasize germane pair wise interactions and improve the performance of original knowledge-based potential functions on protein-ligand decoy detection. One of the most important potential function classes is the force field (FF) potential with one example being the Amber collection of FFs, which are widely available in the AMBER suite of simulation programs. However, for use in RF modeling studies, one needs pair wise energies that are hard to directly extract from Amber. To address this issue, FFENCODER-PL was constructed to calculate the pair wise energies based on the FF14SB and GAFF2 FFs in Amber. FFENCODER-PL was validated using 275 ligand and 21 protein-ligand structures. RF models were built by combining an RF classification algorithm with the pair wise energies calculated from FFENCODER-PL. CASF-2016 (Su et al. J. Chem. Inf. Model. 2019, 59, 895-913) was employed to test the performance of the resultant RF models, which outperformed 33 scoring functions on accuracy and native ranking tests. For the best decoy RMSD test, RF models give a best decoy with an RMSD of around 2 A ̊ from the native pose after including the best decoy-decoy comparisons in the RF model. The relative importance of the RF algorithm and force field potentials was also tested with the results suggesting that both the RF algorithm and force field potentials are important and combining them is the only way to achieve high accuracy. Finally, FFENCODER-PL makes force field-based pair wise energies available for further development of machine learning-based scoring functions. The codes and data used in this paper can be found at https://github.com/JunPei000/Amber&#95;protein&#95;ligand&#95;encoding.

We report a new multi-GPU capable ab initio Hartree-Fock/density functional theory implementation integrated into the open source QUantum Interaction Computational Kernel (QUICK) program. Details on the load balancing algorithms for electron repulsion integrals and exchange correlation quadrature across multiple GPUs are described. Benchmarking studies carried out on up to four GPU nodes, each containing four NVIDIA V100-SXM2 type GPUs demonstrate that our implementation is capable of achieving excellent load balancing and high parallel efficiency. For representative medium to large size protein/organic molecular systems, the observed parallel efficiencies remained above 82% for the Kohn-Sham matrix formation and above 90% for nuclear gradient calculations. The accelerations on NVIDIA A100, P100, and K80 platforms also have realized parallel efficiencies higher than 68% in all tested cases, paving the way for large-scale ab initio electronic structure calculations with QUICK.

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.

The quantum mechanics/molecular mechanics (QM/MM) approach is an essential and well-established tool in computational chemistry that has been widely applied in a myriad of biomolecular problems in the literature. In this publication, we report the integration of the QUantum Interaction Computational Kernel (QUICK) program as an engine to perform electronic structure calculations in QM/MM simulations with AMBER. This integration is available through either a file-based interface (FBI) or an application programming interface (API). Since QUICK is an open-source GPU-accelerated code with multi-GPU parallelization, users can take advantage of "free of charge" GPU-acceleration in their QM/MM simulations. In this work, we discuss implementation details and give usage examples. We also investigate energy conservation in typical QM/MM simulations performed at the microcanonical ensemble. Finally, benchmark results for two representative systems in bulk water, the N -methylacetamide (NMA) molecule and the photoactive yellow protein (PYP), show the performance of QM/MM simulations with QUICK and AMBER using a varying number of CPU cores and GPUs. Our results highlight the acceleration obtained from a single or multiple GPUs; we observed speedups of up to 53× between a single GPU vs a single CPU core and of up to 2.6× when comparing four GPUs to a single GPU. Results also reveal speedups of up to 3.5× when the API is used instead of FBI.

While accurately modeling the conformational ensemble is required for predicting properties of flexible molecules, the optimal method of obtaining the conformational ensemble appears as varied as their applications. Ensemble structures have been modeled by generation, refinement, and clustering of conformations with a sufficient number of samples. We present a conformational clustering algorithm intended to automate the conformational clustering step through the Louvain algorithm, which requires minimal hyperparameters and importantly no predefined number of clusters or threshold values. The conformational graphs produced by this method for O -succinyl-l-homoserine, oxidized nicotinamide adenine dinucleotide, and 200 representative metabolites each preserved the geometric/energetic correlation expected for points on the potential energy surface. Clustering based on these graphs provides partitions informed by the potential energy surface. Automating conformational clustering in a workflow with AutoGraph may mitigate human biases introduced by guess and check over hyperparameter selection while allowing flexibility to the result by not imposing predefined criteria other than optimizing the model's loss function. Associated codes are available at https://github.com/TanemuraKiyoto/AutoGraph.

Commonly seen in rare-earth chemistry and materials science, highly charged metal ions play key roles in many chemical processes. Computer simulations have become an important tool for scientific research nowadays. Meaningful simulations require reliable parameters. In the present work, we parametrized 18 M(III) and 6 M(IV) metal ions for four new water models (OPC3, OPC, TIP3P-FB, TIP4P-FB) in conjunction with each of the 12-6 and 12-6-4 nonbonded models. Similar to what was observed previously, issues with the 12-6 model can be fixed by using the 12-6-4 model. Moreover, the four new water models showed comparable performance or considerable improvement over the previous water models (TIP3P, SPC/E, and TIP4P EW ) in the same category (3-point or 4-point water models, respectively). Finally, we reported a study of a metalloprotein system demonstrating the capability of the 12-6-4 model to model metalloproteins. The reported parameters will facilitate accurate simulations of highly charged metal ions in aqueous solution.

Monovalent ions play significant roles in various biological and material systems. Recently, four new water models (OPC3, OPC, TIP3P-FB, and TIP4P-FB), with significantly improved descriptions of condensed phase water, have been developed. The pairwise interaction between the metal ion and water necessitates the development of ion parameters specifically for these water models. Herein, we parameterized the 12-6 and the 12-6-4 nonbonded models for 12 monovalent ions with the respective four new water models. These monovalent ions contain eight cations including alkali metal ions (Li + , Na + , K + , Rb + , Cs + ), transition-metal ions (Cu + and Ag + ), and Tl + from the boron family, along with four halide anions (F - , Cl - , Br - , I - ). Our parameters were designed to reproduce the target hydration free energies (the 12-6 hydration free energy (HFE) set), the ion-oxygen distances (the 12-6 ion-oxygen distance (IOD) set), or both of them (the 12-6-4 set). The 12-6-4 parameter set provides highly accurate structural features overcoming the limitations of the routinely used 12-6 nonbonded model for ions. Specifically, we note that the 12-6-4 parameter set is able to reproduce experimental hydration free energies within 1 kcal/mol and experimental ion-oxygen distances within 0.01 Å simultaneously. We further reproduced the experimentally determined activity derivatives for salt solutions, validating the ion parameters for simulations of ion pairs. The improved performance of the present water models over our previous parameter sets for the TIP3P, TIP4P, and SPC/E water models (Li, P. et al J. Chem. Theory Comput. 2015 11 1645 1657) highlights the importance of the choice of water model in conjunction with the metal ion parameter set.

The MCPB.py program greatly facilitates force field parameterization for metal sites in metalloproteins and organometallic compounds. Herein we present an example of MCPB.py to the parameterization of the dioxygen binding metal site of peptidylglycine-alphahydroxylating monooxygenase (PHM), which contains a copper ion. In this example, we also extend the functionality of MCPB.py to support molecular dynamics (MD) simulations in GROMACS through a python script. Illustrative MD simulations were performed using GROMACS and the results were analyzed. Notes about the program were also provided in this chapter, to assist MCPB.py users for metal site parameterizations.

Modern scientometric techniques, applied at scale, can provide valuable information that complements qualitative investigation of the accumulation of knowledge in a field. We discuss a trio of articles from computational chemistry selected from an analysis of 181 million tri-cited articles.

Combined quantum mechanical/molecular mechanical (QM/MM) models using semiempirical and ab initio methods have been extensively reported on over the past few decades. These methods have been shown to be capable of providing unique insights into a range of problems, but they are still limited to relatively short time scales, especially QM/MM models using ab initio methods. An intermediate approach between a QM based model and classical mechanics could help fill this time-scale gap and facilitate the study of a range of interesting problems. Reactive force fields represent the intermediate approach explored in this paper. A widely used reactive model is ReaxFF, which has largely been applied to materials science problems and is generally used as a stand-alone (i.e., the full system is modeled using ReaxFF). We report a hybrid ReaxFF/AMBER molecular dynamics (MD) tool, which introduces ReaxFF capabilities to capture bond breaking and formation within the AMBER MD software package. This tool enables us to study local reactive events in large systems at a fraction of the computational costs of QM/MM models. We describe the implementation of ReaxFF/AMBER, validate this implementation using a benzene molecule solvated in water, and compare its performance against a range of similar approaches. To illustrate the predictive capabilities of ReaxFF/AMBER, we carried out a Claisen rearrangement study in aqueous solution. In a first for ReaxFF, we were able to use AMBER's potential of mean force (PMF) capabilities to perform a PMF study on this organic reaction. The ability to capture local reaction events in large systems using combined ReaxFF/AMBER opens up a range of problems that can be tackled using this model to address both chemical and biological processes.

Protein-protein interactions (PPIs) are ubiquitous and functionally of great importance in biological systems. Hence, the accurate prediction of PPIs by protein-protein docking and scoring tools is highly desirable in order to characterize their structure and biological function. Ab initio docking protocols are divided into the sampling of docking poses to produce at least one near-native structure, and then to evaluate the vast candidate structures by scoring. Concurrent development in both sampling and scoring is crucial for the deployment of protein-protein docking software. In the present work, we apply a machine learning model on pairwise potentials to refine the task of protein quaternary structure native structure detection among decoys. A decoy set was featurized using the Knowledge and Empirical Combined Scoring Algorithm 2 (KECSA2) pairwise potential. The highly unbalanced decoy set was then balanced using a comparison concept between native and decoy structures. The resultant comparison descriptors were used to train a logistic regression (LR) classifier. The LR model yielded the optimal performance for native detection among decoys compared with conventional scoring functions, while exhibiting lesser performance for the detection of low root mean square deviation decoy structures. Its deployment on an independent benchmark set confirms that the scoring function performs competitively relative to other scoring functions. The scripts used are available at https://github.com/TanemuraKiyoto/PPI-native-detection-via-LR., (© 2020 Wiley Periodicals LLC.)

The goal of the present manuscript is to succinctly trace the key technological steps in the evolution of alchemical free energy methods (AFEMs) from a purely theoretical construct to a method that is now widely used in the biotechnological and pharmaceutical industries. More specifically, we focus on relative binding free energy (RBFE) computations which are more routinely applied in computer aided drug design (CADD) campaigns rather than the more computationally intensive absolute binding free energy (ABFE) computations. We have not been exhaustive in the development of our timeline but rather try to weave a story about how theoretical ideas were ultimately converted into contemporary free energy capabilities. Necessarily this story-telling approach limits us from citing all work on AFEMs, and we apologize for this shortcoming. However, for those interested in a broad delineation of all the work done in this area they are directed to the many excellent reviews that are extant.

Herein we provide high-precision validation tests of the latest GPU-accelerated free energy code in AMBER. We demonstrate that consistent free energy results are obtained in both the gas phase and in solution. We first show, in the context of thermodynamic integration (TI), that the results are invariant with respect to "split" (e.g., stepwise decharge-vdW-recharge) versus "unified" protocols. This brought to light a subtle inconsistency in previous versions of AMBER that was traced to the improper treatment of 1-4 vdW and electrostatic interactions involving atoms across the softcore boundary. We illustrate that under the assumption that the ensembles produced by different legs of the alchemical transformation between molecules A and B in the gas phase and aqueous phase are very small, the inconsistency in the relative hydration free energy ΔΔ G hydr [A → B] = Δ G aq [A → B] - Δ G gas [A → B] is minimal. However, for general cases where the ensembles are shown to be substantially different, as expected in ligand-protein binding applications, these errors can be large. Finally, we demonstrate that results for relative hydration free energy simulations are independent of TI or multistate Bennett's acceptance ratio (MBAR) analysis, invariant to the specific choice of the softcore region, and in agreement with results derived from absolute hydration free energy values.

MRP.py is a Python-based parametrization program for covalently modified amino acid residues for molecular dynamics simulations. Charge derivation is performed via an RESP charge fit, and force constants are obtained through rewriting of either protein or GAFF database parameters. This allows for the description of interfacial interactions between the modifed residue and protein. MRP.py is capable of working with a variety of protein databases. MRP.py's highly general and systematic method of obtaining parameters allows the user to circumvent the process of parametrizing the modified residue-protein interface. Two examples, a covalently bound inhibitor and covalent adduct consisting of modified residues, are provided in the Supporting Information.

To obtain accurate and converged free energy calculations for ligand binding to biomolecular systems requires validated force fields and extensive sampling of the energy landscape, which requires exhaustive and effective conformational searching methods. Herein, we introduce the consecutive histograms Monte Carlo (CHMC) sampling protocol that generates receptor-ligand binding modes within a series of continuously distributed sampling units ranging from placement near the geometric center of the receptor's binding site to fully unbound states. This protocol employs independent energy-state sampling for calculating the ensemble energy within every predefined location along the receptor-ligand dissociation pathway, without the need to traverse the energy barriers as in molecular dynamic simulations during the dissociation procedure. We applied this method to a set of selected receptor targets with their corresponding ligands providing detailed studies of molecular binding free energy predictions. The results show that the CHMC gives an excellent accounting of the free energy surfaces and binding free energies at a reasonable computational cost.

Atom pairwise potential functions make up an essential part of many scoring functions for protein decoy detection. With the development of machine learning (ML) tools, there are multiple ways to combine potential functions to create novel ML models and methods. Potential function parameters can be easily extracted; however, it is usually hard to directly obtain the calculated atom pairwise energies from scoring functions. Amber, as one of the most popular suites of modeling programs, has an extensive history and library of force field potential functions. In this work, we directly used the force field parameters in ff94 and ff14SB from Amber and encoded them to calculate atom pairwise energies for different interactions. Two sets of structures (single amino acid set and a dipeptide set) were used to evaluate the performance of our encoded Amber potentials. From the comparison results between energy terms obtained from our encoding and Amber, we find energy difference within ±0.06 kcal/mol for all tested structures. Previously we have shown that the Random Forest (RF) model can help to emphasize more important atom pairwise interactions and ignore insignificant ones [Pei, J.; Zheng, Z.; Merz, K. M. J. Chem. Inf. Model. 2019, 59, 1919-1929]. Here, as an example of combining ML methods with traditional potential functions, we followed the same work flow to combine the RF models with force field potential functions from Amber. To determine the performance of our RF models with force field potential functions, 224 different protein native-decoy systems were used as our training and testing sets We find that the RF models with ff94 and ff14SB force field parameters outperformed all other scoring functions (RF models with KECSA2, RWplus, DFIRE, dDFIRE, and GOAP) considered in this work for native structure detection, and they performed similarly in detecting the best decoy. Through inclusion of best decoy to decoy comparisons in building our RF models, we were able to generate models that outperformed the score functions tested herein both on accuracy and best decoy detection, again showing the performance and flexibility of our RF models to tackle this problem. Finally, the importance of the RF algorithm and force field parameters were also tested and the comparison results suggest that both the RF algorithm and force field potentials are important with the ML scoring function achieving its best performance only by combining them together. All code and data used in this work are available at https://github.com/JunPei000/FFENCODER&#95;for&#95;Protein&#95;Folding&#95;Pose&#95;Selection.

A major challenge for metabolomic analysis is to obtain an unambiguous identification of the metabolites detected in a sample. Among metabolomics techniques, NMR spectroscopy is a sophisticated, powerful, and generally applicable spectroscopic tool that can be used to ascertain the correct structure of newly isolated biogenic molecules. However, accurate structure prediction using computational NMR techniques depends on how much of the relevant conformational space of a particular compound is considered. It is intrinsically challenging to calculate NMR chemical shifts using high-level DFT when the conformational space of a metabolite is extensive. In this work, we developed NMR chemical shift calculation protocols using a machine learning model in conjunction with standard DFT methods. The pipeline encompasses the following steps: (1) conformation generation using a force field (FF)-based method, (2) filtering the FF generated conformations using the ASE-ANI machine learning model, (3) clustering of the optimized conformations based on structural similarity to identify chemically unique conformations, (4) DFT structural optimization of the unique conformations, and (5) DFT NMR chemical shift calculation. This protocol can calculate the NMR chemical shifts of a set of molecules using any available combination of DFT theory, solvent model, and NMR-active nuclei, using both user-selected reference compounds and/or linear regression methods. Our protocol reduces the overall computational time by 2 orders of magnitude over methods that optimize the conformations using fully ab initio methods, while still producing good agreement with experimental observations. The complete protocol is designed in such a manner that makes the computation of chemical shifts tractable for a large number of conformationally flexible metabolites.