We analyze convergence of energies and forces for the AMOEBA classical polarizable model when evaluated as a many-body expansion (MBE) against the corresponding N-body parent potential in the context of a condensed-phase water simulation. This is in contrast to most MBE formulations based on quantum mechanics, which focus only on convergence of energies for gas-phase clusters. Using a single water molecule as a definition of a body, we find that truncation of the MBE at third order, 3-AMOEBA, captures direct polarization exactly and yields apparent good convergence of the mutual polarization energy. However, it renders large errors in the magnitude of polarization forces and requires at least fourth-order terms in the MBE to converge toward the parent potential gradient values. We can improve the convergence of polarization forces for 3-AMOEBA by embedding the polarization response of dimers and trimers within a complete representation of the fixed electrostatics of the entire system. We show that the electrostatic embedding formalism helps identify the specific configurations involving linear hydrogen-bonding arrangements that are poorly convergent at the 3-body level. By extending the definition of a body to be a large water cluster, we can reduce errors in forces to yield an approximate polarization model that is up to 10 times faster than the parent potential. The 3-AMOEBA model offers new ways to investigate how the properties of bulk water depend on the degree of connectivity in the liquid.

In this work we present a replicated data hybrid OpenMP/MPI implementation of a hierarchical progression of approximate classical polarizable models that yields speedups of up to ∼10 compared to the standard OpenMP implementation of the exact parent AMOEBA polarizable model. In addition, our parallel implementation exhibits reasonable weak and strong scaling. The resulting parallel software will prove useful for those who are interested in how molecular properties converge in the condensed phase with respect to the MBE, it provides a fruitful test bed for exploring different electrostatic embedding schemes, and offers an interesting possibility for future exascale computing paradigms.

We have adapted a hybrid extended Lagrangian self-consistent field (EL/SCF) approach, developed for time reversible Born Oppenheimer molecular dynamics for quantum electronic degrees of freedom, to the problem of classical polarization. In this context, the initial guess for the mutual induction calculation is treated by auxiliary induced dipole variables evolved via a time-reversible velocity Verlet scheme. However, we find numerical instability, which is manifested as an accumulation in the auxiliary velocity variables, that in turn results in an unacceptable increase in the number of SCF cycles to meet even loose convergence tolerances for the real induced dipoles over the course of a 1 ns trajectory of the AMOEBA14 water model. By diagnosing the numerical instability as a problem of resonances that corrupt the dynamics, we introduce a simple thermostating scheme, illustrated using Berendsen weak coupling and Nose-Hoover chain thermostats, applied to the auxiliary dipole velocities. We find that the inertial EL/SCF (iEL/SCF) method provides superior energy conservation with less stringent convergence thresholds and a correspondingly small number of SCF cycles, to reproduce all properties of the polarization model in the NVT and NVE ensembles accurately. Our iEL/SCF approach is a clear improvement over standard SCF approaches to classical mutual induction calculations and would be worth investigating for application to ab initio molecular dynamics as well.

Computational modeling at the atomistic and mesoscopic levels has undergone dramatic development in the past 10 years to meet the challenge of adequately accounting for the many-body nature of intermolecular interactions. At the heart of this challenge is the ability to identify the strengths and specific limitations of pairwise-additive interactions, to improve classical models to explicitly account for many-body effects, and consequently to enhance their ability to describe a wider range of reference data and build confidence in their predictive capacity. However, the corresponding computational cost of these advanced classical models increases significantly enough that statistical convergence of condensed phase observables becomes more difficult to achieve. Here we review a hierarchy of potential energy surface models used in molecular simulations for systems with many degrees of freedom that best meet the trade-off between accuracy and computational speed in order to define a sweet spot for a given scientific problem of interest.

