Development of a Many-Body Force Field for Aqueous Alkali Metal and Halogen Ions: An Energy Decomposition Analysis Guided Approach

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.

What to Make of Zero: Resolving the Statistical Noise from Conformational Reorganization in Alchemical Binding Free Energy Estimates with Metadynamics Sampling

We introduce the self-relative binding free energy (self-RBFE) approach to evaluate the intrinsic statistical variance of dual-topology alchemical binding free energy estimators. The self-RBFE is the relative binding free energy between a ligand and a copy of the same ligand, and its true value is zero. Nevertheless, because the two copies of the ligand move independently, the self-RBFE value produced by a finite-length simulation fluctuates and can be used to measure the variance of the model. The results of this validation provide evidence that a significant fraction of the errors observed in benchmark studies reflect the statistical fluctuations of unconverged estimates rather than the models' accuracy. Furthermore, we find that ligand reorganization is a significant contributing factor to the statistical variance of binding free energy estimates and that metadynamics-accelerated conformational sampling of the torsional degrees of freedom of the ligand can drastically reduce the time to convergence.

Performance and Analysis of the Alchemical Transfer Method for Binding-Free-Energy Predictions of Diverse Ligands

The Alchemical Transfer Method (ATM) is herein validated against the relative binding-free energies (RBFEs) of a diverse set of protein-ligand complexes. We employed a streamlined setup workflow, a bespoke force field, and AToM-OpenMM software to compute the RBFEs of the benchmark set prepared by Schindler and collaborators at Merck KGaA. This benchmark set includes examples of standard small R-group ligand modifications as well as more challenging scenarios, such as large R-group changes, scaffold hopping, formal charge changes, and charge-shifting transformations. The novel coordinate perturbation scheme and a dual-topology approach of ATM address some of the challenges of single-topology alchemical RBFE methods. Specifically, ATM eliminates the need for splitting electrostatic and Lennard-Jones interactions, atom mapping, defining ligand regions, and postcorrections for charge-changing perturbations. Thus, ATM is simpler and more broadly applicable than conventional alchemical methods, especially for scaffold-hopping and charge-changing transformations. Here, we performed well over 500 RBFE calculations for eight protein targets and found that ATM achieves accuracy comparable to that of existing state-of-the-art methods, albeit with larger statistical fluctuations. We discuss insights into the specific strengths and weaknesses of the ATM method that will inform future deployments. This study confirms that ATM can be applied as a production tool for RBFE predictions across a wide range of perturbation types within a unified, open-source framework.

The Energetic Origins of Pi-Pi Contacts in Proteins

Accurate potential energy models of proteins must describe the many different types of noncovalent interactions that contribute to a protein's stability and structure. Pi-pi contacts are ubiquitous structural motifs in all proteins, occurring between aromatic and nonaromatic residues and play a nontrivial role in protein folding and in the formation of biomolecular condensates. Guided by a geometric criterion for isolating pi-pi contacts from classical molecular dynamics simulations of proteins, we use quantum mechanical energy decomposition analysis to determine the molecular interactions that stabilize different pi-pi contact motifs. We find that neutral pi-pi interactions in proteins are dominated by Pauli repulsion and London dispersion rather than repulsive quadrupole electrostatics, which is central to the textbook Hunter-Sanders model. This results in a notable lack of variability in the interaction profiles of neutral pi-pi contacts even with extreme changes in the dielectric medium, explaining the prevalence of pi-stacked arrangements in and between proteins. We also find interactions involving pi-containing anions and cations to be extremely malleable, interacting like neutral pi-pi contacts in polar media and like typical ion-pi interactions in nonpolar environments. Like-charged pairs such as arginine-arginine contacts are particularly sensitive to the polarity of their immediate surroundings and exhibit canonical pi-pi stacking behavior only if the interaction is mediated by environmental effects, such as aqueous solvation.

General Charge Transfer Dipole Model for AMOEBA-Like Force Fields

