Comparison of two simulation methods to compute solvation free energies and partition coefficients

The thermodynamic integration (TI) and expanded ensemble (EE) methods are used here to calculate the hydration free energy in water, the solvation free energy in 1‐octanol, and the octanol‐water partition coefficient for a six compounds of varying functionality using the optimized potentials for liquid simulations (OPLS) all‐atom (AA) force field parameters and atomic charges. Both methods use the molecular dynamics algorithm as a primary component of the simulation protocol, and both have found wide applications in fields such as the calculation of activity coefficients, phase behavior, and partition coefficients. Both methods result in solvation free energies and 1‐octanol/water partition coefficients with average absolute deviations (AAD) from experimental data to within 4 kJ/mol and 0.5 log units, respectively. Here, we find that in simulations the OPLS‐AA force field parameters (with fixed charges) can reproduce solvation free energies of solutes in 1‐octanol with AAD of about half that for the solute hydration free energies using a extended simple point charge (SPC/E) model of water. The computational efficiency of the two simulation methods are compared based on the time (in nanoseconds) required to obtain similar standard deviations in the solvation free energies and 1‐octanol/water partition coefficients. By this analysis, the EE method is found to be a factor of nine more efficient than the TI algorithm. For both methods, solvation free energy calculations in 1‐octanol consume roughly an order of magnitude more CPU hours than the hydration free energy calculations. © 2012 Wiley Periodicals, Inc.     

Free-energy analysis of the molecular binding into lipid membrane with the method of energy representation

A statistical-mechanical treatment of the molecular binding into lipid membrane is presented in combination with molecular simulation. The membrane solution is viewed as an inhomogeneous, mixed solvent system, and the free energy of solvation of a solute in membrane is computed with a realistic set of potential functions by the method of energy representation. Carbon monoxide, carbon dioxide, benzene, and ethylbenzene are adopted as model solutes to analyze the binding into 1,2-dimyristoyl-sn-glycero-3-phosphatidylcholine (DMPC) membrane. It is shown that the membrane inside is more favorable than bulk water and that the solute distribution is diffuse throughout the membrane inside. The membrane-water partition coefficient is then constructed with the help of the Kirkwood-Buff theory from the solvation free energy obtained separately in the hydrophobic, glycerol, headgroup, and aqueous regions. To discuss the role of repulsive and attractive interactions, the solvation free energy is partitioned into the DMPC and water contributions and the effect of water to stabilize the benzene and ethylbenzene solutes within the membrane is pointed out.     

Coupled semiempirical quantum mechanics and molecular mechanics (QM/MM) calculations on the aqueous solvation free energies of ionized molecules

The aqueous solvation free energies of ionized molecules were computed using a coupled quantum mechanical and molecular mechanical (QM/MM) model based on the AM1, MNDO, and PM3 semiempirical molecular orbital methods for the solute molecule and the TIP3P molecular mechanics model for liquid water. The present work is an extension of our model for neutral solutes where we assumed that the total free energy is the sum of components derived from the electrostatic/polarization terms in the Hamiltonian plus an empirical “nonpolar” term. The electrostatic/polarization contributions to the solvation free energies were computed using molecular dynamics (MD) simulation and thermodynamic integration techniques, while the nonpolar contributions were taken from the literature. The contribution to the electrostatic/polarization component of the free energy due to nonbonded interactions outside the cutoff radii used in the MD simulations was approximated by a Born solvation term. The experimental free energies were reproduced satisfactorily using variational parameters from the vdW terms as in the original model, in addition to a parameter from the one-electron integral terms. The new one-electron parameter was required to account for the short-range effects of overlapping atomic charge densities. The radial distribution functions obtained from the MD simulations showed the expected H-bonded structures between the ionized solute molecule and solvent molecules. We also obtained satisfactory results by neglecting both the empirical nonpolar term and the electronic polarization of the solute, i.e., by implementing a nonpolarization model. ©1999 John Wiley & Sons, Inc. J Comput Chem 20: 1028–1038, 1999     

Pathway dependence of the efficiency of calculating free energy and entropy of solute–solute association in water

