Ion solvation thermodynamics from simulation with a polarizable force field

Thermodynamic measurements of the solvation of salts and electrolytes are relatively straightforward, but it is not possible to separate total solvation free energies into distinct cation and anion contributions without reference to an additional extrathermodynamic assumption. The present work attempts to resolve this difficulty using molecular dynamics simulations with the AMOEBA polarizable force field and perturbation techniques to directly compute absolute solvation free energies for potassium, sodium, and chloride ions in liquid water and formamide. Corresponding calculations are also performed with two widely used nonpolarizable force fields. The simulations with the polarizable force field accurately reproduce in vacuo quantum mechanical results, experimental ion-cluster solvation enthalpies, and experimental solvation free energies for whole salts, while the other force fields do not. The results indicate that calculations with a polarizable force field can capture the thermodynamics of ion solvation and that the solvation free energies of the individual ions differ by several kilocalories from commonly cited values.     

Evaluation of solvation free energies for small molecules with the AMOEBA polarizable force field

The effects of electronic polarization in biomolecular interactions will differ depending on the local dielectric constant of the environment, such as in solvent, DNA, proteins, and membranes. Here the performance of the AMOEBA polarizable force field is evaluated under nonaqueous conditions by calculating the solvation free energies of small molecules in four common organic solvents. Results are compared with experimental data and equivalent simulations performed with the GAFF pairwise‐additive force field. Although AMOEBA results give mean errors close to “chemical accuracy,” GAFF performs surprisingly well, with statistically significantly more accurate results than AMOEBA in some solvents. However, for both models, free energies calculated in chloroform show worst agreement to experiment and individual solutes are consistently poor performers, suggesting non‐potential‐specific errors also contribute to inaccuracy. Scope for the improvement of both potentials remains limited by the lack of high quality experimental data across multiple solvents, particularly those of high dielectric constant. © 2016 The Authors. Journal of Computational Chemistry Published by Wiley Periodicals, Inc.     

Effective polarization energy of the naphthalene molecular crystal: a study on the polarizable force field

polarization energy of the localized charge in organic solids consists of electronic polarization energy, permanent electrostatic interactions, and inter/intra molecular relaxation energies. The effective electronic polarization energies for an electron/hole carrier were successfully estimated by AMOEBA polarizable force field in naphthalene molecular crystals. Both electronic polarization energy and permanent electrostatic interaction were in agreement with the preview experimental values. In addition, the influence of the multipoles from different distributed mutipole analysis (DMA) fitting options on the electrostatic interactions are discussed in this paper. We found that the multipoles obtained from Gauss-Hermite quadrature without diffuse function or grid-based quadrature with 0.325 Å H atomic radius will give reasonable electronic polarization energies and permanent interactions for electron and hole carriers.     

Molecular dynamics study on hydrophobic effects in aqueous urea solutions.

Hydrophobic effects in aqueous urea were analyzed by molecular dynamics simulations. The contribution of solvents to the potential of mean force between two methane molecules was calculated by using molecular dynamics simulations and was compared with the solubility data of hydrocarbons in aqueous urea. Both the simulation results and the solubility data indicated that urea stabilizes methane-methane association. The stabilization was due to increasing the solvation free energies of small hydrocarbons such as methane by addition of urea. The solvation free energies of larger hydrocarbons, on the other hand, are decreased by addition of urea. This effect of the solute size on hydrophobic free energies in aqueous urea was also analyzed by using molecular dynamics simulations by means of division of the solvation process into two parts: the cavity formation and the introduction of the solute-solvent attractive interactions. In the cavity formation, urea increased hydrophobic free energies, and in the introduction of the solute-solvent attractive interactions, urea decreased hydrophobic free energies. The influence of urea on hydrophobic free energies was determined by the balance of effects of the two parts of the solvation process.     

Solvation of monovalent anions in acetonitrile and N,N-dimethylformamide: Parameterization of the IEF-PCM model

