United polarizable multipole water model for molecular mechanics simulation

Computational exploration of naturally derived peptides inhibitory mechanisms against ACE enzyme, from interactions to structural-dynamics

Naturally derived peptides have gained significant attention because of their potential to reduce blood pressure. These peptides can be derived from various sources, including snake venom, marine organisms, cow's milk, seahorses, and plants. Investigating the underlying mechanisms of these peptides in lowering blood pressure is crucial for replacing synthetic drugs currently in use. In this regard, we conducted in silico studies, such as molecular docking, classical molecular dynamics (MD) simulations, and QM/MM simulation methods, on two naturally derived peptides (PAGPRGPA and WALKGYK) and Angiotensin II (ANGII) as the enzyme substrates for binding to sACE. The results of 500 ns MD calculations showed that the WALKGYK peptide occupies the ACE binding site, similar to the ANGII and BPPb peptides, two distinguished C-domain-specific inhibitors. Furthermore, QM/MM calculations demonstrated that no peptide bond cleavage was mediated by Zn2+ at the catalytic site of the peptides. However, during the 500 ns MD simulation, the backbone oxygens of LYS4 and GLY5 of the WALKGYK peptides were tightly coordinated to the Zn2+ ion. Free energy calculations also confirmed the higher affinity of the WALKGYK peptides for binding to sACE. In addition, structural analysis correlation showed a different pattern in PAGPRGPA compared to WALKGYK and ANGII. Despite the similarity of the peptide PAGPRGPA to typical ACE peptide inhibitors with hydrophobic ends, the electrostatic composition of the WALKGYK peptide showed a higher tendency towards ACE inhibition. Therefore, peptide residue compositions, such as WALKGYK, can be considered for designing new inhibitors with fewer side effects for sACE C-terminal inhibitors.

Influence of the Lennard-Jones Combination Rules on the Simulated Properties of Organic Liquids at Optimal Force-Field Parametrization

We recently introduced the CombiFF scheme [Oliveira et al., J. Chem. Theory Comput. 2020, 16, 7525], an approach for the automated refinement of force-field parameters against experimental condensed-phase data for large compound families. Using this scheme, once the time-consuming task of target-data selection and curation has been performed, the force-field optimization itself is both straightforward and fast. As a result, CombiFF provides an ideal framework for evaluating the influence of functional-form decisions on the accuracy of a force field at an optimal level of parametrization. We already used this approach to assess the effect of using an all-atom representation compared to united-atom representations in the force field [Oliveira et al., J. Chem. Theory Comput. 2022, 18, 6757]. Here, CombiFF is applied to assess the effect of three Lennard-Jones combination rules, geometric-mean (GM), Lorentz-Berthelot (LB), or Waldman-Hagler (WH), on the simulated properties of organic liquids. The comparison is performed in terms of the experimental liquid density ρ_liq, vaporization enthalpy ΔH_vap, surface-tension coefficient γ, static relative dielectric permittivity ϵ, and self-diffusion coefficient D. The calibrations of the three force-field variants are carried out independently against 2044 experimental values for ρ_liq, and ΔH_vap concerning 1516 compounds. The resulting root-mean-square deviations from experiment are 30.0, 26.9, and 36.7 kg m^-3 for ρ_liq and 2.8, 2.8, and 2.9 kJ mol^-1 for ΔH_vap, when applying the GM, LB, and WH combination rules, respectively. In terms of these (and the other) properties, the three combination rules perform comparatively well, with the GM and LB results being more similar to each other and slightly more accurate compared to experiment. In contrast, the use of distinct combination rules for the parameter calibration and property calculation leads to much larger errors.

Comparison of the United- and All-Atom Representations of (Halo)alkanes Based on Two Condensed-Phase Force Fields Optimized against the Same Experimental Data Set