Molecular simulations of water using classical, molecular mechanic potential energy functions have enjoyed a 50-year history of development, and much has been learned regarding their parameterization and the essential physics that must be captured in order to reproduce water properties across the phase diagram and across system sizes, from the dimer to the condensed phase. While pairwise-additive force fields using fixed, point charge-based electrostatics have dominated this history owing to computational cost, their limitations in transferability are being recognized, owing particularly to the lack of many-body effects, as well as an inherent difficulty in capturing quantum mechanical effects that become important at short intermolecular separation. This has spurred an impressive development of novel functional forms and parameterization schemes to account for such effects, especially the leading many-body effect of polarization. This review discusses recent efforts in the development of advanced models of water, particularly with regard to important details of their parameterization from quantum mechanical or experimental data, the development of novel functional forms including machine learning-based models, and algorithms that reduce the computational cost of polarization dramatically, permitting them to potentially become competitive with pairwise-additive models as the standby of condensed-phase simulation. These technical developments are appraised based on their ability to impact numerical calculations on water, particularly the condensed phase, and it is hoped that this article provides a clear connection between the essential physics captured by the model and their fitness across a range of environments. WIREs Comput Mol Sci 2018, 8:e1355. doi: 10.1002/wcms.1355. This article is categorized under: Computer and Information Science > Computer Algorithms and Programming Molecular and Statistical Mechanics > Molecular Interactions Software > Simulation Methods.

AMOEBA is a molecular mechanics force field that addresses some of the shortcomings of a fixed partial charge model, by including permanent atomic point multipoles through quadrupoles, as well as many-body polarization through the use of point inducible dipoles. In this work, we investigate how well AMOEBA formulates its non-bonded interactions, and how it implicitly incorporates quantum mechanical effects such as charge penetration (CP) and charge transfer (CT), for water-water and water-ion interactions. We find that AMOEBA's total interaction energies, as a function of distance and over angular scans for the water dimer and for a range of water-monovalent cations, agree well with an advanced density functional theory (DFT) model, whereas the water-halides and water-divalent cations show significant disagreement with the DFT result, especially in the compressed region when the two fragments overlap. We use a second-generation energy decomposition analysis (EDA) scheme based on absolutely localized molecular orbitals (ALMOs) to show that in the best cases AMOEBA relies on cancellation of errors by softening of the van der Waals (vdW) wall to balance permanent electrostatics that are too unfavorable, thereby compensating for the missing CP effect. CT, as another important stabilizing effect not explicitly taken into account in AMOEBA, is also found to be incorporated by the softened vdW interaction. For the water-halides and water-divalent cations, this compensatory approach is not as well executed by AMOEBA over all distances and angles, wherein permanent electrostatics remains too unfavorable and polarization is overdamped in the former while overestimated in the latter. We conclude that the DFT-based EDA approach can help refine a next-generation AMOEBA model that either realizes a better cancellation of errors for problematic cases like those illustrated here, or serves to guide the parametrization of explicit functional forms for short-range contributions from CP and/or CT.

In this work, we evaluate the accuracy of the classical AMOEBA model for representing many-body interactions, such as polarization, charge transfer, and Pauli repulsion and dispersion, through comparison against an energy decomposition method based on absolutely localized molecular orbitals (ALMO-EDA) for the water trimer and a variety of ion-water systems. When the 2- and 3-body contributions according to the many-body expansion are analyzed for the ion-water trimer systems examined here, the 3-body contributions to Pauli repulsion and dispersion are found to be negligible under ALMO-EDA, thereby supporting the validity of the pairwise-additive approximation in AMOEBA's 14-7 van der Waals term. However AMOEBA shows imperfect cancellation of errors for the missing effects of charge transfer and incorrectness in the distance dependence for polarization when compared with the corresponding ALMO-EDA terms. We trace the larger 2-body followed by 3-body polarization errors to the Thole damping scheme used in AMOEBA, and although the width parameter in Thole damping can be changed to improve agreement with the ALMO-EDA polarization for points about equilibrium, the correct profile of polarization as a function of intermolecular distance cannot be reproduced. The results suggest that there is a need for re-examining the damping and polarization model used in the AMOEBA force field and provide further insights into the formulations of polarizable force fields in general.