The present work reports the parameterization of the polarizable continuum model for predicting the free energies of solvation for monovalent anions in acetonitrile and N,N-dimethylformamide. The parameterization of the model for acetonitrile employed the experimental free energies of solvation for a set of 12 charged solutes, containing H, C, N, O, S, F, Cl, Br, and I atoms. For the N,N-dimethylformamide solutions, experimental solvation free energies for 11 monovalent anions were used. A mean absolute error of 0.7 kcal/mol in the solvation free energies has been achieved for the 12 anions in acetonitrile, whereas the mean absolute error for the 11 anions corresponds to 0.5 kcal/mol in N,N-dimethylformamide. These results indicate that the polarizable continuum model is a suitable methodology for the study of thermodynamic effects in solutions of monovalent anions in both solvents.     

Multipole electrostatics in hydration free energy calculations

Hydration free energy (HFE) is generally used for evaluating molecular solubility, which is an important property for pharmaceutical and chemical engineering processes. Accurately predicting HFE is also recognized as one fundamental capability of molecular mechanics force field. Here, we present a systematic investigation on HFE calculations with AMOEBA polarizable force field at various parameterization and simulation conditions. The HFEs of seven small organic molecules have been obtained alchemically using the Bennett Acceptance Ratio method. We have compared two approaches to derive the atomic multipoles from quantum mechanical calculations: one directly from the new distributed multipole analysis and the other involving fitting to the electrostatic potential around the molecules. Wave functions solved at the MP2 level with four basis sets (6-311G*, 6-311++G(2d,2p), cc-pVTZ, and aug-cc-pVTZ) are used to derive the atomic multipoles. HFEs from all four basis sets show a reasonable agreement with experimental data (root mean square error 0.63 kcal/mol for aug-cc-pVTZ). We conclude that aug-cc-pVTZ gives the best performance when used with AMOEBA, and 6-311++G(2d,2p) is comparable but more efficient for larger systems. The results suggest that the inclusion of diffuse basis functions is important for capturing intermolecular interactions. The effect of long-range correction to van der Waals interaction on the hydration free energies is about 0.1 kcal/mol when the cutoff is 12?, and increases linearly with the number of atoms in the solute/ligand. In addition, we also discussed the results from a hybrid approach that combines polarizable solute with fixed-charge water in the HFE calculation.     

Free energy of solvation from molecular dynamics simulation applying Voronoi-Delaunay triangulation to the cavity creation

The free energy of solvation for a large number of representative solutes in various solvents has been calculated from the polarizable continuum model coupled to molecular dynamics computer simulation. A new algorithm based on the Voronoi-Delaunay triangulation of atom-atom contact points between the solute and the solvent molecules is presented for the estimation of the solvent-accessible surface surrounding the solute. The volume of the inscribed cavity is used to rescale the cavitational contribution to the solvation free energy for each atom of the solute atom within scaled particle theory. The computation of the electrostatic free energy of solvation is performed using the Voronoi-Delaunay surface around the solute as the boundary for the polarizable continuum model. Additional short-range contributions to the solvation free energy are included directly from the solute-solvent force field for the van der Waals-type interactions. Calculated solvation free energies for neutral molecules dissolved in benzene, water, CCl4, and octanol are compared with experimental data. We found an excellent correlation between the experimental and computed free energies of solvation for all the solvents. In addition, the employed algorithm for the cavity creation by Voronoi-Delaunay triangulation is compared with the GEPOL algorithm and is shown to predict more accurate free energies of solvation, especially in solvents composed by molecules with nonspherical molecular shapes.     

Toward accurate solvation dynamics of lanthanides and actinides in water using polarizable force fields: from gas-phase energetics to hydration free energies

In this contribution, we focused on the use of polarizable force fields to model the structural, energetic, and thermodynamical properties of lanthanides and actinides in water. In a first part, we chose the particular case of the Th(IV) cation to demonstrate the capabilities of the AMOEBA polarizable force field to reproduce both reference ab initio gas-phase energetics and experimental data including coordination numbers and radial distribution functions. Using such model, we predicted the first polarizable force field estimate of Th(IV) solvation free energy, which accounts for −1,638 kcal/mol. In addition, we proposed in a second part of this work a full extension of the SIBFA (Sum of Interaction Between Fragments Ab initio computed) polarizable potential to lanthanides (La(III) and Lu(III)) and to actinides (Th(IV)) in water. We demonstrate its capabilities to reproduce all ab initio contributions as extracted from energy decomposition analysis computations, including many-body charge transfer and discussed its applicability to extended molecular dynamics and its parametrization on high-level post-Hartree–Fock data.     