The level of accuracy that can be achieved by a force field is influenced by choices made in the interaction-function representation and in the relevant simulation parameters. These choices, referred to here as functional-form variants (FFVs), include for example the model resolution, the charge-derivation procedure, the van der Waals combination rules, the cutoff distance, and the treatment of the long-range interactions. Ideally, assessing the effect of a given FFV on the intrinsic accuracy of the force-field representation requires that only the specific FFV is changed and that this change is performed at an optimal level of parametrization, a requirement that may prove extremely challenging to achieve in practice. Here, we present a first attempt at such a comparison for one specific FFV, namely the choice of a united-atom (UA) versus an all-atom (AA) resolution in a force field for saturated acyclic (halo)alkanes. Two force-field versions (UA vs AA) are optimized in an automated way using the CombiFF approach against 961 experimental values for the pure-liquid densities ρ_liq and vaporization enthalpies ΔH_vap of 591 compounds. For the AA force field, the torsional and third-neighbor Lennard-Jones parameters are also refined based on quantum-mechanical rotational-energy profiles. The comparison between the UA and AA resolutions is also extended to properties that have not been included as parameterization targets, namely the surface-tension coefficient γ, the isothermal compressibility κ_T, the isobaric thermal-expansion coefficient α_P, the isobaric heat capacity c_P, the static relative dielectric permittivity ϵ, the self-diffusion coefficient D, the shear viscosity η, the hydration free energy ΔG_wat, and the free energy of solvation ΔG_che in cyclohexane. For the target properties ρ_liq and ΔH_vap, the UA and AA resolutions reach very similar levels of accuracy after optimization. For the nine other properties, the AA representation leads to more accurate results in terms of η; comparably accurate results in terms of γ, κ_T, α_P, ϵ, D, and ΔG_che; and less accurate results in terms of c_P and ΔG_wat. This work also represents a first step toward the calibration of a GROMOS-compatible force field at the AA resolution.

Automation of AMOEBA polarizable force field for small molecules: Poltype 2

A next-generation protocol (Poltype 2) has been developed which automatically generates AMOEBA polarizable force field parameters for small molecules. Both features and computational efficiency have been drastically improved. Notable advances include improved database transferability using SMILES, robust torsion fitting, non-aromatic ring torsion parameterization, coupled torsion-torsion parameterization, Van der Waals parameter refinement using ab initio dimer data and an intelligent fragmentation scheme that produces parameters with dramatically reduced ab initio computational cost. Additional improvements include better local frame assignment for atomic multipoles, automated formal charge assignment, Zwitterion detection, smart memory resource defaults, parallelized fragment job submission, incorporation of Psi4 quantum package, ab initio error handling, ionization state enumeration, hydration free energy prediction and binding free energy prediction. For validation, we have applied Poltype 2 to ~1000 FDA approved drug molecules from DrugBank. The ab initio molecular dipole moments and electrostatic potential values were compared with Poltype 2 derived AMOEBA counterparts. Parameters were further substantiated by calculating hydration free energy (HFE) on 40 small organic molecules and were compared with experimental data, resulting in an RMSE error of 0.59 kcal/mol. The torsion database has expanded to include 3543 fragments derived from FDA approved drugs. Poltype 2 provides a convenient utility for applications including binding free energy prediction for computational drug discovery. Further improvement will focus on automated parameter refinement by experimental liquid properties, expansion of the Van der Waals parameter database and automated parametrization of modified bio-fragments such as amino and nucleic acids.

O-Acetyl Migration within the Sialic Acid Side Chain: A Mechanistic Study Using the Ab Initio Nanoreactor

