Anisotropic, Polarizable Molecular Mechanics Studies of Inter- and Intramolecular Interactions and Ligand-Macromolecule Complexes. A Bottom-Up Strategy

Advancing Force Fields Parameterization: A Directed Graph Attention Networks Approach

Force fields (FFs) are an established tool for simulating large and complex molecular systems. However, parametrizing FFs is a challenging and time-consuming task that relies on empirical heuristics, experimental data, and computational data. Recent efforts aim to automate the assignment of FF parameters using pre-existing databases and on-the-fly ab initio data. In this study, we propose a graph-based force field (GB-FFs) model to directly derive parameters for the Generalized Amber Force Field (GAFF) from chemical environments and research into the influence of functional forms. Our end-to-end parametrization approach predicts parameters by aggregating the basic information in directed molecular graphs, eliminating the need for expert-defined procedures and enhances the accuracy and transferability of GAFF across a broader range of molecular complexes. Simulation results are compared to the original GAFF parametrization. In practice, our results demonstrate an improved transferability of the model, showcasing its improved accuracy in modeling intermolecular and torsional interactions, as well as improved solvation free energies. The optimization approach developed in this work is fully applicable to other nonpolarizable FFs as well as to polarizable ones.

Beyond isotropic repulsion: Classical anisotropic repulsion by inclusion of p orbitals

Accurate modeling of intermolecular repulsion is an integral component in force field development. Although repulsion can be explicitly calculated by applying the Pauli exclusion principle, this approach is computationally viable only for systems of limited sizes. Instead, it has previously been shown that repulsion can be reformulated in a "classical" picture: the Pauli exclusion principle prohibits electrons from occupying the same state, leading to a depletion of electronic charge between atoms, giving rise to an enhanced nuclear-nuclear electrostatic repulsion. This classical picture is called the isotropic S2/R approximation, where S is the overlap and R is the interatomic distance. This approximation accurately captures the repulsion of isotropic atoms such as noble gas dimers; however, a key deficiency is that it fails to capture the angular dependence of the repulsion of anisotropic molecules. To include directionality, the wave function must at least be a linear combination of s and p orbitals. We derive a new anisotropic S2/R repulsion model through the inclusion of the anisotropic p orbital term in the total wave function. Because repulsion is pairwise and decays rapidly, it can be truncated at a short range, making it amenable for efficient calculation of energy and forces in complex biomolecular systems. We present a parameterization of the S101 dimer database against the ab initio benchmark symmetry-adapted perturbation theory, which yields an rms error of only 0.9 kcal/mol. The importance of the anisotropic term is demonstrated through angular scans of water-water dimers and dimers involving halobenzene. Simulation of liquid water shows that the model can be computed efficiently for realistic system sizes.

Assessment of Amino Acid Electrostatic Parametrizations of the Polarizable Gaussian Multipole Model

Accurate parametrization of amino acids is pivotal for the development of reliable force fields for molecular modeling of biomolecules such as proteins. This study aims to assess amino acid electrostatic parametrizations with the polarizable Gaussian Multipole (pGM) model by evaluating the performance of the pGM-perm (with atomic permanent dipoles) and pGM-ind (without atomic permanent dipoles) variants compared to the traditional RESP model. The 100-conf-combterm fitting strategy on tetrapeptides was adopted, in which (1) all peptide bond atoms (-CO-NH-) share identical set of parameters and (2) the total charges of the two terminal N-acetyl (ACE) and N-methylamide (NME) groups were set to neutral. The accuracy and transferability of electrostatic parameters across peptides with varying lengths and real-world examples were examined. The results demonstrate the enhanced performance of the pGM-perm model in accurately representing the electrostatic properties of amino acids. This insight underscores the potential of the pGM-perm model and the 100-conf-combterm strategy for the future development of the pGM force field.

Multiscale computational simulation of pollutant behavior at water interfaces

The investigation of pollutant behavior at water interfaces is critical to understand pollution in aquatic systems. Computational methods allow us to overcome the limitations of experimental analysis, delivering valuable insights into the chemical mechanisms and structural characteristics of pollutant behavior at interfaces across a range of scales, from microscopic to mesoscopic. Quantum mechanics, all-atom molecular dynamics simulations, coarse-grained molecular dynamics simulations, and dissipative particle dynamics simulations represent diverse molecular interaction calculation methods that can effectively model pollutant behavior at environmental interfaces from atomic to mesoscopic scales. These methods provide a rich variety of information on pollutant interactions with water surfaces. This review synthesizes the advancements in applying typical computational methods to the formation, adsorption, binding, and catalytic conversion of pollutants at water interfaces. By drawing on recent advancements, we critically examine the current challenges and offer our perspective on future directions. This review seeks to advance our understanding of computational techniques for elucidating pollutant behavior at water interfaces, a critical aspect of water research.

Enforcing Local DNA Kinks by Sequence-Selective Trisintercalating Oligopeptides of a Tricationic Porphyrin: A Polarizable Molecular Dynamics Study

Bisacridinyl-bisarginyl porphyrin (BABAP) is a trisintercalating derivative of a tricationic porphyrin, formerly designed and synthesized in order to selectively target and photosensitize the ten-base pair palindromic sequence d(CGGGCGCCCG)2 . We resorted to the previously derived (Far et al., 2004) lowest energy-minimized (EM) structure of the BABAP complex with this sequence as a starting point. We performed polarizable molecular dynamics (MD) on this complex. It showed, over a 150 ns duration, the persistent binding of the Arg side-chain on each BABAP arm to the two G bases upstream from the central porphyrin intercalation site. We subsequently performed progressive shortenings of the connector chain linking the Arg-Gly backbone to the acridine, from n=6 methylenes to 4, followed by removal of the Gly backbone and further connector shortenings, from n=4 to n=1. These resulted into progressive deformations ('kinks') of the DNA backbone. In its most accented kinked structure, the DNA backbone was found to have a close overlap with that of DNA bound to Cre recombinase, with, at the level of one acridine intercalation site, negative roll and positive tilt values consistent with those experimentally found for this DNA at its own kinked dinucleotide sequence. Thus, in addition to their photosensitizing properties, some BABAP derivatives could induce sequence-selective, controlled DNA deformations, which are targets for cleavage by endonucleases or for repair enzymes.

Force-field-enhanced neural network interactions: from local equivariant embedding to atom-in-molecule properties and long-range effects

We introduce FENNIX (Force-Field-Enhanced Neural Network InteraXions), a hybrid approach between machine-learning and force-fields. We leverage state-of-the-art equivariant neural networks to predict local energy contributions and multiple atom-in-molecule properties that are then used as geometry-dependent parameters for physically-motivated energy terms which account for long-range electrostatics and dispersion. Using high-accuracy ab initio data (small organic molecules/dimers), we trained a first version of the model. Exhibiting accurate gas-phase energy predictions, FENNIX is transferable to the condensed phase. It is able to produce stable Molecular Dynamics simulations, including nuclear quantum effects, for water predicting accurate liquid properties. The extrapolating power of the hybrid physically-driven machine learning FENNIX approach is exemplified by computing: (i) the solvated alanine dipeptide free energy landscape; (ii) the reactive dissociation of small molecules.

