Many biological processes are based on molecular recognition between highly charged molecules such as nucleic acids, inorganic ions, charged amino acids, etc. For such cases, it has been demonstrated that molecular simulations with fixed partial charges often fail to achieve experimental accuracy. Although incorporation of more advanced electrostatic models (such as multipoles, mutual polarization, etc.) can significantly improve simulation accuracy, it increases computational expense by a factor of 5–20×. Indirect free energy (IFE) methods can mitigate this cost by modeling intermediate states at fixed-charge resolution. For example, an efficient “reference” model such as a pairwise Amber, CHARMM, or OPLS-AA force field can be used to derive an initial estimate, followed by thermodynamic corrections to a more advanced “target” potential such as the polarizable AMOEBA model. Unfortunately, all currently described IFE methods encounter difficulties reweighting more than ∼50 atoms between resolutions due to extensive scaling of both the magnitude of the thermodynamic corrections and their statistical uncertainty. We present an approach called “simultaneous bookending” (SB) that is fundamentally different from existing IFE methods based on a tunable sampling approximation, which permits scaling to thousands of atoms. SB is demonstrated on the relative binding affinity of Mg2+/Ca2+ to a set of metalloproteins with up to 2972 atoms, finding no statistically significant difference between direct AMOEBA results and those from correcting Amber to AMOEBA. The ability to change the resolution of thousands of atoms during reweighting suggests the approach may be applicable in the future to protein–protein binding affinities or nucleic acid thermodynamics.

The transformation of a pharmaceutical solid from an anhydrous crystal into a hydrated form during drug development represents a risk to a product’s safety and efficacy due to the potential impact on stability, bioavailability, and manufacturability. In this work, we examine 10 classical free energy simulation protocols to evaluate the thermodynamic stability of hydrated crystals relative to their anhydrous forms. Molecular dynamics simulations are used to compute the Gibbs free energies of the crystals of three pharmaceutically relevant systems using two fixed-charge potentials, GAFF and OPLS, as well as the polarizable AMOEBA model. In addition, we explore a variety of water models, including TIP3P, TIP4P, and AMOEBA, for both the interstitial water and the effects of ambient humidity. The AMOEBA model predicts free energy values most consistent with experimental measurements among the models examined. The benefits of a fully polarizable water model relative to fixed-charged models appear to derive predominantly from a better treatment of water’s dipole moment in the crystalline phase. Despite this improved physical treatment, we find that no single model produces reliable predictions of the phase boundary between hydrated and anhydrous crystals from theory alone. However, we show that accurate phase diagrams can be constructed from the simulations by introducing a single experimentally determined coexistence point. With this single experimental data point as input, the phase boundary is correctly predicted within 10% relative humidity on the temperature range of 15 to 75 °C for all three systems examined. Furthermore, we demonstrate that with this known coexistence point as an input, the differences between the various potentials and the water models become insignificant, as all models yield accurate phase boundaries regardless of whether polarization is included due to significant temperature-dependent error cancellation between models.

Computational protein design, ab initio protein/RNA folding, and protein-ligand screening can be too computationally demanding for explicit treatment of solvent. For these applications, implicit solvent offers a compelling alternative, which we describe here for the polarizable atomic multipole AMOEBA force field based on three treatments of continuum electrostatics: numerical solutions to the Poisson-Boltzmann equation (PBE), the domain-decomposition Conductor-like Screening Model (ddCOSMO) approximation to the PBE, and the analytic generalized Kirkwood (GK) approximation. The continuum electrostatic models are combined with a nonpolar estimator based on novel cavitation and dispersion terms. Electrostatic model parameters are numerically optimized using a least squares style target function based on a library of 103 small molecule solvation free energy differences. Mean signed errors for the APBS, ddCOSMO, and GK models are 0.05, 0.00, and 0.00 kcal/mol, respectively, while the mean unsigned errors are 0.70, 0.63, and 0.51 kcal/mol, respectively. Validation of the electrostatic response of the resulting implicit solvents, which are available in the Tinker (or Tinker-HP), OpenMM, and Force Field X software packages, is based on comparisons to explicit solvent simulations for a series of proteins and nucleic acids. Overall, the emergence of performative implicit solvent models for polarizable force fields will open the door to their use for folding and design applications.

Computational protein design, ab initio protein/RNA folding, and protein-ligand screening can be too computationally demanding for explicit treatment of solvent. For these applications, implicit solvent offers a compelling alternative, which we describe here for the polarizable atomic multipole AMOEBA force field based on three treatments of continuum electrostatics: numerical solutions to the Poisson-Boltzmann equation (PBE), the domain-decomposition Conductor-like Screening Model (ddCOSMO) approximation to the PBE, and the analytic generalized Kirkwood (GK) approximation. The continuum electrostatic models are combined with a nonpolar estimator based on novel cavitation and dispersion terms. Electrostatic model parameters are numerically optimized using a least squares style target function based on a library of 103 small molecule solvation free energy differences. Mean signed errors for the APBS, ddCOSMO, and GK models are 0.05, 0.00, and 0.00 kcal/mol, respectively, while the mean unsigned errors are 0.70, 0.63, and 0.51 kcal/mol, respectively. Validation of the electrostatic response of the resulting implicit solvents, which are available in the Tinker (or Tinker-HP), OpenMM, and Force Field X software packages, is based on comparisons to explicit solvent simulations for a series of proteins and nucleic acids. Overall, the emergence of performative implicit solvent models for polarizable force fields will open the door to their use for folding and design applications.<br>