Many disease-causing viruses target sialic acids on the surface of host cells. Some viruses bind preferentially to sialic acids with O-acetyl modification at the hydroxyl group of C7, C8, or C9 on the glycerol-like side chain. Studies of proteins binding to sialosides containing O-acetylated sialic acids are crucial in understanding the related diseases but experimentally difficult due to the lability of the ester group. We recently showed that O-acetyl migration among hydroxyl groups of C7, C8, and C9 in sialic acids occurs in all directions in a pH-dependent manner. In the current study, we elucidate a full mechanistic pathway for the migration of O-acetyl among C7, C8, and C9. We used an ab initio nanoreactor to explore potential reaction pathways and density functional theory, pK_a calculations, and umbrella sampling to investigate elementary steps of interest. We found that when a base is present, migration is easy in any direction and involves three key steps: deprotonation of the hydroxyl group, cyclization between the two carbons, and the migration of the O-acetyl group. This dynamic equilibrium may play a defensive role against pathogens that evolve to gain entry to the cell by binding selectively to one acetylation state.

Comparison of molecular dynamics simulations of water with neutron and X-ray scattering experiments

The atomistic structure and dynamics obtained from molecular dynamics (MD) simulations with the example of TIP3P (rigid and flexible) and TIP4P/2005 (rigid) water is compared to neutron and X-ray scattering data at ambient conditions. Neutron and X-ray diffractograms are calculated from the simulations for four isotopic substitutions as well as the incoherent intermediate scattering function for neutrons. The resulting curves are compared to each other and to published experimental data. Differences between simulated and measured intermediate scattering functions are quantified by fitting an analytic model to the computed values. The sensitivity of the scattering curves to the parameters of the MD simulations is demonstrated on the example of two parameters, bond length and angle.

The automated optimisation of a coarse-grained force field using free energy data

Atomistic models provide a detailed representation of molecular systems, but are sometimes inadequate for simulations of large systems over long timescales. Coarse-grained models enable accelerated simulations by reducing the number of degrees of freedom, at the cost of reduced accuracy. New optimisation processes to parameterise these models could improve their quality and range of applicability. We present an automated approach for the optimisation of coarse-grained force fields, by reproducing free energy data derived from atomistic molecular simulations. To illustrate the approach, we implemented hydration free energy gradients as a new target for force field optimisation in ForceBalance and applied it successfully to optimise the un-charged side-chains and the protein backbone in the SIRAH protein coarse-grain force field. The optimised parameters closely reproduced hydration free energies of atomistic models and gave improved agreement with experiment.

Systematic optimization of a fragment-based force field against experimental pure-liquid properties considering large compound families: application to oxygen and nitrogen compounds

The CombiFF approach is a workflow for the automated refinement of force-field parameters against experimental condensed-phase data, considering entire classes of organic molecules constructed using a fragment library via combinatorial isomer enumeration. One peculiarity of this approach is that it relies on an electronegativity-equalization scheme to account for induction effects within molecules, with values of the atomic hardness and electronegativity as electrostatic parameters, rather than the partial charges themselves. In a previous article [M. P. Oliveira, M. Andrey, S. R. Rieder, L. Kern, D. F. Hahn, S. Riniker, B. A. C. Horta and P. H. Hünenberger, J. Chem. Theory. Comput. 2020, 16, 7525], CombiFF was introduced and applied to calibrate a GROMOS-compatible united-atom force field for the saturated acyclic (halo-)alkane family. Here, this scheme is employed for the construction of a corresponding force field for saturated acyclic compounds encompassing eight common chemical functional groups involving oxygen and/or nitrogen atoms, namely: ether, aldehyde, ketone, ester, alcohol, carboxylic acid, amine, and amide. Monofunctional as well as homo-polyfunctional compounds are considered. A total of 1712 experimental liquid densities ρliq and vaporization enthalpies ΔHvap concerning 1175 molecules are used for the calibration (339 molecules) and validation (836 molecules) of the 102 non-bonded interaction parameters of the force field. Using initial parameter values based on the GROMOS 2016H66 parameter set, convergence is reached after five iterations. Given access to one processor per simulated system, this operation only requires a few days of wall-clock computing time. After optimization, the root-mean-square deviations from experiment are 29.9 (22.4) kg m-3 for ρliq and 4.1 (5.5) kJ mol-1 for ΔHvap for the calibration (validation) set. Thus, a very good level of agreement with experiment is achieved in terms of these two properties, although the errors are inhomogeneously distributed across the different chemical functional groups.