The development of highly accurate force fields is always an importance aspect in molecular modeling. In this work, we introduce a general damping-based charge transfer dipole (D-CTD) model to describe the charge transfer energy and the corresponding charge flow for H, C, N, O, P, S, F, Cl, and Br elements in common bio-organic systems. Then, two effective schemes to evaluate the charge flow from the corresponding induced dipole moment between the interacting molecules were also proposed and discussed. The potential applicability of the D-CTD model in ion-containing systems was also demonstrated in a series of ion-water complexes including Li⁺, Na⁺, K⁺, Mg²⁺, Ca²⁺, Fe²⁺, Zn²⁺, Pt²⁺, F^-, Cl^-, Br^-, and I^- ions. In general, the D-CTD model demonstrated good accuracy and good transferability in both charge transfer energy and the corresponding charge flow for a wide range of model systems. By distinguishing the intermolecular charge redistribution (charge transfer) under the influence of an external electric field from the accompanying intramolecular charge redistribution (polarization), the D-CTD model is theoretically consistent with current induced dipole-based polarizable dipole models and hence can be easily implemented and parameterized. Along with our previous work in charge penetration-corrected electrostatics, a bottom-up approach constructed water model was also proposed and demonstrated. The structure-maker and structure-breaker roles of cations and anions were also correctly reproduced using Na⁺, K⁺, Cl^-, and I^- ions in the new water model, respectively. This work demonstrates a cost-effective approach to describe the charge transfer phenomena. The water and ion models also show the feasibility of a modulated development approach for future force fields.

Many-Body Effects in Aqueous Systems: Synergies Between Interaction Analysis Techniques and Force Field Development

Interaction analysis techniques, including the many-body expansion (MBE), symmetry-adapted perturbation theory, and energy decomposition analysis, allow for an intuitive understanding of complex molecular interactions. We review these methods by first providing a historical context for the study of many-body interactions and discussing how nonadditivities emerge from Hamiltonians containing strictly pairwise-additive interactions. We then elaborate on the synergy between these interaction analysis techniques and the development of advanced force fields aimed at accurately reproducing the Born-Oppenheimer potential energy surface. In particular, we focus on ab initio-based force fields that aim to explicitly reproduce many-body terms and are fitted to high-level electronic structure results. These force fields generally incorporate many-body effects through (a) parameterization of distributed multipoles, (b) explicit fitting of the MBE, (c) inclusion of many-atom features in a neural network, and (d) coarse-graining of many-body terms into an effective two-body term. We also discuss the emerging use of the MBE to improve the accuracy and speed of ab initio molecular dynamics.

Force Decomposition Analysis: A Method to Decompose Intermolecular Forces into Physically Relevant Component Contributions

Computational quantum chemistry can be more than just numerical experiments when methods are specifically adapted to investigate chemical concepts. One important example is the development of energy decomposition analysis (EDA) to reveal the physical driving forces behind intermolecular interactions. In EDA, typically the interaction energy from a good-quality density functional theory (DFT) calculation is decomposed into multiple additive components that unveil permanent and induced electrostatics, Pauli repulsion, dispersion, and charge-transfer contributions to noncovalent interactions. Herein, we formulate, implement, and investigate decomposing the forces associated with intermolecular interactions into the same components. The resulting force decomposition analysis (FDA) is potentially useful as a complement to the EDA to understand chemistry, while also providing far more information than an EDA for data analysis purposes such as training physics-based force fields. We apply the FDA based on absolutely localized molecular orbitals (ALMOs) to analyze interactions of water with sodium and chloride ions as well as in the water dimer. We also analyze the forces responsible for geometric changes in carbon dioxide upon adsorption onto (and activation by) gold and silver anions. We also investigate how the force components of an EDA-based force field for water clusters, namely MB-UCB, compare to those from force decomposition analysis.

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.

Classical Exchange Polarization: An Anisotropic Variable Polarizability Model

Polarizability, or the tendency of the electron distribution to distort under an electric field, often depends on the local chemical environment. For example, the polarizability of a chloride ion is larger in gas phase compared to a chloride ion solvated in water. This effect is due to the restriction the Pauli exclusion principle places on the allowed electron states. Because no two electrons can occupy the same state, when a highly polarizable atom comes in close contact with other atoms or molecules, the space of allowed states can dramatically decrease. This constraint suggests that an accurate molecular mechanics polarizability model should depend on the radial distance between neighboring atoms. This paper introduces a variable polarizability model within the framework of the HIPPO (Hydrogen-like Intermolecular Polarizable Potential) force field, by damping the polarizability as a function of the orbital overlap of two atoms. This effectively captures the quantum mechanical exchange polarization effects, without explicit utilization of antisymmetrized wave functions. We show that the variable polarizability model remarkably improves the two-body polarization energies and three-body energies of ion-ion and ion-water systems. Under this model, no manual tuning of atomic polarizabilities for monatomic ions is required; the gas-phase polarizability can be used because an appropriate damping function is able to correct the polarizability at short range.