Amber free energy tools: Interoperable software for free energy simulations using generalized quantum mechanical/molecular mechanical and machine learning potentials.

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.

Binding of Phosphate Species to Ca2+ and Mg2+ in Aqueous Solution.

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 H2PO4-, HPO42-, PO43- coordinating with Ca2+ and Mg2+. 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 Ca2+ and Mg2+ 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.

LigandDiff: de Novo Ligand Design for 3D Transition Metal Complexes with Diffusion Models.

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.

Modeling Zinc Complexes Using Neural Networks.

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.

Modeling Fe(II) Complexes Using Neural Networks.

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.

Thermodynamics of Metal-Acetate Interactions.

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 C4 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.

Molecular Gas-Phase Conformational Ensembles.

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.

Geometry Optimization: A Comparison of Different Open-Source Geometry Optimizers.

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.

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.

Getting zinc into and out of cells.

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.

Rapid and Automated Ab Initio Metabolite Collisional Cross Section Prediction from SMILES Input.

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.

Consequences of Overfitting the van der Waals Radii of Ions.

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 Rmin/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.

AMBER Free Energy Tools: A New Framework for the Design of Optimized Alchemical Transformation Pathways.

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.