Solvent-scaling as an alternative to coarse-graining in adaptive-resolution simulations: The adaptive solvent-scaling (AdSoS) scheme

A new approach termed Adaptive Solvent-Scaling (AdSoS) is introduced for performing simulations of a solute embedded in a fine-grained (FG) solvent region itself surrounded by a coarse-grained (CG) solvent region, with a continuous FG ↔ CG switching of the solvent resolution across a buffer layer. Instead of relying on a distinct CG solvent model, the AdSoS scheme is based on CG models defined by a dimensional scaling of the FG solvent by a factor s, accompanied by an s-dependent modulation of the atomic masses and interaction parameters. The latter changes are designed to achieve an isomorphism between the dynamics of the FG and CG models, and to preserve the dispersive and dielectric solvation properties of the solvent with respect to a solute at FG resolution. This scaling approach offers a number of advantages compared to traditional coarse-graining: (i) the CG parameters are immediately related to those of the FG model (no need to parameterize a distinct CG model); (ii) nearly ideal mixing is expected for CG variants with similar s-values (ideal mixing holding in the limit of identical s-values); (iii) the solvent relaxation timescales should be preserved (no dynamical acceleration typical for coarse-graining); (iv) the graining level N_G (number of FG molecules represented by one CG molecule) can be chosen arbitrarily (in particular, N_G = s³ is not necessarily an integer); and (v) in an adaptive-resolution scheme, this level can be varied continuously as a function of the position (without requiring a bundling mechanism), and this variation occurs at a constant number of particles per molecule (no occurrence of fractional degrees of freedom in the buffer layer). By construction, the AdSoS scheme minimizes the thermodynamic mismatch between the different regions of the adaptive-resolution system, leading to a nearly homogeneous scaled solvent density s³ρ. Residual density artifacts in and at the surface of the boundary layer can easily be corrected by means of a grid-based biasing potential constructed in a preliminary pure-solvent simulation. This article introduces the AdSoS scheme and provides an initial application to pure atomic liquids (no solute) with Lennard-Jones plus Coulomb interactions in a slab geometry.

Abiraterone and D4, 3-keto Abiraterone binding to CYP17A1, a structural comparison study by molecular dynamic simulation

The importance of computer-aided drug design and development is clear nowadays. These approaches smooth the way of designing some efficient candidates based on drugs in use. At this place, we studied the mechanism of D4-abiraterone (D4A), the active metabolite of Abiraterone (Abi), binding to CYP17A1 compared with Abi. The molecular dynamics simulation results reveal that the metabolite, which lacks the key 3β-OH group, has a varied H-bond forming pattern. The critical H-bond between 3β-OH of Abi with Asn_202 turns to 3 Keto-O of D4A with Arg_239 in the substrate-binding site. This interaction causes a remarkable distance of 0.63 nm between D4A nitrogen and Fe in heme, which reduces its 17,20 lyase selectivity. The D4A keto moiety presents an immense number of H-bond with surrounding solvent molecules compared with the Abi hydroxyl group. As a result, D4A develops a weaker H-bond network with the enzyme. Otherwise, the heterocyclic nature of inhibitors helps for noticeable van der Waals interaction formation with CYP17A1. However, Abi stabilized position in the binding site helps more van der Waals interactions deposition than D4A. These results convinced the importance of the conserved H-bond for acquiring the proper position by the substrate or inhibitor in the binding site.

Water in Nanopores and Biological Channels: A Molecular Simulation Perspective