Prediction of hydration free energies for the SAMPL4 data set with the AMOEBA polarizable force field

Hydration free energy calculations are often used to validate molecular simulation methodologies and molecular mechanics force fields. We use the free-energy perturbation method together with the AMOEBA polarizable force field and the Poltype parametrization protocol to predict the hydration free energies of 52 molecules as part of the SAMPL4 blind challenge. For comparison, similar calculations are performed using the non-polarizable General Amber force field. Against our expectations, the latter force field gives the better results compared to experiment. One possible explanation is the sensitivity of the AMOEBA results to the conformation used for parametrization.     

Improving the performance of the coupled reference interaction site model-hyper-netted chain (RISM-HNC)/simulation method for free energy of solvation

The coupled reference interaction site model-hyper-netted chain (RISM-HNC)/ simulation methodology determines solvation free energies as a function of the set of all radial distribution functions of solvent atoms about atomic solute sites. These functions are determined from molecular dynamics (MD) or Monte Carlo (MC) simulations rather than from solving the RISM and HNC equations iteratively. Previous applications of the method showed that it can predict relative free energies of solvation for small solutes accurately. However, the errors scale with the system size. In this study, we propose the use of the hard-sphere free energy as the reference and a linear response approximation to improve the performance, i.e., accuracy and robustness, of the method, particularly removing the size dependency of the error. The details of the new formalism are presented. To validate the proposed formalism, solvation free energies of N-methylacetamide and methylamine are computed using the new RISM-HNC-based expressions in addition to a linear response expression, which are compared to previous thermodynamic integration and thermodynamic perturbation results performed with the same force field. Additionally, free energies of solvation for cyclohexane, pyridine, benzene and derivatives, and other small organic molecules are calculated and compared to experimental values.     

A valence bond model for aqueous Cu(II) and Zn(II) ions in the AMOEBA polarizable force field

A general molecular mechanics (MM) model for treating aqueous Cu²⁺ and Zn²⁺ ions was developed based on valence bond (VB) theory and incorporated into the atomic multipole optimized energetics for biomolecular applications (AMOEBA) polarizable force field. Parameters were obtained by fitting MM energies to that computed by ab initio methods for gas‐phase tetra‐ and hexa‐aqua metal complexes. Molecular dynamics (MD) simulations using the proposed AMOEBA‐VB model were performed for each transition metal ion in aqueous solution, and solvent coordination was evaluated. Results show that the AMOEBA‐VB model generates the correct square‐planar geometry for gas‐phase tetra‐aqua Cu²⁺ complex and improves the accuracy of MM model energetics for a number of ligation geometries when compared to quantum mechanical (QM) computations. On the other hand, both AMOEBA and AMOEBA‐VB generate results for Zn²⁺–water complexes in good agreement with QM calculations. Analyses of the MD trajectories revealed a six‐coordination first solvation shell for both Cu²⁺ and Zn²⁺ ions in aqueous solution, with ligation geometries falling in the range reported by previous studies. © 2012 Wiley Periodicals, Inc.     

Molecular dynamics simulations of methane hydrate using polarizable force fields

Molecular dynamics simulations of methane hydrate have been carried out using the polarizable AMOEBA and COS/G2 force fields. Properties calculated include the temperature dependence of the lattice constant, the OC and OO radial distribution functions, and the vibrational spectra. Both the AMOEBA and COS/G2 force fields are found to successfully account for the available experimental data, with overall somewhat better agreement with experiment being found for the AMOEBA model. Comparison is made with previous results obtained using TIP4P and SPC/E effective two-body force fields and the polarizable TIP4P-FQ force field, which allows for in-plane polarization only. Significant differences are found between the properties calculated using the TIP4P-FQ model and those obtained using the other models, indicating an inadequacy of restricting explicit polarization to in-plane only.     