Carbon dioxide (CO2) is a detrimental greenhouse gas and is the main contributor to global warming. In addressing this environmental challenge, a promising approach emerges through the utilization of deep eutectic solvents (DESs) as an ecofriendly and sustainable medium for effective CO2 capture. Chemically reactive DESs, which form chemical bonds with the CO2, are superior to nonreactive, physically based DESs for CO2 absorption. However, there are no accurate computational models that provide accurate predictions of the CO2 solubility in chemically reactive DESs. Here, we develop machine learning (ML) models to predict the solubility of CO2 in chemically reactive DESs. As training data, we collected 214 data points for the CO2 solubility in 149 different chemically reactive DESs at different temperatures, pressures, and DES molar ratios from published work. The physics-driven input features for the ML models include σ-profile descriptors that quantify the relative probability of a molecular surface segment having a certain screening charge density and were calculated with the first-principle quantum chemical method COSMO-RS. We show here that, although COSMO-RS does not explicitly calculate chemical reaction profiles, the COSMO-RS-derived σ-profile features can be used to predict bond formation. Of the models trained, an artificial neural network (ANN) provides the most accurate CO2 solubility prediction with an average absolute relative deviation of 2.94% on the testing sets. Overall, this work provides ML models that can predict CO2 solubility precisely and thus accelerate the design and application of chemically reactive DESs.

Quantum chemistry in the form of density functional theory (DFT) calculations is a powerful numerical experiment for predicting intermolecular interaction energies. However, no chemical insight is gained in this way beyond predictions of observables. Energy decomposition analysis (EDA) can quantitatively bridge this gap by providing values for the chemical drivers of the interactions, such as permanent electrostatics, Pauli repulsion, dispersion, and charge transfer. These energetic contributions are identified by performing DFT calculations with constraints that disable components of the interaction. This review describes the second-generation version of the absolutely localized molecular orbital EDA (ALMO-EDA-II). The effects of different physical contributions on changes in observables such as structure and vibrational frequencies upon complex formation are characterized via the adiabatic EDA. Example applications include red- versus blue-shifting hydrogen bonds; the bonding and frequency shifts of CO, N₂, and BF bound to a [Ru(II)(NH₃)₅]^{2 +} moiety; and the nature of the strongly bound complexes between pyridine and the benzene and naphthalene radical cations. Additionally, the use of ALMO-EDA-II to benchmark and guide the development of advanced force fields for molecular simulation is illustrated with the recent, very promising, MB-UCB potential.

Many proteins have small-molecule binding pockets that are not easily detectable in the ligand-free structures. These cryptic sites require a conformational change to become apparent; a cryptic site can therefore be defined as a site that forms a pocket in a holo structure, but not in the apo structure. Because many proteins appear to lack druggable pockets, understanding and accurately identifying cryptic sites could expand the set of drug targets. Previously, cryptic sites were identified experimentally by fragment-based ligand discovery and computationally by long molecular dynamics simulations and fragment docking. Here, we begin by constructing a set of structurally defined apo-holo pairs with cryptic sites. Next, we comprehensively characterize the cryptic sites in terms of their sequence, structure, and dynamics attributes. We find that cryptic sites tend to be as conserved in evolution as traditional binding pockets but are less hydrophobic and more flexible. Relying on this characterization, we use machine learning to predict cryptic sites with relatively high accuracy (for our benchmark, the true positive and false positive rates are 73% and 29%, respectively). We then predict cryptic sites in the entire structurally characterized human proteome (11,201 structures, covering 23% of all residues in the proteome). CryptoSite increases the size of the potentially "druggable" human proteome from ~40% to ~78% of disease-associated proteins. Finally, to demonstrate the utility of our approach in practice, we experimentally validate a cryptic site in protein tyrosine phosphatase 1B using a covalent ligand and NMR spectroscopy. The CryptoSite Web server is available at http://salilab.org/cryptosite.