This Review explores the dynamic behavior of water within nanopores and biological channels in lipid bilayer membranes. We focus on molecular simulation studies, alongside selected structural and other experimental investigations. Structures of biological nanopores and channels are reviewed, emphasizing those high-resolution crystal structures, which reveal water molecules within the transmembrane pores, which can be used to aid the interpretation of simulation studies. Different levels of molecular simulations of water within nanopores are described, with a focus on molecular dynamics (MD). In particular, models of water for MD simulations are discussed in detail to provide an evaluation of their use in simulations of water in nanopores. Simulation studies of the behavior of water in idealized models of nanopores have revealed aspects of the organization and dynamics of nanoconfined water, including wetting/dewetting in narrow hydrophobic nanopores. A survey of simulation studies in a range of nonbiological nanopores is presented, including carbon nanotubes, synthetic nanopores, model peptide nanopores, track-etched nanopores in polymer membranes, and hydroxylated and functionalized nanoporous silica. These reveal a complex relationship between pore size/geometry, the nature of the pore lining, and rates of water transport. Wider nanopores with hydrophobic linings favor water flow whereas narrower hydrophobic pores may show dewetting. Simulation studies over the past decade of the behavior of water in a range of biological nanopores are described, including porins and β-barrel protein nanopores, aquaporins and related polar solute pores, and a number of different classes of ion channels. Water is shown to play a key role in proton transport in biological channels and in hydrophobic gating of ion channels. An overall picture emerges, whereby the behavior of water in a nanopore may be predicted as a function of its hydrophobicity and radius. This informs our understanding of the functions of diverse channel structures and will aid the design of novel nanopores. Thus, our current level of understanding allows for the design of a nanopore which promotes wetting over dewetting or vice versa. However, to design a novel nanopore, which enables fast, selective, and gated flow of water de novo would remain challenging, suggesting a need for further detailed simulations alongside experimental evaluation of more complex nanopore systems.

Three-site and five-site fixed-charge water models compatible with AMOEBA force field

In a typical biomolecular simulation using Atomic Multipole Optimized Energetics for Biomolecular Applications (AMOEBA) force field, the vast majority molecules in the simulation box consist of water, and these water molecules consume the most CPU power due to the explicit mutual induction effect. To improve the computational efficiency, we here develop two new nonpolarizable water models (with flexible bonds and fixed charges) that are compatible with AMOEBA solute: the 3-site AW3C and 5-site AW5C. To derive the force-field parameters for AW3C and AW5C, we fit to six experimental liquid thermodynamic properties: liquid density, enthalpy of vaporization, dielectric constant, isobaric heat capacity, isothermal compressibility and thermal expansion coefficient, at a broad range of temperatures from 261.15 to 353.15 K under 1.0 atm pressure. We further validate our AW3C and AW5C water models by showing that they can well reproduce the radial distribution function g(r), self-diffusion constant D, and hydration free energy from the AMOEBA03 water model and the experimental observations. Furthermore, we show that our AW3C and AW5C water models can greatly accelerate (>5 times) the bulk water as well as biomolecular simulations when compared to AMOEBA water. Specifically, we demonstrate that the applications of AW3C and AW5C water models to simulate a DNA duplex lead to a threefold acceleration, and in the meanwhile well maintain the structural properties as the fully polarizable AMOEBA water. We expect that our AW3C and AW5C water models hold great promise to be widely applied to simulate complex bio-molecules using the AMOEBA force field.

Support Vector Regression-Based Monte Carlo Simulation of Flexible Water Clusters

Molecular simulations based on classical force fields are computationally efficient but lack accuracy due to the empirical formulation of non-bonded interactions. Quantum mechanical (QM) methods, albeit accurate, have inhibitory computational costs for large molecules and clusters. Hence, to overcome the bottleneck, machine learning (ML)-based methods have been employed in the recent years. We had earlier reported a combined scheme of many-body expansion (MBE) and ML to predict the interaction energies of rigid water clusters. In this work, we proceed toward building a flexible water model using the ML-MBE scheme. This ML-MBE scheme has an error of <1% for interaction energy prediction in comparison to the parent QM method for flexible water decamers. Machine learning-based Monte Carlo simulations (MLMC) are performed with this water model, and the structural properties of these configurations are compared with those obtained from ab initio molecular dynamics (AIMD) and the TIP3P classical force field. The radial distribution functions, tetrahedral order parameters, and number of hydrogen bonds in AIMD and MLMC have a similar qualitative and quantitative trend, whereas the classical force fields show a significant deviation.