Tinker‐OpenMM: Absolute and relative alchemical free energies using AMOEBA on GPUs

The capabilities of the polarizable force fields for alchemical free energy calculations have been limited by the high computational cost and complexity of the underlying potential energy functions. In this work, we present a GPU‐based general alchemical free energy simulation platform for polarizable potential AMOEBA. Tinker‐OpenMM, the OpenMM implementation of the AMOEBA simulation engine has been modified to enable both absolute and relative alchemical simulations on GPUs, which leads to a ∼200‐fold improvement in simulation speed over a single CPU core. We show that free energy values calculated using this platform agree with the results of Tinker simulations for the hydration of organic compounds and binding of host–guest systems within the statistical errors. In addition to absolute binding, we designed a relative alchemical approach for computing relative binding affinities of ligands to the same host, where a special path was applied to avoid numerical instability due to polarization between the different ligands that bind to the same site. This scheme is general and does not require ligands to have similar scaffolds. We show that relative hydration and binding free energy calculated using this approach match those computed from the absolute free energy approach. © 2017 Wiley Periodicals, Inc.     

Automation of AMOEBA polarizable force field parameterization for small molecules

A protocol to generate parameters for the AMOEBA polarizable force field for small organic molecules has been established, and polarizable atomic typing utility, Poltype, which fully automates this process, has been implemented. For validation, we have compared with quantum mechanical calculations of molecular dipole moments, optimized geometry, electrostatic potential, and conformational energy for a variety of neutral and charged organic molecules, as well as dimer interaction energies of a set of amino acid side chain model compounds. Furthermore, parameters obtained in gas phase are substantiated in liquid-phase simulations. The hydration free energy (HFE) of neutral and charged molecules have been calculated and compared with experimental values. The RMS error for the HFE of neutral molecules is less than 1 kcal/mol. Meanwhile, the relative error in the predicted HFE of salts (cations and anions) is less than 3% with a correlation coefficient of 0.95. Overall, the performance of Poltype is satisfactory and provides a convenient utility for applications such as drug discovery. Further improvement can be achieved by the systematic study of various organic compounds, particularly ionic molecules, and refinement and expansion of the parameter database.     

Calculation of acidic dissociation constants in water: solvation free energy terms. Their accuracy and impact

Three polarizable continuum models, DPCM, CPCM, and IEFPCM, have been applied to calculate free energy differences for nine neutral compounds and their anions. On the basis of solvation free energies, the pK_a values were obtained for the compounds in question by using three thermodynamic cycles: one, involving the combined experimental and calculated data, as well as two other cycles solely with calculated data. This paper deals with the influence of factors such as the SCRF model applied, choice of a particular thermodynamic cycle, atomic radii used to build a cavity in the solvent (water), optimization of geometry in water, inclusion of electron correlation, and the dimension of the basis set on the solvation free energies and on the calculated pK_a values. Electronic supplementary material The online version of this article (doi: ) contains supplementary material, which is available to authorized users.     

Calculation of the free energy of solvation for neutral analogs of amino acid side chains

The ability of the GROMOS96 force field to reproduce partition constants between water and two less polar solvents (cyclohexane and chloroform) for analogs of 18 of the 20 naturally occurring amino acids has been investigated. The estimations of the solvation free energies in water, in cyclohexane solution, and chloroform solution are based on thermodynamic integration free energy calculations using molecular dynamics simulations. The calculations show that while the force field reproduces the experimental solvation free energies of nonpolar analogs with reasonable accuracy the solvation free energies of polar analogs in water are systematically overestimated (too positive). The dependence of the calculated free energies on the atomic partial charges was also studied.     

Polarizable atomic multipole solutes in a Poisson-Boltzmann continuum

