Determination of partial molar volumes from free energy perturbation theory

Enthalpies and entropies of hydration from Monte Carlo simulations

The changes in free energy, enthalpy, and entropy for transfer of a solute from the gas phase into solution are the fundamental thermodynamic quantities that characterize the solvation process. Owing to the development of methods based on free-energy perturbation theory, computation of free energies of solvation has become routine in conjunction with Monte Carlo (MC) statistical mechanics and molecular dynamics (MD) simulations. Computation of the enthalpy change and by inference the entropy change is more challenging. Two methods are considered in this work corresponding to direct averaging for the solvent and solution and to computing the temperature derivative of the free energy in the van't Hoff approach. The application is for neutral organic solutes in TIP4P water using long MC simulations to improve precision. Definitive results are also provided for pure TIP4P water. While the uncertainty in computed free energies of hydration is ca. 0.05 kcal mol-1, it is ca. 0.4 kcal mol-1 for the enthalpy changes from either van't Hoff plots or the direct method with sampling for 5 billion MC configurations. Partial molar volumes of hydration are also computed by the direct method; they agree well with experimental data with an average deviation of 3 cm3 mol-1. In addition, the results permit breakdown of the errors in the free energy changes from the OPLS-AA force field into their enthalpic and entropic components. The excess hydrophobicity of organic solutes is enthalpic in origin.

Hydrophilic Versus Hydrophobic Coupling in the Pressure Dependence of the Chemical Potential of Alkali Metal and Halide Ions in Water

We computed the chemical potential for some alkali metal ions (K⁺, Rb⁺, and Cs⁺) and two halide ions (Br^- and I^-) in aqueous solution at ambient T and various pressures in the range 1-8000 atm. Results were obtained from classic Monte Carlo simulations in the NPT ensemble by means of the free energy perturbation method. Here, the chemical potential is computed as the sum of a term relative to a Lennard-Jones solute and a term relative to the process in which this solute is transformed into the ion. Hydrophobic and hydrophilic features of these two components of the chemical potential show opposite behaviors under isothermal compression. The increase in pressure determines an increase in the hydrophobic component, which becomes more positive with a stronger effect for larger ions. Correspondingly, the values of the hydrophilic component become more negative for alkali ions, whereas they are only slightly affected by compression for halide ions. Hydrophobic-hydrophilic quasi-compensation in the slopes is observed for Rb⁺. For a smaller ion, such as K⁺, the dependence on pressure of the hydrophilic component is slightly dominant. For a larger ion, as observed in the cases of Cs⁺, Br^-, and I^-, the hydrophobic component assumes the determinant role. Pressure dependence of the chemical potential is little affected by corrections introduced for molecular potential truncation. This view can change for possible boundary artifacts that could have affected the static electrostatic potential. Some inference is obtained from comparison with experimental data at 1 atm on the free energy of hydration. Discrepancies show the characteristic asymmetry between cations and anions. The further addition of a correction based on the static potential significantly reduces these discrepancies with important error cancellation on the sum of chemical potentials of ions of opposite charge. The correction is applied also at higher pressures, and results are compared with those obtained by adding an alternative correction that is based on the water number density. Regardless of the ion, changes of the chemical potential induced by an increase in pressure appear to be dominated by the hydrophobic component, in particular when using the alternative correction. For bromide and iodide electrolytes, the two corrections give chemical potentials in good agreement.

Adapting free energy perturbation simulations for large macrocyclic ligands: how to dissect contributions from direct binding and free ligand flexibility

Large and flexible ligands gain increasing interest in the development of bioactive agents. They challenge the applicability of computational ligand optimization strategies originally developed for small molecules. Free energy perturbation (FEP) is often used for predicting binding affinities of small molecule ligands, however, its use for more complex ligands remains limited. Herein, we report the structure-based design of peptide macrocycles targeting the protein binding site of human adaptor protein 14-3-3. We observe a surprisingly strong dependency of binding affinities on relatively small variations in substituent size. FEP was performed to rationalize observed trends. To account for insufficient convergence of FEP, restrained calculations were performed and complemented with extensive REST MD simulations of the free ligands. These calculations revealed that changes in affinity originate both from altered direct interactions and conformational changes of the free ligand. In addition, MD simulations provided the basis to rationalize unexpected trends in ligand lipophilicity. We also verified the anticipated interaction site and binding mode for one of the high affinity ligands by X-ray crystallography. The introduced fully-atomistic simulation protocol can be used to rationalize the development of structurally complex ligands which will support future ligand maturation efforts.

Computational methodology for solubility prediction: Application to sparingly soluble organic/inorganic materials

The solubility of a crystalline material can be estimated from the absolute free energy of the solid and the excess solvation free energy. In the earlier work, we presented a general-purpose molecular-dynamics-based methodology enabling solubility predictions of crystalline compounds, yielding accurate estimates of the aqueous solubilities of naphthalene at various pressures and temperatures. In the present work, we investigate a number of prototypical complex materials, including phenanthrene, calcite, and aragonite over a range of temperatures and pressures. Our results provide stronger evidence for the power of the methodology for universal solubility predictions.

Computational methodology for solubility prediction: Application to the sparingly soluble solutes