In this study we investigate two alternative pathways to compute the free energy and the entropy of small molecule association (ΔF_assoc and ΔS_assoc) in water. The first route (direct pathway) uses thermodynamic integration as function of the distance R between the solutes. The mean force and the mean covariance of the force with the energy in solution are calculated from molecular dynamics simulation followed by integration of these quantities with respect to the reaction coordinate R. The alternative approach examined (solvation pathway) would first remove the solutes from the solution using thermodynamic integration as function of a solvation coupling parameter λ, change the solute–solute distance in vacuo and then solvate back the solute pair at the new separation distance. The system studied was a pair of CH₄ molecules in water. We investigate the influence of the CH₄–water interaction strength on the obtained ΔF_assoc and ΔS_assoc values by changing van der Waals and Coulomb interaction and evaluated the accuracy and efficiency for the two pathways. We find that the direct route seems more suitable for the calculation of free energies of hydrophobic solutes while the solvation pathway performs better when calculating entropy changes for solutes that have a stronger interaction with the solvent.     

New approach to free energy of solvation applying continuum models to molecular dynamics simulation

A new approach to the calculation of the free energy of solvation from trajectories obtained by molecular dynamics simulation is presented. The free energy of solvation is computed as the sum of three contributions originated at the cavitation of the solute by the solvent, the solute-solvent nonpolar (repulsion and dispersion) interactions, and the electrostatic solvation of the solute. The electrostatic term is calculated based on ideas developed for the broadly used continuum models, the cavitational contribution from the excluded volume by the Claverie-Pierotti model, and the Van der Waals term directly from the molecular dynamics simulation. The proposed model is tested for diluted aqueous solutions of simple molecules containing a variety of chemically important functions: methanol, methylamine, water, methanethiol, and dichloromethane. These solutions were treated by molecular dynamics simulations using SPC/E water and the OPLS force field for the organic molecules. Obtained free energies of solvation are in very good agreement with experimental data.     

Electrostatic solvation free energy of amino acid side chain analogs: implications for the validity of electrostatic linear response in water

Electrostatic free energies of solvation for 15 neutral amino acid side chain analogs are computed. We compare three methods of varying computational complexity and accuracy for three force fields: free energy simulations, Poisson-Boltzmann (PB), and linear response approximation (LRA) using AMBER, CHARMM, and OPLS-AA force fields. We find that deviations from simulation start at low charges for solutes. The approximate PB and LRA produce an overestimation of electrostatic solvation free energies for most of molecules studied here. These deviations are remarkably systematic. The variations among force fields are almost as large as the variations found among methods. Our study confirms that success of the approximate methods for electrostatic solvation free energies comes from their ability to evaluate free energy differences accurately.     

Breaking the polar‐nonpolar division in solvation free energy prediction

Implicit solvent models divide solvation free energies into polar and nonpolar additive contributions, whereas polar and nonpolar interactions are inseparable and nonadditive. We present a feature functional theory (FFT) framework to break this ad hoc division. The essential ideas of FFT are as follows: (i) representability assumption: there exists a microscopic feature vector that can uniquely characterize and distinguish one molecule from another; (ii) feature‐function relationship assumption: the macroscopic features, including solvation free energy, of a molecule is a functional of microscopic feature vectors; and (iii) similarity assumption: molecules with similar microscopic features have similar macroscopic properties, such as solvation free energies. Based on these assumptions, solvation free energy prediction is carried out in the following protocol. First, we construct a molecular microscopic feature vector that is efficient in characterizing the solvation process using quantum mechanics and Poisson–Boltzmann theory. Microscopic feature vectors are combined with macroscopic features, that is, physical observable, to form extended feature vectors. Additionally, we partition a solvation dataset into queries according to molecular compositions. Moreover, for each target molecule, we adopt a machine learning algorithm for its nearest neighbor search, based on the selected microscopic feature vectors. Finally, from the extended feature vectors of obtained nearest neighbors, we construct a functional of solvation free energy, which is employed to predict the solvation free energy of the target molecule. The proposed FFT model has been extensively validated via a large dataset of 668 molecules. The leave‐one‐out test gives an optimal root‐mean‐square error (RMSE) of 1.05 kcal/mol. FFT predictions of SAMPL0, SAMPL1, SAMPL2, SAMPL3, and SAMPL4 challenge sets deliver the RMSEs of 0.61, 1.86, 1.64, 0.86, and 1.14 kcal/mol, respectively. Using a test set of 94 molecules and its associated training set, the present approach was carefully compared with a classic solvation model based on weighted solvent accessible surface area. © 2017 Wiley Periodicals, Inc.     