Modeling the change in the electrostatics of organic molecules upon moving from vacuum into solvent, due to polarization, has long been an interesting problem. In vacuum, experimental values for the dipole moments and polarizabilities of small, rigid molecules are known to high accuracy; however, it has generally been difficult to determine these quantities for a polar molecule in water. A theoretical approach introduced by Onsager [J. Am. Chem. Soc. 58, 1486 (1936)] used vacuum properties of small molecules, including polarizability, dipole moment, and size, to predict experimentally known permittivities of neat liquids via the Poisson equation. Since this important advance in understanding the condensed phase, a large number of computational methods have been developed to study solutes embedded in a continuum via numerical solutions to the Poisson-Boltzmann equation. Only recently have the classical force fields used for studying biomolecules begun to include explicit polarization in their functional forms. Here the authors describe the theory underlying a newly developed polarizable multipole Poisson-Boltzmann (PMPB) continuum electrostatics model, which builds on the atomic multipole optimized energetics for biomolecular applications (AMOEBA) force field. As an application of the PMPB methodology, results are presented for several small folded proteins studied by molecular dynamics in explicit water as well as embedded in the PMPB continuum. The dipole moment of each protein increased on average by a factor of 1.27 in explicit AMOEBA water and 1.26 in continuum solvent. The essentially identical electrostatic response in both models suggests that PMPB electrostatics offers an efficient alternative to sampling explicit solvent molecules for a variety of interesting applications, including binding energies, conformational analysis, and pK(a) prediction. Introduction of 150 mM salt lowered the electrostatic solvation energy between 2 and 13 kcalmole, depending on the formal charge of the protein, but had only a small influence on dipole moments.     

Combining the polarizable Drude force field with a continuum electrostatic Poisson–Boltzmann implicit solvation model

In this work, we have combined the polarizable force field based on the classical Drude oscillator with a continuum Poisson–Boltzmann/solvent‐accessible surface area (PB/SASA) model. In practice, the positions of the Drude particles experiencing the solvent reaction field arising from the fixed charges and induced polarization of the solute must be optimized in a self‐consistent manner. Here, we parameterized the model to reproduce experimental solvation free energies of a set of small molecules. The model reproduces well‐experimental solvation free energies of 70 molecules, yielding a root mean square difference of 0.8 kcal/mol versus 2.5 kcal/mol for the CHARMM36 additive force field. The polarization work associated with the solute transfer from the gas‐phase to the polar solvent, a term neglected in the framework of additive force fields, was found to make a large contribution to the total solvation free energy, comparable to the polar solute–solvent solvation contribution. The Drude PB/SASA also reproduces well the electronic polarization from the explicit solvent simulations of a small protein, BPTI. Model validation was based on comparisons with the experimental relative binding free energies of 371 single alanine mutations. With the Drude PB/SASA model the root mean square deviation between the predicted and experimental relative binding free energies is 3.35 kcal/mol, lower than 5.11 kcal/mol computed with the CHARMM36 additive force field. Overall, the results indicate that the main limitation of the Drude PB/SASA model is the inability of the SASA term to accurately capture non‐polar solvation effects. © 2018 Wiley Periodicals, Inc.     

Predicting hydrophobic solvation by molecular simulation: 1. Testing united‐atom alkane models

We present a systematic test of the performance of three popular united‐atom force fields—OPLS‐UA, GROMOS and TraPPE—at predicting hydrophobic solvation, more precisely at describing the solvation of alkanes in alkanes. Gibbs free energies of solvation were calculated for 52 solute/solvent pairs from Molecular Dynamics simulations and thermodynamic integration making use of the IBERCIVIS volunteer computing platform. Our results show that all force fields yield good predictions when both solute and solvent are small linear or branched alkanes (up to pentane). However, as the size of the alkanes increases, all models tend to increasingly deviate from experimental data in a systematic fashion. Furthermore, our results confirm that specific interaction parameters for cyclic alkanes in the united‐atom representation are required to account for the additional excluded volume within the ring. Overall, the TraPPE model performs best for all alkanes, but systematically underpredicts the magnitude of solvation free energies by about 6% (RMSD of 1.2 kJ/mol). Conversely, both GROMOS and OPLS‐UA systematically overpredict solvation free energies (by ∼13% and 15%, respectively). The systematic trends suggest that all models can be improved by a slight adjustment of their Lennard‐Jones parameters. © 2016 Wiley Periodicals, Inc.