Combining Force Fields and Neural Networks for an Accurate Representation of Chemically Diverse Molecular Interactions

A key goal of molecular modeling is the accurate reproduction of the true quantum mechanical potential energy of arbitrary molecular ensembles with a tractable classical approximation. The challenges are that analytical expressions found in general purpose force fields struggle to faithfully represent the intermolecular quantum potential energy surface at close distances and in strong interaction regimes; that the more accurate neural network approximations do not capture crucial physics concepts, e.g., nonadditive inductive contributions and application of electric fields; and that the ultra-accurate narrowly targeted models have difficulty generalizing to the entire chemical space. We therefore designed a hybrid wide-coverage intermolecular interaction model consisting of an analytically polarizable force field combined with a short-range neural network correction for the total intermolecular interaction energy. Here, we describe the methodology and apply the model to accurately determine the properties of water, the free energy of solvation of neutral and charged molecules, and the binding free energy of ligands to proteins. The correction is subtyped for distinct chemical species to match the underlying force field, to segment and reduce the amount of quantum training data, and to increase accuracy and computational speed. For the systems considered, the hybrid ab initio parametrized Hamiltonian reproduces the two-body dimer quantum mechanics (QM) energies to within 0.03 kcal/mol and the nonadditive many-molecule contributions to within 2%. Simulations of molecular systems using this interaction model run at speeds of several nanoseconds per day.

Perspectives on Ligand/Protein Binding Kinetics Simulations: Force Fields, Machine Learning, Sampling, and User-Friendliness

Computational techniques applied to drug discovery have gained considerable popularity for their ability to filter potentially active drugs from inactive ones, reducing the time scale and costs of preclinical investigations. The main focus of these studies has historically been the search for compounds endowed with high affinity for a specific molecular target to ensure the formation of stable and long-lasting complexes. Recent evidence has also correlated the in vivo drug efficacy with its binding kinetics, thus opening new fascinating scenarios for ligand/protein binding kinetic simulations in drug discovery. The present article examines the state of the art in the field, providing a brief summary of the most popular and advanced ligand/protein binding kinetics techniques and evaluating their current limitations and the potential solutions to reach more accurate kinetic models. Particular emphasis is put on the need for a paradigm change in the present methodologies toward ligand and protein parametrization, the force field problem, characterization of the transition states, the sampling issue, and algorithms' performance, user-friendliness, and data openness.

Fitting Potential Energy Surfaces by Learning the Charge Density Matrix with Permutationally Invariant Polynomials

The electronic energy in the Hartree-Fock (HF) theory is the trace of the product of the charge density matrix (CDM) with the one-electron and two-electron matrices represented in an atomic orbital basis, where the two-electron matrix is also a function of the same CDM. In this work, we examine a formalism of analytic representation of a generic molecular potential energy surface (PES) as a sum of a linearly parameterized HF and a correction term, the latter formally representing the electron correlation energy, also linearly parameterized, by expressing the elements of CDM using permutationally invariant polynomials (PIPs). We show on a variety of numerical examples, ranging from exemplary two-electron systems HeH+ and H3+ to the more challenging cases of methanium (CH5+) fragmentation and high-energy tautomerization of formamide to formimidic acid that such a formulation requires significantly fewer, 10-20% of PIPs, to accomplish the same accuracy of the fit as the conventional representation at practically the same computational cost.

Scalable hybrid deep neural networks/polarizable potentials biomolecular simulations including long-range effects

Deep-HP is a scalable extension of the Tinker-HP multi-GPU molecular dynamics (MD) package enabling the use of Pytorch/TensorFlow Deep Neural Network (DNN) models. Deep-HP increases DNNs' MD capabilities by orders of magnitude offering access to ns simulations for 100k-atom biosystems while offering the possibility of coupling DNNs to any classical (FFs) and many-body polarizable (PFFs) force fields. It allows therefore the introduction of the ANI-2X/AMOEBA hybrid polarizable potential designed for ligand binding studies where solvent-solvent and solvent-solute interactions are computed with the AMOEBA PFF while solute-solute ones are computed by the ANI-2X DNN. ANI-2X/AMOEBA explicitly includes AMOEBA's physical long-range interactions via an efficient Particle Mesh Ewald implementation while preserving ANI-2X's solute short-range quantum mechanical accuracy. The DNN/PFF partition can be user-defined allowing for hybrid simulations to include key ingredients of biosimulation such as polarizable solvents, polarizable counter ions, etc.… ANI-2X/AMOEBA is accelerated using a multiple-timestep strategy focusing on the model's contributions to low-frequency modes of nuclear forces. It primarily evaluates AMOEBA forces while including ANI-2X ones only via correction-steps resulting in an order of magnitude acceleration over standard Velocity Verlet integration. Simulating more than 10 μs, we compute charged/uncharged ligand solvation free energies in 4 solvents, and absolute binding free energies of host-guest complexes from SAMPL challenges. ANI-2X/AMOEBA average errors are discussed in terms of statistical uncertainty and appear in the range of chemical accuracy compared to experiment. The availability of the Deep-HP computational platform opens the path towards large-scale hybrid DNN simulations, at force-field cost, in biophysics and drug discovery.

Routine Molecular Dynamics Simulations Including Nuclear Quantum Effects: From Force Fields to Machine Learning Potentials

We report the implementation of a multi-CPU and multi-GPU massively parallel platform dedicated to the explicit inclusion of nuclear quantum effects (NQEs) in the Tinker-HP molecular dynamics (MD) package. The platform, denoted Quantum-HP, exploits two simulation strategies: the Ring-Polymer Molecular Dynamics (RPMD) that provides exact structural properties at the cost of a MD simulation in an extended space of multiple replicas and the adaptive Quantum Thermal Bath (adQTB) that imposes the quantum distribution of energy on a classical system via a generalized Langevin thermostat and provides computationally affordable and accurate (though approximate) NQEs. We discuss some implementation details, efficient numerical schemes, and parallelization strategies and quickly review the GPU acceleration of our code. Our implementation allows an efficient inclusion of NQEs in MD simulations for very large systems, as demonstrated by scaling tests on water boxes with more than 200,000 atoms (simulated using the AMOEBA polarizable force field). We test the compatibility of the approach with Tinker-HP's recently introduced Deep-HP machine learning potentials module by computing water properties using the DeePMD potential with adQTB thermostatting. Finally, we show that the platform is also compatible with the alchemical free energy estimation capabilities of Tinker-HP and fast enough to perform simulations. Therefore, we study how NQEs affect the hydration free energy of small molecules solvated with the recently developed Q-AMOEBA water force field. Overall, the Quantum-HP platform allows users to perform routine quantum MD simulations of large condensed-phase systems and will help to shed new light on the quantum nature of important interactions in biological matter.

Accurate and Interpretable Dipole Interaction Model-Based Machine Learning for Molecular Polarizability