Theory of solvation in polar nematics

We develop a linear response theory of solvation of ionic and dipolar solutes in anisotropic, axially symmetric polar solvents. The theory is applied to solvation in polar nematic liquid crystals. The formal theory constructs the solvation response function from projections of the solvent dipolar susceptibility on rotational invariants. These projections are obtained from Monte Carlo simulations of a fluid of dipolar spherocylinders which can exist both in the isotropic and nematic phases. Based on the properties of the solvent susceptibility from simulations and the formal solution, we have obtained a formula for the solvation free energy which incorporates the experimentally available properties of nematics and the length of correlation between the dipoles in the liquid crystal. The theory provides a quantitative framework for analyzing the steady-state and time-resolved optical spectra and makes several experimentally testable predictions. The equilibrium free energy of solvation, anisotropic in the nematic phase, is given by a quadratic function of cosine of the angle between the solute dipole and the solvent nematic director. The sign of solvation anisotropy is determined by the sign of dielectric anisotropy of the solvent: solvation anisotropy is negative in solvents with positive dielectric anisotropy and vice versa. The solvation free energy is discontinuous at the point of isotropic-nematic phase transition. The amplitude of this discontinuity is strongly affected by the size of the solute becoming less pronounced for larger solutes. The discontinuity itself and the magnitude of the splitting of the solvation free energy in the nematic phase are mostly affected by microscopic dipolar correlations in the nematic solvent. Illustrative calculations are presented for the equilibrium Stokes shift and the Stokes shift time correlation function of coumarin-153 in 4-n-pentyl-4'-cyanobiphenyl and 4,4-n-heptyl-cyanopiphenyl solvents as a function of temperature in both the nematic and isotropic phases.     

Coupled reference interaction site model/simulation approach for thermochemistry of solvation: theory and prospects

We present a new methodology for computing solvation free energy, which is based upon the reference interaction site model (RISM)/hypernetted chain (HNC) solvation free energy expression, but which substitutes radial distribution functions taken from simulations for those calculated by simultaneous solution of the RISM and HNC equations. Consequently, solvation free energy can be obtained from a single molecular dynamics or Monte Carlo simulation. Here we describe in detail the coupled RISM/simulation approach, and offer some error analysis. Finally we give the results of its application to a set of small test molecules in aqueous solution. The success shown in some of our results demonstrates that the coupled RISM/simulation approach is worth considering further as a potentially useful tool in studies of solvated systems, such as aqueous molecular biosystems.     

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.     

Dynamics of water trapped between hydrophobic solutes

We describe the model dynamical behavior of the solvent between two nanoscopic hydrophobic solutes. The dynamics of the vicinal water in various sized traps is found to be significantly different from bulk behavior. We consider the dynamics at normal temperature and pressure at three intersolute distances corresponding to the three solvent separated minima in the free energy profile between the solutes with attractions. These three states correspond to one, two, and three intervening layers of water molecules. Results are obtained from a molecular dynamics simulation at constant temperature and pressure (NPT) ensemble. Translational diffusion of water molecules trapped between the two solutes has been analyzed from the velocity correlation function as well as from the mean square displacement of the water molecules. The rotational behavior has been analyzed through the reorientational dynamics of the dipole moment vector of the water molecule by calculating both first and second rank dipole-dipole correlation functions. Both the translational and reorientational mobilities of water are found to be much slower at the smaller separation and increases as the separation between solutes becomes larger. The occupation time distribution functions calculated from the trajectories also show that the relaxation is much slower for the smallest intersolute separation as compared to other wider separations. The sublinear trend in mean square displacement and the stretched exponential decay of the relaxation of dipolar correlation and occupation distribution function indicate that the dynamical behavior of water in the confined region between two large hydrophobic solutes departs from usual Brownian behavior. This behavior is reminiscent of the behavior of water in the vicinity of protein surface clefts or trapped between two domains of a protein.     

Universal solvation model based on conductor‐like screening model