Bottom-up derived flexible water model with dipole and quadrupole moments for coarse-grained molecular simulations

A bottom-up CG water model is developed to capture the electrostatic multipoles, structural correlation and thermodynamics of water.

Torsional Energy Barriers of Biaryls Could Be Predicted by Electron Richness/Deficiency of Aromatic Rings; Advancement of Molecular Mechanics toward Atom-Type Independence

Biaryl molecules are ubiquitous pharmacophores found in natural products and pharmaceuticals. In spite of this, existing molecular mechanics force fields are unable to accurately reproduce their torsional energy profiles, except for a few well-parametrized cases. This effectively limits the ability of structure-based drug design methods to correctly identify hits involving biaryls with confidence (e.g., during virtual screening, employing docking and/or molecular dynamics simulations). Continuing in our endeavor to quantify organic chemistry principles, we showed that the torsional energy profile of biaryl compounds could be computed on-the-fly based on the electron richness/deficiency of the aromatic rings. This method, called H-TEQ 4.0, was developed using a set of 131 biaryls. It was subsequently validated on a separate set of 100 diverse biaryls, including multisubstituted, bicyclic and tricyclic druglike molecules, and produced an average root-mean-square error (RMSE) of 0.95 kcal·mol^-1. For comparison, GAFF2 produced an RMSE of 3.88 kcal·mol^-1, owing to problems associated with the transferability of torsion parameters. The success of H-TEQ 4.0 provided further evidence that force fields could transition to become atom-type independent, providing that the correct chemical principles are used. Overall, this method solved the problem of transferability of biaryl torsion parameters, while simultaneously improving the overall accuracy of the force field.

AMOEBA+ Classical Potential for Modeling Molecular Interactions

Classical potentials based on isotropic and additive atomic charges have been widely used to model molecules in computers for the past few decades. The crude approximations in the underlying physics are hindering both their accuracy and transferability across chemical and physical environments. Here we present a new classical potential, AMOEBA+, to capture essential intermolecular forces, including permanent electrostatics, repulsion, dispersion, many-body polarization, short-range charge penetration, and charge transfer, by extending the polarizable multipole-based AMOEBA (Atomic Multipole Optimized Energetics for Biomolecular Applications) model. For a set of common organic molecules, we show that AMOEBA+ with general parameters can reproduce both quantum mechanical interactions and energy decompositions according to Symmetry-Adapted Perturbation Theory (SAPT). Additionally, a new water model based on the AMOEBA+ framework captures various liquid-phase properties in molecular dynamics simulations while remaining consistent with SAPT energy decompositions, utilizing both ab initio data and experimental liquid properties. Our results demonstrate that it is possible to improve the physical basis of classical force fields to advance their accuracy and general applicability.

Polarizable Force Fields for Biomolecular Simulations: Recent Advances and Applications

Realistic modeling of biomolecular systems requires an accurate treatment of electrostatics, including electronic polarization. Due to recent advances in physical models, simulation algorithms, and computing hardware, biomolecular simulations with advanced force fields at biologically relevant timescales are becoming increasingly promising. These advancements have not only led to new biophysical insights but also afforded opportunities to advance our understanding of fundamental intermolecular forces. This article describes the recent advances and applications, as well as future directions, of polarizable force fields in biomolecular simulations.

A QM/MM Derived Polarizable Water Model for Molecular Simulation