Polarizabilities play significant roles in describing dispersive and inductive interactions of the atom and molecular systems. However, an accurate prediction of molecular polarizabilities from first principles is computationally prohibitive. Although physical models or statistical machine learning models have been proposed, either a lack of accurate description of local chemical environments or demanding a large number of samples for training has limited their practical applications. In this study, we combine a physically inspired dipole interaction model and an accurate neural network method for predicting the polarizability tensors of molecules. With the local chemical environment precisely described and the requirement of rotational covariance naturally fulfilled, this hybrid model is proven to give an accurate molecular polarizability prediction, essentially reducing the number of training samples. The atomic polarizabilities are physically interpretable and transferable to larger molecules unseen in the training set. This promising method may find its wide range of applications, such as spectroscopic simulations and the construction of polarizable force fields.

Transferability of the Electrostatic Parameters of the Polarizable Gaussian Multipole Model

Accuracy and transferability are the two highly desirable properties of molecular mechanical force fields. Compared with the extensively used point-charge additive force fields that apply fixed atom-centered point partial charges to model electrostatic interactions, polarizable force fields are thought to have the advantage of modeling the atomic polarization effects. Previous works have demonstrated the accuracy of the recently developed polarizable Gaussian multipole (pGM) models. In this work, we assessed the transferability of the electrostatic parameters of the pGM models with (pGM-perm) and without (pGM-ind) atomic permanent dipoles in terms of reproducing the electrostatic potentials surrounding molecules/oligomers absent from electrostatic parameterizations. Encouragingly, both the pGM-perm and pGM-ind models show significantly improved transferability than the additive model in the tests (1) from water monomer to water oligomer clusters; (2) across different conformations of amino acid dipeptides and tetrapeptides; (3) from amino acid tetrapeptides to longer polypeptides; and (4) from nucleobase monomers to Watson-Crick base pair dimers and tetramers. Furthermore, we demonstrated that the double-conformation fittings using amino acid tetrapeptides in the αR and β conformations can result in good transferability not only across different tetrapeptide conformations but also from tetrapeptides to polypeptides with lengths ranging from 1 to 20 repetitive residues for both the pGM-ind and pGM-perm models. In addition, the observation that the pGM-ind model has significantly better accuracy and transferability than the point-charge additive model, even though they have an identical number of parameters, strongly suggest the importance of intramolecular polarization effects. In summary, this and previous works together show that the pGM models possess both accuracy and transferability, which are expected to serve as foundations for the development of next-generation polarizable force fields for modeling various polarization-sensitive biological systems and processes.

Treating Polarization Effects in Charged and Polar Bio-Molecules Through Variable Electrostatic Parameters

Polarization plays important roles in charged and hydrogen bonding containing systems. Much effort ranging from the construction of physics-based models to quantum mechanism (QM)-based and machine learning (ML)-assisted models have been devoted to incorporating the polarization effect into the conventional force fields at different levels, such as atomic and coarse grained (CG). The application of polarizable force fields or polarization models was limited by two aspects, namely, computational cost and transferability. Different from physics-based models, no predetermining parameters were required in the QM-based approaches. Taking advantage of both the accuracy of QM calculations and efficiency of molecular mechanism (MM) and ML, polarization effects could be treated more efficiently while maintaining the QM accuracy. The computational cost could be reduced with variable electrostatic parameters, such as the charge, dipole, and electronic dielectric constant with the help of linear scaling fragmentation-based QM calculations and ML models. Polarization and entropy effects on the prediction of partition coefficient of druglike molecules are demonstrated by using both explicit or implicit all-atom molecular dynamics simulations and machine learning-assisted models. Directions and challenges for future development are also envisioned.

Computer Aided Development of Nucleic Acid Applications in Nanotechnologies

Utilization of nucleic acids (NAs) in nanotechnologies and nanotechnology-related applications is a growing field with broad application potential, ranging from biosensing up to targeted cell delivery. Computer simulations are useful techniques that can aid design and speed up development in this field. This review focuses on computer simulations of hybrid nanomaterials composed of NAs and other components. Current state-of-the-art molecular dynamics simulations, empirical force fields (FFs), and coarse-grained approaches for the description of deoxyribonucleic acid and ribonucleic acid are critically discussed. Challenges in combining biomacromolecular and nanomaterial FFs are emphasized. Recent applications of simulations for modeling NAs and their interactions with nano- and biomaterials are overviewed in the fields of sensing applications, targeted delivery, and NA templated materials. Future perspectives of development are also highlighted.

Improving Condensed-Phase Water Dynamics with Explicit Nuclear Quantum Effects: The Polarizable Q-AMOEBA Force Field

We introduce a new parametrization of the AMOEBA polarizable force field for water denoted Q-AMOEBA, for use in simulations that explicitly account for nuclear quantum effects (NQEs). This study is made possible thanks to the recently introduced adaptive Quantum Thermal Bath (adQTB) simulation technique which computational cost is comparable to classical molecular dynamics. The flexible Q-AMOEBA model conserves the initial AMOEBA functional form, with an intermolecular potential including an atomic multipole description of electrostatic interactions (up to quadrupole), a polarization contribution based on the Thole interaction model and a buffered 14-7 potential to model van der Waals interactions. It has been obtained by using a ForceBalance fitting strategy including high-level quantum chemistry reference energies and selected condensed-phase properties targets. The final Q-AMOEBA model is shown to accurately reproduce both gas-phase and condensed-phase properties, notably improving the original AMOEBA water model. This development allows the fine study of NQEs on water liquid phase properties such as the average H-O-H angle compared to its gas-phase equilibrium value, isotope effects, and so on. Q-AMOEBA also provides improved infrared spectroscopy prediction capabilities compared to AMOEBA03. Overall, we show that the impact of NQEs depends on the underlying model functional form and on the associated strength of hydrogen bonds. Since adQTB simulations can be performed at near classical computational cost using the Tinker-HP package, Q-AMOEBA can be extended to organic molecules, proteins, and nucleic acids opening the possibility for the large-scale study of the importance of NQEs in biophysics.

M-Chem: a Modular Software Package for Molecular Simulation that Spans Scientific Domains

We present a new software package called M-Chem that is designed from scratch in C++ and parallelized on shared-memory multi-core architectures to facilitate efficient molecular simulations. Currently, M-Chem is a fast molecular dynamics (MD) engine that supports the evaluation of energies and forces from two-body to many-body all-atom potentials, reactive force fields, coarse-grained models, combined quantum mechanics molecular mechanics (QM/MM) models, and external force drivers from machine learning, augmented by algorithms that are focused on gains in computational simulation times. M-Chem also includes a range of standard simulation capabilities including thermostats, barostats, multi-timestepping, and periodic cells, as well as newer methods such as fast extended Lagrangians and high quality electrostatic potential generation. At present M-Chem is a developer friendly environment in which we encourage new software contributors from diverse fields to build their algorithms, models, and methods in our modular framework. The long-term objective of M-Chem is to create an interdisciplinary platform for computational methods with applications ranging from biomolecular simulations, reactive chemistry, to materials research.

Accurate Reproduction of Quantum Mechanical Many-Body Interactions in Peptide Main-Chain Hydrogen-Bonding Oligomers by the Polarizable Gaussian Multipole Model