The solubility of a crystalline substance in the solution can be estimated from its absolute solid free energy and excess solvation free energy. Here, we present a numerical method, which enables convenient solubility estimation of general molecular crystals at arbitrary thermodynamic conditions where solid and solution can coexist. The methodology is based on standard alchemical free energy methods, such as thermodynamic integration and free energy perturbation, and consists of two parts: (1) systematic extension of the Einstein crystal method to calculate the absolute solid free energies of molecular crystals at arbitrary temperatures and pressures and (2) a flexible cavity method that can yield accurate estimates of the excess solvation free energies. As an illustration, via classical Molecular Dynamic simulations, we show that our approach can predict the solubility of OPLS-AA-based (Optimized Potentials for Liquid Simulations All Atomic) naphthalene in SPC (Simple Point Charge) water in good agreement with experimental data at various temperatures and pressures. Because the procedure is simple and general and only makes use of readily available open-source software, the methodology should provide a powerful tool for universal solubility prediction.

Excess Volumes from the Pressure Derivative of the Excess Chemical Potential: Testing Simple Models for Cavity Formation in Water

Excess volumes and excess compressibilities for hard spheres in water were computed by pressure derivatives of the excess chemical potential, which is equivalent to the work of cavity formation. This is relevant to the application of continuum solvation methods at various pressures. The excess chemical potential was modeled within phenomenological expressions for curved surfaces plus a pressure-volume term, for which two approaches were adopted, differing for the radius of the spherical volume. This implies a different dependence on pressure of parameters. In all cases, in the surface term, for the pressure derivative of parameters of the curvature function, use was made of the previously proposed expressions for the first two moments obtained from the density and radial distribution of oxygens in liquid water. Only for the parameter which has the dimension of surface tension (γ̃) was explicit dependence on pressure considered and results are affected by the specific polynomial used. In agreement with what inferred from simulation results obtained for cavities in TIP4P water, negative and positive adsorptions at the contact radius were extrapolated for a very large cavity at 1 and 8000 atm, respectively. The expressions here employed for the excess chemical potential predict the zero value of asymptotic adsorption to be at a pressure between 500 and 800 atm, which can be compared to results from the revised scaled particle theory. In the same range, for a nanometer-sized cavity, a change of behavior occurs regarding the ratio between the excess Helmholtz free energy and the product between pressure and excess volume.

1.14*CM1A-LBCC: Localized Bond-Charge Corrected CM1A Charges for Condensed-Phase Simulations

The quality of the 1.14*CM1A and 1.20*CM5 charge models was evaluated for calculations of free energies of hydration. For a set of 426 neutral molecules, 1.14*CM1A and 1.20*CM5 yield MADs of 1.26 and 1.21 kcal/mol, respectively. The 1.14*CM1A charges, which can be readily obtained for large systems, exhibit large deviations only for a subset of functional groups. The results for these cases were systematically improved using localized bond-charge corrections (LBCC) by which offsetting adjustments are made to the partial charges for atoms in specified bond types. Only 19 LBCCs were needed to yield 1.14*CM1A-LBCC charges that reduce the errors for the 426 ΔG_hyd values to only 0.61 kcal/mol. The modified charge method was also tested in computation of heats of vaporization and densities for pure organic liquids, yielding average errors of 1.40 kcal/mol and 0.024 g/cm³, similar to those for 1.14*CM1A.

Application of a BOSS-Gaussian interface for QM/MM simulations of Henry and methyl transfer reactions

Hybrid quantum mechanics and molecular mechanics (QM/MM) computer simulations have become an indispensable tool for studying chemical and biological phenomena for systems too large to treat with QM alone. For several decades, semiempirical QM methods have been used in QM/MM simulations. However, with increased computational resources, the introduction of ab initio and density function methods into on-the-fly QM/MM simulations is being increasingly preferred. This adaptation can be accomplished with a program interface that tethers independent QM and MM software packages. This report introduces such an interface for the BOSS and Gaussian programs, featuring modification of BOSS to request QM energies and partial atomic charges from Gaussian. A customizable C-shell linker script facilitates the interprogram communication. The BOSS-Gaussian interface also provides convenient access to Charge Model 5 (CM5) partial atomic charges for multiple purposes including QM/MM studies of reactions. In this report, the BOSS-Gaussian interface is applied to a nitroaldol (Henry) reaction and two methyl transfer reactions in aqueous solution. Improved agreement with experiment is found by determining free-energy surfaces with MP2/CM5 QM/MM simulations than previously reported investigations using semiempirical methods.

Evaluation of CM5 Charges for Nonaqueous Condensed-Phase Modeling

Partial atomic charges for neutral molecules from quantum mechanical calculations are typically scaled for use in molecular modeling of liquid-phase systems. Optimal scale factors of 1.14 for CM1A and 1.27 for CM5 charges were previously determined for minimizing errors in free energies of hydration. The adequacy of the 1.14*CM1A and 1.27*CM5 models are evaluated here in pure liquid simulations in combination with the OPLS-AA force field. For 22 organic liquids, the 1.14*CM1A and 1.27*CM5 models yield mean unsigned errors (MUEs) of ca. 1.40 kcal/mol for heats of vaporization. Not surprisingly, this reflects overpolarization with the scale factors derived for aqueous media. Prediction of pure liquid properties using CM5 charges is optimized using a scale factor of 1.14, which reduces the MUE for heats of vaporization to 0.89 kcal/mol. However, due to the impracticality of using different scale factors in different explicit-solvent condensed-phase simulations, a universal scale factor of 1.20 emerged for CM5 charges. This provides a balance between errors in computed pure liquid properties and free energies of hydration. Computation of free energies of hydration by the GB/SA method further found that 1.20 is equally suited for use in explicit or implicit treatments of aqueous solvation. With 1.20*CM5 charges, a variety of condensed-phase simulations can be pursued while maintaining average errors of 1.0 kcal/mol in key thermodynamic properties.