In this work, we propose an improved QM/MM-based strategy to determine condensed-phase polarizabilities and we use this approach to optimize a new and simple polarizable four-site water model for classical molecular simulation. For the determination of the model value for the polarizability from QM/MM, we show that our proposed consensus-fitting strategy significantly reduces the uncertainty in calculated polarizabilities in cases where the size of the local external electric field is small. By fitting electrostatic, polarization and dispersion properties of our water model based on quantum and/or combined QM/MM calculations, only a single model parameter (describing exchange repulsion) is left for empirical calibration. The resulting model performs well in describing relevant pure-liquid thermodynamic and transport properties, which illustrates the merit of our approach to minimize the number of free variables in our model.

Assimilating Radial Distribution Functions To Build Water Models with Improved Structural Properties

The structural properties of three- and four-site water models are improved by extending the ForceBalance parametrization code to include a new methodology allowing for the targeting of any radial distribution function (RDF) during the parametrization of a force field. The mean squared difference (MSD) between the experimental and simulated RDFs contributes to an objective function, allowing for the systematic optimization of force field parameters to reach closer overall agreement with experiment. RDF fitting is applied to develop modified versions of the TIP3P and TIP4P/2005 water models in which the Lennard-Jones potential is replaced by a Buckingham potential. The optimized TIP3P-Buckingham and TIP4P-Buckingham potentials feature 93 and 98% lower MSDs in the OO RDF compared to the TIP3P and TIP4P/2005 models respectively, with marked decreases in the height of the first peak. Additionally, these Buckingham models predict the entropy of water more accurately, reducing the error in the entropy of TIP3P from 11 to 3% and the error in the entropy of TIP4P/2005 from 11 to 2%. These new Buckingham models have improved predictive power for many nonfitted properties particularly in the case of TIP3P. Our work directly demonstrates how the Buckingham potential can improve the description of water's structural properties beyond the Lennard-Jones potential. Moreover, adding a Buckingham potential is a favorable alternative to adding interaction sites in terms of computational speed on modern GPU hardware.

An empirical extrapolation scheme for efficient treatment of induced dipoles

Many cutting edge force fields include polarization, to enhance their accuracy and range of applicability. In this work, we develop efficient strategies for the induced dipole polarization method. By fitting various orders of perturbation theory (PT) dipoles to a diverse training set, we arrive at a family of fully analytic methods - whose nth order is referred to OPTn - that span the full spectrum of polarization methods from the fast zeroth-order approach that neglects mutual dipole coupling, approaching the fully variational approach at high order. Our training set contains many difficult cases where the PT series diverges, and we demonstrate that our OPTn methods still deliver excellent results in these cases. Our tests show that the OPTn methods exhibit rapid convergence towards the exact answer with each increasing PT order. The fourth order OPT4 method, whose costs are commensurate with three iterations of the leading conjugate gradient method, is a particularly promising candidate to be used as a drop-in replacement for existing solvers without further parameterization.

Evaluating Parametrization Protocols for Hydration Free Energy Calculations with the AMOEBA Polarizable Force Field

Hydration free energy (HFE) calculations are often used to assess the performance of biomolecular force fields and the quality of assigned parameters. The AMOEBA polarizable force field moves beyond traditional pairwise additive models of electrostatics and may be expected to improve upon predictions of thermodynamic quantities such as HFEs over and above fixed-point-charge models. The recent SAMPL4 challenge evaluated the AMOEBA polarizable force field in this regard but showed substantially worse results than those using the fixed-point-charge GAFF model. Starting with a set of automatically generated AMOEBA parameters for the SAMPL4 data set, we evaluate the cumulative effects of a series of incremental improvements in parametrization protocol, including both solute and solvent model changes. Ultimately, the optimized AMOEBA parameters give a set of results that are not statistically significantly different from those of GAFF in terms of signed and unsigned error metrics. This allows us to propose a number of guidelines for new molecule parameter derivation with AMOEBA, which we expect to have benefits for a range of biomolecular simulation applications such as protein-ligand binding studies.