A key advantage of polarizable force fields is their ability to model the atomic polarization effects that play key roles in the atomic many-body interactions. In this work, we assessed the accuracy of the recently developed polarizable Gaussian Multipole (pGM) models in reproducing quantum mechanical (QM) interaction energies, many-body interaction energies, as well as the nonadditive and additive contributions to the many-body interactions for peptide main-chain hydrogen-bonding conformers, using glycine dipeptide oligomers as the model systems. Two types of pGM models were considered, including that with (pGM-perm) and without (pGM-ind) permanent atomic dipoles. The performances of the pGM models were compared with several widely used force fields, including two polarizable (Amoeba13 and ff12pol) and three additive (ff19SB, ff15ipq, and ff03) force fields. Encouragingly, the pGM models outperform all other force fields in terms of reproducing QM interaction energies, many-body interaction energies, as well as the nonadditive and additive contributions to the many-body interactions, as measured by the root-mean-square errors (RMSEs) and mean absolute errors (MAEs). Furthermore, we tested the robustness of the pGM models against polarizability parameterization errors by employing alternative polarizabilities that are either scaled or obtained from other force fields. The results show that the pGM models with alternative polarizabilities exhibit improved accuracy in reproducing QM many-body interaction energies as well as the nonadditive and additive contributions compared with other polarizable force fields, suggesting that the pGM models are robust against the errors in polarizability parameterizations. This work shows that the pGM models are capable of accurately modeling polarization effects and have the potential to serve as templates for developing next-generation polarizable force fields for modeling various biological systems.

Non-covalent interactions from a Quantum Chemical Topology perspective

About half a century after its little-known beginnings, the quantum topological approach called QTAIM has grown into a widespread, but still not mainstream, methodology of interpretational quantum chemistry. Although often confused in textbooks with yet another population analysis, be it perhaps an elegant but somewhat esoteric one, QTAIM has been enriched with about a dozen other research areas sharing its main mathematical language, such as Interacting Quantum Atoms (IQA) or Electron Localisation Function (ELF), to form an overarching approach called Quantum Chemical Topology (QCT). Instead of reviewing the latter’s role in understanding non-covalent interactions, we propose a number of ideas emerging from the full consequences of the space-filling nature of topological atoms, and discuss how they (will) impact on interatomic interactions, including non-covalent ones. The architecture of a force field called FFLUX, which is based on these ideas, is outlined. A new method called Relative Energy Gradient (REG) is put forward, which is able, by computation, to detect which fragments of a given molecular assembly govern the energetic behaviour of this whole assembly. This method can offer insight into the typical balance of competing atomic energies both in covalent and non-covalent case studies. A brief discussion on so-called bond critical points is given, highlighting concerns about their meaning, mainly in the arena of non-covalent interactions.

Development of the Quantum-Inspired SIBFA Many-Body Polarizable Force Field: Enabling Condensed-Phase Molecular Dynamics Simulations

We present the extension of the Sum of Interactions Between Fragments Ab initio Computed (SIBFA) many-body polarizable force field to condensed-phase molecular dynamics (MD) simulations. The quantum-inspired SIBFA procedure is grounded on simplified integrals obtained from localized molecular orbital theory and achieves full separability of its intermolecular potential. It embodies long-range multipolar electrostatics (up to quadrupole) coupled to a short-range penetration correction (up to charge-quadrupole), exchange repulsion, many-body polarization, many-body charge transfer/delocalization, exchange dispersion, and dispersion (up to C10). This enables the reproduction of all energy contributions of ab initio symmetry-adapted perturbation theory (SAPT(DFT)) gas-phase reference computations. The SIBFA approach has been integrated within the Tinker-HP massively parallel MD package. To do so, all SIBFA energy gradients have been derived and the approach has been extended to enable periodic boundary conditions simulations using smooth particle mesh Ewald. This novel implementation also notably includes a computationally tractable simplification of the many-body charge transfer/delocalization contribution. As a proof of concept, we perform a first computational experiment defining a water model fitted on a limited set of SAPT(DFT) data. SIBFA is shown to enable a satisfactory reproduction of both gas-phase energetic contributions and condensed-phase properties highlighting the importance of its physically motivated functional form.

Accurate Deep Learning-Aided Density-Free Strategy for Many-Body Dispersion-Corrected Density Functional Theory

Using a deep neuronal network (DNN) model trained on the large ANI-1 data set of small organic molecules, we propose a transferable density-free many-body dispersion (DNN-MBD) model. The DNN strategy bypasses the explicit Hirshfeld partitioning of the Kohn-Sham electron density required by MBD models to obtain the atom-in-molecules volumes used by the Tkatchenko-Scheffler polarizability rescaling. The resulting DNN-MBD model is trained with minimal basis iterative Stockholder atomic volumes and, coupled to density functional theory (DFT), exhibits comparable (if not greater) accuracy to other approaches based on different partitioning schemes. Implemented in the Tinker-HP package, the DNN-MBD model decreases the overall computational cost compared to MBD models where the explicit density partitioning is performed. Its coupling with the recently introduced Stochastic formulation of the MBD equations (J. Chem. Theory Comput. 2022, 18 (3), 1633-1645) enables large routine dispersion-corrected DFT calculations at preserved accuracy. Furthermore, the DNN electron density-free features extend the MBD model's applicability beyond electronic structure theory within methodologies such as force fields and neural networks.

Comparison of approximate intermolecular potentials for ab initio fragment calculations on medium sized N-heterocycles

The ground state intermolecular potential of bimolecular complexes of N‐heterocycles is analyzed for the impact of individual terms in the interaction energy as provided by various, conceptually different theories. Novel combinations with several formulations of the electrostatic, Pauli repulsion, and dispersion contributions are tested at both short‐ and long‐distance sides of the potential energy surface, for various alignments of the pyrrole dimer as well as the cytosine–uracil complex. The integration of a DFT/CCSD density embedding scheme, with dispersion terms from the effective fragment potential (EFP) method is found to provide good agreement with a reference CCSD(T) potential overall; simultaneously, a quantum mechanics/molecular mechanics approach using CHELPG atomic point charges for the electrostatic interaction, augmented by EFP dispersion and Pauli repulsion, comes also close to the reference result. Both schemes have the advantage of not relying on predefined force fields; rather, the interaction parameters can be determined for the system under study, thus being excellent candidates for ab initio modeling.

O(N) Stochastic Evaluation of Many-Body van der Waals Energies in Large Complex Systems

We propose a new strategy to solve the key equations of the many-body dispersion (MBD) model by Tkatchenko, DiStasio Jr., and Ambrosetti. Our approach overcomes the original O(N3) computational complexity that limits its applicability to large molecular systems within the context of O(N) density functional theory. First, to generate the required frequency-dependent screened polarizabilities, we introduce an efficient solution to the Dyson-like self-consistent screening equations. The scheme reduces the number of variables and, coupled to a direct inversion of the iterative subspace extrapolation, exhibits linear-scaling performances. Second, we apply a stochastic Lanczos trace estimator resolution to the equations evaluating the many-body interaction energy of coupled quantum harmonic oscillators. While scaling linearly, it also enables communication-free pleasingly parallel implementations. As the resulting O(N) stochastic massively parallel MBD approach is found to exhibit minimal memory requirements, it opens up the possibility of computing accurate many-body van der Waals interactions of millions-atoms' complex materials and solvated biosystems with computational times in the range of minutes.