Atomic surface tensions are parameterized for use with solvation models in which the electrostatic part of the calculation is based on the conductor‐like screening model (COSMO) and the semiempirical molecular orbital methods AM1, PM3, and MNDO/d. The convergence of the calculated polarization free energies with respect to the numerical parameters of the electrostatic calculations is first examined. The accuracy and precision of the calculated values are improved significantly by adjusting two parameters that control the segmentation of the solvent‐accessible surface that is used for the calculations. The accuracy of COSMO calculations is further improved by adopting an optimized set of empirical electrostatic atomic radii. Finally, the electrostatic calculation is combined with SM5‐type atomic surface tension functionals that are used to compute the nonelectrostatic portions of the solvation free energy. All parameterizations are carried out using rigid (R) gas‐phase geometries; this combination (SM5‐type surface tensions, COSMO electrostatics, and rigid geometries) is called SM5CR. Six air–water and 76 water–solvent partition coefficients are added to the training set of air–solvent data points previously used to parameterize the SM5 suite of solvation models, thereby bringing the total number of data points in the training set to 2266. The model yields free energies of solvation and transfer with mean unsigned errors of 0.63, 0.59, and 0.61 kcal/mol for AM1, PM3, and MNDO/d, respectively, over all 2217 data points for neutral solutes in the training set and mean unsigned errors of 3.0, 2.7, and 3.1 kcal/mol, respectively, for 49 data points for the ions. © 2000 John Wiley & Sons, Inc. J Comput Chem 21: 340–366, 2000     

Thermodynamics of solvation and solvophobic effect in formamide

Using semi-adiabatic calorimetry, we measured the enthalpies of solution for various low-polar compounds including alkanes, aromatic hydrocarbons and their halogenated derivatives in formamide at temperature of 298 K. For the same compounds, the values of limiting activity coefficients in formamide were determined using GC headspace analysis at 298 K, and Gibbs free energies of solution and solvation were calculated. Based on these data and the available literature values of the Gibbs free energy of solvation in formamide for a number of other low-polar solutes, a study of the solvophobic effect in this solvent is performed, and its resemblance to the hydrophobic effect in aqueous solutions is demonstrated. It is shown that the contribution of the solvophobic effect into the solvation Gibbs free energy in formamide is much higher than that in aliphatic alcohols, but lower than that in water. Like in water, the magnitude of this contribution for different solutes linearly increases with the solute molecular volume. Solvophobic effect also significantly affects the enthalpies of dissolution in formamide, causing them to be more negative in the case of alkanes and more positive in the case of arenes.     

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.     

Free-energy analysis of the electron-density fluctuation in the quantum-mechanical/molecular-mechanical simulation combined with the theory of energy representation

The relationship is investigated for QM/MM (quantum-mechanical/molecular-mechanical) systems between the fluctuations of the electronic state of the QM subsystem and of the solvation effect due to the QM-MM interaction. The free-energy change due to the electron-density fluctuation around its average is highlighted, and is evaluated through an approximate functional formulated in terms of distribution functions of the many-body coupling (pairwise non-additive) part of the QM-MM interaction energy. A set of QM/MM simulations are conducted in MM water solvent for QM water solute in ambient and supercritical conditions and for QM glycine solute in the neutral and zwitterionic forms. The variation of the electronic distortion energy of the QM solute in the course of QM/MM simulation is then shown to be compensated by the corresponding variation of the free energy of solvation. The solvation free energy conditioned by the electronic distortion energy is further analyzed with its components. It is found that the many-body contribution is essentially equal between the free energy and the average sum of solute-solvent interaction energy.     

Electron transfer in a radical ion pair: quantum calculations of the solvent reorganization energy

Results are presented for an investigation of intermolecular electron transfer (ET) in solution by means of quantum calculations. The two molecules that are involved in the ET reaction form a solvent-separated radical ion pair. The solvent plays an important role in the ET between the two molecules. In particular, it can give rise to specific solute-solvent interactions with the solutes. An example of specific interactions is the formation of a hydrogen bond between a protic solvent and one of the molecules involved in the ET. We address the study of this system by means of quantum calculations on the solutes immersed in a continuum solvent. However, when the solvent can give rise to hydrogen bond formation with the negatively charged ion after ET, we explicitly consider solvent molecules in the solute cavity, determining the hydrogen bond energetic contribution to the overall interaction energy. Solute-solvent pair distribution functions, showing the different arrangement of solvent molecules before and after ET in the first solvation shell, are reported. We provide results of the solvent reorganization energy from quantum calculations for both the two isolated fragments and the ion pair in solution. Results are in agreement with available experimental data.