Reliable Prediction of the Protein-Ligand Binding Affinity Using a Charge Penetration Corrected AMOEBA Force Field: A Case Study of Drug Resistance Mutations in Abl Kinase

Protein mutations that directly impair drug binding are related to therapeutic resistance, and accurate prediction of their impact on drug binding would benefit drug design and clinical practice. Here, we have developed a scoring strategy that predicts the effect of the mutations on the protein-ligand binding affinity. In view of the critical importance of electrostatics in protein-ligand interactions, the charge penetration corrected AMOEBA force field (AMOEBA_CP model) was employed to improve the accuracy of the calculated electrostatic energy. We calculated the electrostatic energy using an energy decomposition analysis scheme based on the generalized Kohn-Sham (GKS-EDA). The AMOEBA_CP model was validated by a protein-fragment-ligand complex data set (Abl236) constructed from the co-crystal structures of the cancer target Abl kinase with six inhibitors. To predict ligand binding affinity changes upon protein mutation of Abl kinase, we used sampling protocol with multistep simulated annealing to search conformations of mutant proteins. The scoring strategy based on AMOEBA_CP model has achieved considerable performance in predicting resistance for 8 kinase inhibitors across 144 clinically identified point mutations. Overall, this study illustrates that the AMOEBA_CP model, which accurately treats electrostatics through penetration correction, enables the accurate prediction of the mutation-induced variation of protein-ligand binding affinity.

An Efficient Gaussian-Accelerated Molecular Dynamics (GaMD) Multilevel Enhanced Sampling Strategy: Application to Polarizable Force Fields Simulations of Large Biological Systems

We introduce a novel multilevel enhanced sampling strategy grounded on Gaussian-accelerated Molecular Dynamics (GaMD). First, we propose a GaMD multi-GPUs-accelerated implementation within the Tinker-HP molecular dynamics package. We introduce the new "dual-water" mode and its use with the flexible AMOEBA polarizable force field. By adding harmonic boosts to the water stretching and bonding terms, it accelerates the solvent-solute interactions while enabling speedups, thanks to the use of fast multiple-time step integrators. To further reduce the time-to-solution, we couple GaMD to Umbrella Sampling (US). The GaMD─US/dual-water approach is tested on the 1D Potential of Mean Force (PMF) of the solvated CD2-CD58 system (168 000 atoms), allowing the AMOEBA PMF to converge within 1 kcal/mol of the experimental value. Finally, Adaptive Sampling (AS) is added, enabling AS-GaMD capabilities but also the introduction of the new Adaptive Sampling-US-GaMD (ASUS-GaMD) scheme. The highly parallel ASUS-GaMD setup decreases time to convergence by, respectively, 10 and 20 times, compared to GaMD-US and US. Overall, beside the acceleration of PMF computations, Tinker-HP now allows for the simultaneous use of Adaptive Sampling and GaMD-"dual water" enhanced sampling approaches increasing the applicability of polarizable force fields to large-scale simulations of biological systems.

Atomic Polarizabilities for Interactive Dipole Induction Models

Thole-style mutual induction models for molecular polarization have been adopted by several popular polarizable force fields (FFs) for their simplicity and transferability. The atomic polarizability parameters of these models are typically derived by fitting to ab initio or/and experimental molecular polarizabilities. In this work, we improve upon Thole polarizability parameters by employing both high-level quantum mechanics molecular polarizabilities and electrostatic potential (ESP) responses on three-dimensional grids. Our results indicate that the two approaches to derive atomic polarizability parameters are both effective, while the ESP approaches can also capture the polarization for the atoms with lone pair electrons. The resulting polarizability parameters have been validated on a set of over 7200 molecules covering the most common elements found in organic molecules (C, H, O, N, P, S, F, Cl, Br, and I). These parameters have also been tested on the experimentally measured molecular polarizabilities of 422 molecules. The final set of parameters derived in this work show notable improvement over the current AMOEBA set. The result is a highly transferable, expanded set of atomic polarizabilities defined by the local chemical environment in the form of SMARTS patterns. These parameters can be used directly in molecular mechanics polarizable potential energy functions such as AMOEBA, AMOEBA+, and other Thole-style models.

Short-Range Imbalances in the AMBER Lennard-Jones Potential for (Deoxy)Ribose···Nucleobase Lone-Pair···π Contacts in Nucleic Acids

The lone-pair···π (lp···π) (deoxy)ribose···nucleobase stacking is a recurring interaction in Z-DNA and RNAs that is characterized by sub-van der Waals lp···π contacts (<3.0 Å). It is a part of the structural signature of CpG Z-step motifs in Z-DNA and r(UNCG) tetraloops that are known to behave poorly in molecular dynamics (MD) simulations. Although the exact origin of the MD simulation issues remains unclear, a significant part of the problem might be due to an imbalanced description of nonbonded interactions, including the characteristic lp···π stacking. To gain insights into the links between lp···π stacking and MD, we present an in-depth comparison between accurate large-basis-set double-hybrid Kohn-Sham density functional theory calculations DSD-BLYP-D3/ma-def2-QZVPP (DHDF-D3) and data obtained with the nonbonded potential of the AMBER force field (AFF) for NpN Z-steps (N = G, A, C, and U). Among other differences, we found that the AFF overestimates the DHDF-D3 lp···π distances by ∼0.1-0.2 Å, while the deviation between the DHDF-D3 and AFF descriptions sharply increases in the short-range region of the interaction. Based on atom-in-molecule polarizabilities and symmetry-adapted perturbation theory analysis, we inferred that the DHDF-D3 versus AFF differences partly originate in identical nucleobase carbon atom Lennard-Jones (LJ) parameters despite the presence/absence of connected electron-withdrawing groups that lead to different effective volumes or vdW radii. Thus, to precisely model the very short CpG lp···π contact distances, we recommend revision of the nucleobase atom LJ parameters. Additionally, we suggest that the large discrepancy between DHDF-D3 and AFF short-range repulsive part of the interaction energy potential may significantly contribute to the poor performances of MD simulations of nucleic acid systems containing Z-steps. Understanding where, and if possible why, the point-charge-type effective potentials reach their limits is vital for developing next-generation FFs and for addressing specific issues in contemporary MD simulations.

Improvement of the Gaussian Electrostatic Model by separate fitting of Coulomb and exchange-repulsion densities and implementation of a new dispersion term

The description of each separable contribution of the intermolecular interaction is a useful approach to develop polarizable force fields (polFFs). The Gaussian Electrostatic Model (GEM) is based on this approach, coupled with the use of density fitting techniques. In this work, we present the implementation and testing of two improvements of GEM: the Coulomb and exchange-repulsion energies are now computed with separate frozen molecular densities and a new dispersion formulation inspired by the Sum of Interactions Between Fragments Ab initio Computed polFF, which has been implemented to describe the dispersion and charge-transfer interactions. Thanks to the combination of GEM characteristics and these new features, we demonstrate a better agreement of the computed structural and condensed properties for water with experimental results, as well as binding energies in the gas phase with the ab initio reference compared with the previous GEM* potential. This work provides further improvements to GEM and the items that remain to be improved and the importance of the accurate reproduction for each separate contribution.

Catalytic Fields as a Tool to Analyze Enzyme Reaction Mechanism Variants and Reaction Steps

Catalytic fields representing the topology of the optimal molecular environment charge distribution that reduces the activation barrier have been used to examine alternative reaction variants and to determine the role of conserved catalytic residues for two consecutive reactions catalyzed by the same enzyme. Until now, most experimental and conventional top-down theoretical studies employing QM/MM or ONIOM methods have focused on the role of enzyme electric fields acting on broken bonds of reactants. In contrast, our bottom-up approach dealing with a small reactant and transition-state model allows the analysis of the opposite effects: how the catalytic field resulting from the charge redistribution during the enzyme reaction acts on conserved amino acid residues and contributes to the reduction of the activation barrier. This approach has been applied to the family of histidyl tRNA synthetases involved in the translation of the genetic code into the protein amino acid sequence. Activation energy changes related to conserved charged amino acid residues for 12 histidyl tRNA synthetases from different biological species allowed to compare on equal footing the catalytic residues involved in ATP aminoacylation and tRNA charging reactions and to analyze different reaction mechanisms proposed in the literature. A scan of the library of atomic multipoles for amino acid side-chain rotamers within the catalytic field pointed out the change in the Glu83 conformation as the critical catalytic effect, providing, at low computational cost, insight into the electrostatic preorganization of the enzyme catalytic site at a level of detail that has not yet been accessible in conventional experimental or theoretical methods. This opens the way for rational reverse biocatalyst design at a very limited computational cost without resorting to empirical methods.

A Quantum Chemical Topology Picture of Intermolecular Electrostatic Interactions and Charge Penetration Energy

Based on the Interacting Quantum Atoms approach, we present herein a conceptual and theoretical framework of short-range electrostatic interactions, whose accurate description is still a challenging problem in molecular modeling. For all the noncovalent complexes in the S66 database, the fragment-based and atomic decomposition of the electrostatic binding energies is performed using both the charge density of the dimers and the unrelaxed densities of the monomers. This energy decomposition together with dispersion corrections gives rise to a pairwise approximation to the total binding energy. It also provides energetic descriptors at varying distance that directly address the atomic and molecular electrostatic interactions as described by point-charge or multipole-based potentials. Additionally, we propose a consistent definition of the charge penetration energy within quantum chemical topology, which is mainly characterized in terms of the intramolecular electrostatic energy. Finally, we discuss some practical implications of our results for the design and validation of electrostatic potentials.

Atomistic Polarizable Embeddings: Energy, Dynamics, Spectroscopy, and Reactivity

The computational modeling of realistic extended systems, relevant in, e.g., Chemistry and Biophysics, is a fundamental problem of paramount importance in contemporary research. Enzymatic catalysis and photoinduced processes in pigment-protein complexes are typical problems targeted by computer-aided approaches, to complement experiments as interpretative tools at a molecular scale. The daunting complexity of this task lies in between the opposite stringent requirements of results' reliability for structural/dynamical properties and related intermolecular interactions, and a mandatory principle of realism in the modeling strategy. Therefore, in practice, a truly realistic computational model of a biologically relevant system can easily fail to meet the accuracy requirement, in order to balance the excessive computational cost necessary to reach the desired precision.To address such an "accuracy vs reality" dualistic requirement, mixed quantum mechanics/classical mechanics approaches within Atomistic (i.e., preserving the discrete particle configuration) Polarizable Embeddings (QM/APEs) methods have been proposed over the years. In this Account, we review recent developments in the design and application of general QM/APE methods, targeting situations where a local intrinsically quantum behavior is coupled to a large molecular system (i.e., an environment), often involving processes with different dynamical time scales, in order to avoid brute-force, unpractical quantum chemistry calculations on the complete system.In the first place, our interest is devoted to the available APEs models presently implemented in computational software, highlighting the quantum chemistry methods that can be used to treat the QM subsystem. We review the coupling strategy between the QM subsystem and the APE, which requires to examine the way the QM/MM mutual interactions are accounted for and how the polarization of the classical environment is considered with respect to (wrt) the quantum variables. Because of the need of reliable molecular and macromolecular structures, a pivotal aspect to address here is the handling of the system dynamics (i.e., gradients wrt nuclear positions are required), especially for large molecular assemblies composed by an overwhelming number of atoms, exploring many conformations on a complex energy landscape.Alongside, we highlight our views on the necessary steps to take toward more accurate general-purposes and transferable explicit embeddings. The main objective to achieve here is to design a more physically grounded multiscale approach. To do so, one should apply advanced new generation classical models to account for refined induction effects that are able to (i) improve the quality of QM/MM interaction energies; (ii) enhance transferability by avoiding the compulsory partial (or total) reparameterization of the classical model. Moreover, the extension of recent developments originating from the field of advanced classical molecular dynamics (MD) to the realm of QM/APE methods is a key direction to improve both speed and efficiency for the phase space exploration of systems of growing size and complexity.Lastly, we point out specific research topics where an advanced QM/APE dynamics can certainly shed some light. For example, we discuss chemical reactions in "harsh" environments and the case of spectroscopic theoretical modeling where the inclusion of refined environment effects is often mandatory.

Development and application of quantum mechanics/molecular mechanics methods with advanced polarizable potentials

Quantum mechanics/molecular mechanics (QM/MM) simulations are a popular approach to study various features of large systems. A common application of QM/MM calculations is in the investigation of reaction mechanisms in condensed-phase and biological systems. The combination of QM and MM methods to represent a system gives rise to several challenges that need to be addressed. The increase in computational speed has allowed the expanded use of more complicated and accurate methods for both QM and MM simulations. Here, we review some approaches that address several common challenges encountered in QM/MM simulations with advanced polarizable potentials, from methods to account for boundary across covalent bonds and long-range effects, to polarization and advanced embedding potentials.

Multimolecular complexes of the phosphodiester anion with Zn(II) or Mg(II) and water molecules-Preliminary validations of a polarizable potential by ab initio quantum chemistry

Dimethyl phosphate (DMP^- ) is a model for the phosphodiester backbone of DNA, RNA, and phospholipids. It is central for the binding of divalent cations and water along the backbone of nucleic acids. Significant polarization and charge-transfer contributions and nonadditivity come into play in the multimolecular complexes organized around phosphate. Prior to large-scale molecular dynamics (MD) with advanced polarizable potentials, it is essential to evaluate how well the values and trends of intermolecular interaction energies (ΔE) from ab initio quantum chemistry (QC) and their individual contributions are reproduced in a diversity of such complexes. These differ by the starting binding modes of a divalent cation, Zn(II), namely direct, bi- or mono-dentate to anionic and/or ester oxygens, versus through-water binding. We present first the results from automated refinements of the individual contributions of the SIBFA potential with respect to their QC counterparts using a Zn(II) or a water probe. This is followed by validations on eight relaxed multimolecular complexes of DMP^- with Zn(II) or Mg(II) and seven waters, then on sixteen complexes of DMP^- with Zn(II) and eight waters in arrangements extracted from MD or energy-minimization on a droplet of sixty-four waters. This monitors the compared evolutions of SIBFA and QC ΔE and their individual contributions in the competing arrangements. Some waters, bridging Zn(II) and DMP^- , were found to have exceptionally large dipole moments, of up to 3.8 Debye. The perspectives of extension to a flexible phosphodiester backbone are discussed in the context of the SIBFA potential for DNA and RNA.

CLIFF: A component-based, machine-learned, intermolecular force field

Computation of intermolecular interactions is a challenge in drug discovery because accurate ab initio techniques are too computationally expensive to be routinely applied to drug-protein models. Classical force fields are more computationally feasible, and force fields designed to match symmetry adapted perturbation theory (SAPT) interaction energies can remain accurate in this context. Unfortunately, the application of such force fields is complicated by the laborious parameterization required for computations on new molecules. Here, we introduce the component-based machine-learned intermolecular force field (CLIFF), which combines accurate, physics-based equations for intermolecular interaction energies with machine-learning models to enable automatic parameterization. The CLIFF uses functional forms corresponding to electrostatic, exchange-repulsion, induction/polarization, and London dispersion components in SAPT. Molecule-independent parameters are fit with respect to SAPT2+(3)δMP2/aug-cc-pVTZ, and molecule-dependent atomic parameters (atomic widths, atomic multipoles, and Hirshfeld ratios) are obtained from machine learning models developed for C, N, O, H, S, F, Cl, and Br. The CLIFF achieves mean absolute errors (MAEs) no worse than 0.70 kcal mol^-1 in both total and component energies across a diverse dimer test set. For the side chain-side chain interaction database derived from protein fragments, the CLIFF produces total interaction energies with an MAE of 0.27 kcal mol^-1 with respect to reference data, outperforming similar and even more expensive methods. In applications to a set of model drug-protein interactions, the CLIFF is able to accurately rank-order ligand binding strengths and achieves less than 10% error with respect to SAPT reference values for most complexes.

Tinker-HP: Accelerating Molecular Dynamics Simulations of Large Complex Systems with Advanced Point Dipole Polarizable Force Fields Using GPUs and Multi-GPU Systems

We present the extension of the Tinker-HP package (Lagardère, Chem. Sci. 2018, 9, 956-972) to the use of Graphics Processing Unit (GPU) cards to accelerate molecular dynamics simulations using polarizable many-body force fields. The new high-performance module allows for an efficient use of single- and multiple-GPU architectures ranging from research laboratories to modern supercomputer centers. After detailing an analysis of our general scalable strategy that relies on OpenACC and CUDA, we discuss the various capabilities of the package. Among them, the multiprecision possibilities of the code are discussed. If an efficient double precision implementation is provided to preserve the possibility of fast reference computations, we show that a lower precision arithmetic is preferred providing a similar accuracy for molecular dynamics while exhibiting superior performances. As Tinker-HP is mainly dedicated to accelerate simulations using new generation point dipole polarizable force field, we focus our study on the implementation of the AMOEBA model. Testing various NVIDIA platforms including 2080Ti, 3090, V100, and A100 cards, we provide illustrative benchmarks of the code for single- and multicards simulations on large biosystems encompassing up to millions of atoms. The new code strongly reduces time to solution and offers the best performances to date obtained using the AMOEBA polarizable force field. Perspectives toward the strong-scaling performance of our multinode massive parallelization strategy, unsupervised adaptive sampling and large scale applicability of the Tinker-HP code in biophysics are discussed. The present software has been released in phase advance on GitHub in link with the High Performance Computing community COVID-19 research efforts and is free for Academics (see https://github.com/TinkerTools/tinker-hp).

Parameterization of Monovalent Ions for the OPC3, OPC, TIP3P-FB, and TIP4P-FB Water Models

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.

Polarizable Multipolar Molecular Dynamics Using Distributed Point Charges

Distributed point charge models (DCM) and their minimal variants (MDCM) have been integrated with tools widely used for condensed-phase simulations, including a virial-based barostat and a slow-growth algorithm for thermodynamic integration. Minimal DCM is further developed in a systematic fashion to reduce fitting errors in the electrostatic interaction energy, and a new fragment-based approach offers considerable speedup of the MDCM fitting process for larger molecules with increased numbers of off-centered charged sites. Finally, polarizable (M)DCM is also introduced in the present work. The developments are used in condensed-phase simulations of popular force fields with commonly applied simulation conditions. (M)DCM equivalents for a range of widely used water force fields and for fluorobenzene (PhF) are developed and applied along with the original models to evaluate the impact of reformulating the electrostatic term. Comparisons of the molecular electrostatic potential (MEP), electrostatic interaction energies, and bulk properties from molecular dynamics simulations for a range of models from simple TIPnP (n = 3-5) to the polarizable, multipolar iAMOEBA models for water and an existing quadrupolar model for PhF confirm that DCMs retain the accuracy of the original models, providing a homogeneous, efficient, and generic point charge alternative to a multipolar electrostatic model for force field development and multilevel simulations.

Effective Fragment Potentials for Flexible Molecules: Transferability of Parameters and Amino Acid Database

An accurate but efficient description of noncovalent interactions is a key to predictive modeling of biological and materials systems. The effective fragment potential (EFP) is an ab initio-based force field that provides a physically meaningful decomposition of noncovalent interactions of a molecular system into Coulomb, polarization, dispersion, and exchange-repulsion components. An EFP simulation protocol consists of two steps, preparing parameters for molecular fragments by a series of ab initio calculations on each individual fragment, and calculation of interaction energy and properties of a total molecular system based on the prepared parameters. As the fragment parameters (distributed multipoles, polarizabilities, localized wave function, etc.) depend on a fragment geometry, straightforward application of the EFP method requires recomputing parameters of each fragment if its geometry changes, for example, during thermal fluctuations of a molecular system. Thus, recomputing fragment parameters can easily become both computational and human bottlenecks and lead to a loss of efficiency of a simulation protocol. An alternative approach, in which fragment parameters are adjusted to different fragment geometries, referred to as "flexible EFP", is explored here. The parameter adjustment is based on translations and rotations of local coordinate frames associated with fragment atoms. The protocol is validated on extensive benchmark of amino acid dimers extracted from molecular dynamics snapshots of a cryptochrome protein. A parameter database for standard amino acids is developed to automate flexible EFP simulations in proteins. To demonstrate applicability of flexible EFP in large-scale protein simulations, binding energies and vertical electron ionization and electron attachment energies of a lumiflavin chromophore of the cryptochrome protein are computed. The results obtained with flexible EFP are in a close agreement with the standard EFP procedure but provide a significant reduction in computational cost.

An Iterative Fragment Scheme for the ACKS2 Electronic Polarization Model: Application to Molecular Dimers and Chains

The treatment of electrostatic interactions is a key ingredient in the force field-based simulation of condensed phase systems. Most approaches used fixed, site-specific point charges. Yet, it is now clear that many applications of force fields (FFs) demand more sophisticated treatments, prompting the implementation of charge equilibration methods in polarizable FFs to allow the redistribution of charge within the system. One approach allowing both, charge redistribution and site-specific polarization, while at the same time solving methodological shortcomings of earlier methods, is the first-principles-derived atom-condensed Kohn-Sham density functional theory method approximated to the second order (ACKS2). In this work, we present two fragment approaches to ACKS2, termed f-ACKS2 and a self-consistent version, scf-ACKS2, that treat condensed phase systems as a collection of electronically polarizable molecular fragments. The fragmentation approach to ACKS2 not only leads to a more transferable and less system-specific collection of electronic response parameters but also opens up the method to large condensed phase systems. We validate the accuracies of f-ACKS2 and scf-ACKS2 by comparing polarization energies and induced dipole moments for a number of charged hydrocarbon dimers against DFT reference calculations. Finally, we also apply both fragmented ACKS2 variants to calculate the polarization energy for electron-hole pair separation along a chain of anthracene molecules and find excellent agreement with reference DFT calculations.

Quantum-Chemistry Based Design of Halobenzene Derivatives With Augmented Affinities for the HIV-1 Viral G4/C16 Base-Pair

The HIV-1 integrase (IN) is a major target for the design of novel anti-HIV inhibitors. Among these, three inhibitors which embody a halobenzene ring derivative (HR) in their structures are presently used in clinics. High-resolution X-ray crystallography of the complexes of the IN-viral DNA transient complex bound to each of the three inhibitors showed in all cases the HR ring to interact within a confined zone of the viral DNA, limited to the highly conserved 5′CpA 3′/5′TpG 3′ step. The extension of its extracyclic CX bond is electron-depleted, owing to the existence of the “sigma-hole.” It interacts favorably with the electron-rich rings of base G₄. We have sought to increase the affinity of HR derivatives for the G₄/C₁₆ base pair. We thus designed thirteen novel derivatives and computed their Quantum Chemistry (QC) intermolecular interaction energies (ΔE) with this base-pair. Most compounds had ΔE values significantly more favorable than those of the HR of the most potent halobenzene drug presently used in clinics, Dolutegravir. This should enable the improvement in a modular piece-wise fashion, the affinities of halogenated inhibitors for viral DNA (vDNA). In view of large scale polarizable molecular dynamics simulations on the entirety of the IN-vDNA-inhibitor complexes, validations of the SIBFA polarizable method are also reported, in which the evolution of each ΔE(SIBFA) contribution is compared to its QC counterpart along this series of derivatives.

Further Optimization and Validation of the Classical Drude Polarizable Protein Force Field

The CHARMM Drude-2013 polarizable force field (FF) was developed to include the explicit treatment of induced electronic polarizability, resulting in a more accurate description of the electrostatic interactions in molecular dynamics (MD) simulations. While the Drude-2013 protein FF has shown success in improving the folding properties of α-helical peptides and to reproduce experimental observables in simulations up to 1 μs, some limitations were noted regarding the stability of β-sheet structures in simulations longer than 100 ns as well as larger deviations from crystal structures in simulations of a number of proteins compared to the additive CHARMM36 protein FF. The origin of the instability has been identified and appears to be primarily due to overestimated atomic polarizabilities and induced dipole-dipole interactions on the Cβ, Cγ, and Cδ side chain atoms. To resolve this and other issues, a number of aspects of the model were revisited, resulting in Drude-2019 protein FF. Backbone parameters were optimized targeting the conformational properties of the (Ala)5 peptide in solution along with gas phase properties of the alanine dipeptide. Dipeptides that contain N-acetylated and N'-methylamidated termini, excluding Gly, Pro, and Ala, were used as models to optimize the atomic polarizabilities and Thole screening factors on selected Cβ, Cγ, and Cδ carbons by targeting quantum mechanical (QM) dipole moments and molecular polarizabilities. In addition, to obtain better conformational properties, side chain χ1 and χ2 dihedral parameters were optimized targeting QM data for the respective side chain dipeptide conformations as well as Protein Data Bank survey data based on the χ1, χ2 sampling from Hamiltonian replica-exchange MD simulations of (Ala)4-X-(Ala)4 in solution, where X is the amino acid of interest. Further improvements include optimizing nonbonded interactions between charged residues to reproduce QM interaction energies of the charged-protein model compounds and experimental osmotic pressures. Validation of the optimized Drude protein FF includes MD simulations of a collection of peptides and proteins including β-sheet structures, as well as transmembrane ion channels. Results showed that the updated Drude-2019 protein FF yields smaller overall root-mean-square differences of proteins as compared to the additive CHARMM36m and Drude-2013 FFs as well as similar or improved agreement with experimental NMR properties, allowing for long time scale simulation studies of proteins and more complex biomolecular systems in conjunction with the remainder of the Drude polarizable FF.

Thermodynamics of Transition Metal Ion Binding to Proteins

Modeling the thermodynamics of a transition metal (TM) ion assembly be it in proteins or in coordination complexes affords us a better understanding of the assembly and function of metalloclusters in diverse application areas including metal organic framework design, TM-based catalyst design, the trafficking of TM ions in biological systems, and drug design in metalloprotein platforms. While the structural details of TM ions bound to metalloproteins are generally well understood via experimental and computational approaches, accurate studies describing the thermodynamics of TM ion binding are rare. Herein, we demonstrate that we can obtain accurate structural and absolute binding free energies of Co²⁺ and Ni²⁺ to the enzyme glyoxalase I using an optimized 12-6-4 (m12-6-4) potential. Critically, this model simultaneously reproduces the solvation free energy of the individual TM ions and reproduces the thermodynamics of TM ion-ligand coordination as well as the thermodynamics of TM ion binding to a protein active site unlike extant models. We find the incorporation of the thermodynamics associated with protonation state changes for the TM ion (un)binding to be crucial. The high accuracy of m12-6-4 potential in this study presents an accurate route to explore more complicated processes associated with TM cluster assembly and TM ion transport.

Embedding Methods for Quantum Chemistry: Applications from Materials to Life Sciences

Quantum mechanical embedding methods hold the promise to transform not just the way calculations are performed, but to significantly reduce computational costs and improve scaling for macro-molecular systems containing hundreds if not thousands of atoms. The field of embedding has grown increasingly broad with many approaches of different intersecting flavors. In this perspective, we lay out the methods into two streams: QM:MM and QM:QM, showcasing the advantages and disadvantages of both. We provide a review of the literature, the underpinning theories including our contributions, and we highlight current applications with select examples spanning both materials and life sciences. We conclude with prospects and future outlook on embedding, and our view on the use of universal test case scenarios for cross-comparisons of the many available (and future) embedding theories.