Activation free energy for proton transfer in solution

Path integral molecular dynamics methods are employed to compute the free energy for proton transfer reactions for strongly hydrogen bonded systems in a polar solvent. The free energy profile is calculated using several different techniques, including: integration of the mean force acting on the proton path with its centroid constrained at different values, the integral form of the free energy calculation in the constrained-reaction-coordinate-dynamics ensemble and direct simulation of the unconstrained dynamics. The results show that estimates of the free energy barrier obtained by harmonic extrapolation are likely to be in error. Both quantum and classical results for the free energy are obtained and compared with simulations using adiabatic quantum dynamics. Comparison of the quantum and classical results show that there are quantum corrections to the solvent contributions to the free energy.     

Blind prediction of distribution in the SAMPL5 challenge with QM based protomer and p<Emphasis Type="Italic">K</Emphasis><Subscript>a</Subscript> corrections

The computation of distribution coefficients between polar and apolar phases requires both an accurate characterization of transfer free energies between phases and proper accounting of ionization and protomerization. We present a protocol for accurately predicting partition coefficients between two immiscible phases, and then apply it to 53 drug-like molecules in the SAMPL5 blind prediction challenge. Our results combine implicit solvent QM calculations with classical MD simulations using the non-Boltzmann Bennett free energy estimator. The OLYP/DZP/SMD method yields predictions that have a small deviation from experiment (RMSD = 2.3 \(\log\) D units), relative to other participants in the challenge. Our free energy corrections based on QM protomer and \({\text{p}}K_{\text{a}}\) calculations increase the correlation between predicted and experimental distribution coefficients, for all methods used. Unfortunately, these corrections are overly hydrophilic, and fail to account for additional effects such as aggregation, water dragging and the presence of polar impurities in the apolar phase. We show that, although expensive, QM-NBB free energy calculations offer an accurate and robust method that is superior to standard MM and QM techniques alone.     

Ab initio calculations of the melting temperatures of refractory bcc metals

We present ab initio calculations of the melting temperatures for bcc metals Nb, Ta and W. The calculations combine phase coexistence molecular dynamics (MD) simulations using classical embedded-atom method potentials and ab initio density functional theory free energy corrections. The calculated melting temperatures for Nb, Ta and W are, respectively, within 3%, 4%, and 7% of the experimental values. We compare the melting temperatures to those obtained from direct ab initio molecular dynamics simulations and see if they are in excellent agreement with each other. The small remaining discrepancies with experiment are thus likely due to inherent limitations associated with exchange-correlation energy approximations within density-functional theory.     

Revisiting the Frenkel-Ladd method to compute the free energy of solids: the Einstein molecule approach

In this paper a new method to evaluate the free energy of solids is proposed. The method can be regarded as a variant of the method proposed by Frenkel and Ladd [J. Chem. Phys. 81, 3188 (1984)]. The main equations of the method can be derived in a simple way. The method can be easily implemented within a Monte Carlo program. We have applied the method to determine the free energy of hard spheres in the solid phase for several system sizes. The obtained free energies agree within the numerical uncertainty with those obtained by Polson et al. [J. Chem. Phys. 112, 5339 (2000)]. The fluid-solid equilibria has been determined for several system sizes and compared to the values published previously by Wilding and Bruce [Phys. Rev. Lett. 85, 5138 (2000)] using the phase switch methodology. It is shown that both the free energies and the coexistence pressures present a strong size dependence and that the results obtained from free energy calculations agree with those obtained using the phase switch method, which constitutes a cross-check of both methodologies. From the results of this work we estimate the coexistence pressure of the fluid-solid transition of hard spheres in the thermodynamic limit to be p*=11.54(4), which is slightly lower than the classical value of Hoover and Ree (p*=11.70) [J. Chem. Phys. 49, 3609 (1968)]. Taking into account the strong size dependence of the free energy of the solid phase, we propose to introduce finite size corrections, which allow us to estimate approximately the free energy of the solid phase in the thermodynamic limit from the known value of the free energy of the solid phase with N molecules. We have also determined the free energy of a Lennard-Jones solid by using both the methodology of this work and the finite size correction. It is shown how a relatively good estimate of the free energy of the system in the thermodynamic limit is obtained even from the free energy of a relatively small system.     

Thermodynamic integration from classical to quantum mechanics

We present a new method for calculating quantum mechanical corrections to classical free energies, based on thermodynamic integration from classical to quantum mechanics. In contrast to previous methods, our method is numerically stable even in the presence of strong quantum delocalization. We first illustrate the method and its relationship to a well-established method with an analysis of a one-dimensional harmonic oscillator. We then show that our method can be used to calculate the quantum mechanical contributions to the free energies of ice and water for a flexible water model, a problem for which the established method is unstable.     

Semiclassical statistical mechanics of two-dimensional hard-body fluids

The problem of calculating the thermodynamic properties of two-dimensional semiclassical hard-body fluids is studied. Explicit expressions are given for the first-order quantum corrections to the free energy, equation of state, and virial coefficients. The numerical results are calculated for the planar hard dumbbell fluid. Significant features are the increase in quantum corrections with increasing eta and increasing L*=L/sigma(0).     

The generalized Born/volume integral implicit solvent model: estimation of the free energy of hydration using London dispersion instead of atomic surface area

A new generalized Born model for estimating the free energy of hydration is presented. The new generalized Born/volume integral (GB/VI) estimates the free energy of hydration as a classical electrostatic energy plus a cavitation energy that is not based upon atomic surface area (SA) used in GB/SA hydration models but on a VI London dispersion energy estimated from quantities already calculated in the classical electrostatic energy. The (relatively few) GB/VI model parameters are fitted to experimental data, and parameterizations for two different atomic partial charge models are presented. Comparison of the calculated and experimental free energies of hydration for 560 small molecules (both neutral and charged) shows good agreement (r(2) = 0.94).     

Confinement free energy of flexible polyelectrolytes in spherical cavities

A weakly charged flexible polyelectrolyte chain in a neutral spherical cavity is analyzed by using self-consistent field theory within an explicit solvent model. Assuming the radial symmetry for the system, it is found that the confinement of the chain leads to creation of a charge density wave along with the development of a potential difference across the center of cavity and the surface. We show that the solvent entropy plays an important role in the free energy of the confined system. For a given radius of the spherical cavity and fixed charge density along the backbone of the chain, solvent and small ion entropies dominate over all other contributions when chain lengths are small. However, with the increase in chain length, chain conformational entropy and polymer-solvent interaction energy also become important. Our calculations reveal that energy due to electrostatic interactions plays a minor role in the free energy. Furthermore, we show that the total free energy under spherical confinement is not extensive in the number of monomers. Results for the osmotic pressure and mean activity coefficient for monovalent salt are presented. We demonstrate that fluctuations at one-loop level lower the free energy and corrections to the osmotic pressure and mean activity coefficient of the salt are discussed. Finite size corrections are shown to widen the range of validity of the fluctuation analysis.     

Offcenter Molecular Charge Distribution Corrections to Born Equation Electrostatic Solvation Energies

Use of simple model charge distributions and classical multipolar dielectric theory has allowed investigation of the corrections appropriate to the Born equation for structured ion solutes. Among the results obtained for centrosymmetric charge distributions, but presumably of general consequence, are the demonstration of nonlinear dependence and great sensitivity of polarization to charge-cavity boundary proximity, rapid falloff of effects upon continued subdivision of charge, and the significance of solvation energy differences for 2- and 3-dimensional charge distributions. That 50% or greater corrections to the Born energy may readily be obtained in real ion systems, and the lack of sensitivity of such corrections to the dielectric of the medium for ε > 10 are also results of general interest. Analytical approximations for the electrostatic work equations appropriate to high dielectric solvation are included, and how uncertainties in charge magnitude and position for real ions might affect the analysis are also briefly considered.     

A comparative analysis of models for hydration free energy calculations using the RISM theory of integral equations

Rigorous calculations and a detailed analysis of the free energies of hydration were performed using the RISM (reference interaction site model) approach for 96 compounds of various chemical natures. A comparison of all the existing models for calculation of hydration free energy, including models that use corrections and semiempirical parameters, was performed for the first time. The applicability ranges of all the models under consideration were determined. It was shown that the PWC model based on jointly using the RISM and chemoinformatics approaches gave most accurate hydration free energy values. This model allows the free energy of hydration to be predicted with high accuracy and, at the same time, does not require the use of substantial computer resources.     

FRee electron model for heteroatomic conjugate molecules

A simple generalization of the free electron model is suggested for describing the heteroatomic conjugate molecules and obtaining analytical equations for the energy levels of these systems. The δ-like potential is used as a potential of a heteroatom in this model, and the energy levels are obtained as corrections for first order perturbation theory.     

The oxidation of tyrosine and tryptophan studied by a molecular dynamics normal hydrogen electrode

The thermochemical constants for the oxidation of tyrosine and tryptophan through proton coupled electron transfer in aqueous solution have been computed applying a recently developed density functional theory (DFT) based molecular dynamics method for reversible elimination of protons and electrons. This method enables us to estimate the solvation free energy of a proton (H(+)) in a periodic model system from the free energy for the deprotonation of an aqueous hydronium ion (H(3)O(+)). Using the computed solvation free energy of H(+) as reference, the deprotonation and oxidation free energies of an aqueous species can be converted to pK(a) and normal hydrogen electrode (NHE) potentials. This conversion requires certain thermochemical corrections which were first presented in a similar study of the oxidation of hydrobenzoquinone [J. Cheng, M. Sulpizi, and M. Sprik, J. Chem. Phys. 131, 154504 (2009)]. Taking a different view of the thermodynamic status of the hydronium ion, these thermochemical corrections are revised in the present work. The key difference with the previous scheme is that the hydronium is now treated as an intermediate in the transfer of the proton from solution to the gas-phase. The accuracy of the method is assessed by a detailed comparison of the computed pK(a), NHE potentials and dehydrogenation free energies to experiment. As a further application of the technique, we have analyzed the role of the solvent in the oxidation of tyrosine by the tryptophan radical. The free energy change computed for this hydrogen atom transfer reaction is very similar to the gas-phase value, in agreement with experiment. The molecular dynamics results however, show that the minimal solvent effect on the reaction free energy is accompanied by a significant reorganization of the solvent.     

Blind prediction of SAMPL4 cucurbit[7]uril binding affinities with the mining minima method

Accurate methods for predicting protein–ligand binding affinities are of central interest to computer-aided drug design for hit identification and lead optimization. Here, we used the mining minima (M2) method to predict cucurbit[7]uril binding affinities from the SAMPL4 blind prediction challenge. We tested two different energy models, an empirical classical force field, CHARMm with VCharge charges, and the Poisson–Boltzmann surface area solvation model; and a semiempirical quantum mechanical (QM) Hamiltonian, PM6-DH+, coupled with the COSMO solvation model and a surface area term for nonpolar solvation free energy. Binding affinities based on the classical force field correlated strongly with the experiments with a correlation coefficient (R²) of 0.74. On the other hand, binding affinities based on the QM energy model correlated poorly with experiments (R² = 0.24), due largely to two major outliers. As we used extensive conformational search methods, these results point to possible inaccuracies in the PM6-DH+ energy model or the COSMO solvation model. Furthermore, the different binding free energy components, solute energy, solvation free energy, and configurational entropy showed significant deviations between the classical M2 and quantum M2 calculations. Comparison of different classical M2 free energy components to experiments show that the change in the total energy, i.e. the solute energy plus the solvation free energy, is the key driving force for binding, with a reasonable correlation to experiment (R² = 0.56); however, accounting for configurational entropy further improves the correlation.     

Using Selectively Applied Accelerated Molecular Dynamics to Enhance Free Energy Calculations

Accelerated molecular dynamics (aMD) has been shown to enhance conformational space sampling relative to classical molecular dynamics; however, the exponential reweighting of aMD trajectories, which is necessary for the calculation of free energies relating to the classical system, is oftentimes problematic, especially for systems larger than small poly peptides. Here, we propose a method of accelerating only the degrees of freedom most pertinent to sampling, thereby reducing the total acceleration added to the system and improving the convergence of calculated ensemble averages, which we term selective aMD. Its application is highlighted in two biomolecular cases. First, the model system alanine dipeptide is simulated with classical MD, all-dihedral aMD, and selective aMD, and these results are compared to the infinite sampling limit as calculated with metadynamics. We show that both forms of aMD enhance the convergence of the underlying free energy landscape by 5-fold relative to classical MD; however, selective aMD can produce improved statistics over all-dihedral aMD due to the improved reweighting. Then we focus on the pharmaceutically relevant case of computing the free energy of the decoupling of oseltamivir in the active site of neuraminidase. Results show that selective aMD greatly reduces the cost of this alchemical free energy transformation, whereas all-dihedral aMD produces unreliable free energy estimates.     

Free energy partitioning analysis of the driving forces that determine ion density profiles near the water liquid-vapor interface

Free energy partitioning analysis is employed to explore the driving forces for ions interacting with the water liquid-vapor interface using recently optimized point charge models for the ions and SPC/E water. The Na(+) and I(-) ions are examined as an example kosmotrope/chaotrope pair. The absolute hydration free energy is partitioned into cavity formation, attractive van der Waals, local electrostatic, and far-field electrostatic contributions. We first compute the bulk hydration free energy of the ions, followed by the free energy to insert the ions at the center of a water slab. Shifts of the ion free energies occur in the slab geometry consistent with the SPC/E surface potential of the water liquid-vapor interface. Then the free energy profiles are examined for ion passage from the slab center to the dividing surface. The profiles show that, for the large chaotropic I(-) ion, the relatively flat total free energy profile results from the near cancellation of several large contributions. The far-field electrostatic part of the free energy, largely due to the water liquid-vapor interface potential, has an important effect on ion distributions near the surface in the classical model. We conclude, however, that the individual forms of the local and far-field electrostatic contributions are expected to be model dependent when comparing classical and quantum results. The substantial attractive cavity free energy contribution for the larger I(-) ion suggests that there is a hydrophobic component important for chaotropic ion interactions with the interface.     

On the Relative Merits of Equilibrium and Non‐Equilibrium Simulations for the Estimation of Free‐Energy Differences

The possibility of estimating equilibrium free‐energy profiles from multiple non‐equilibrium simulations using the fluctuation–dissipation theory or the relation proposed by Jarzynski has attracted much attention. Although the Jarzynski estimator has poor convergence properties for simulations far from equilibrium, corrections have been derived for cases in which the work is Gaussian distributed. Here, we examine the utility of corrections proposed by Gore and collaborators using a simple dissipative system as a test case. The system consists of a single methane‐like particle in explicit water. The Jarzynski equality is used to estimate the change in free energy associated with pulling the methane particle a distance of 3.9 nm at rates ranging from ～0.1 to 100 m s^?1. It is shown that although the corrections proposed by Gore and collaborators have excellent numerical performance, the profiles still converge slowly. Even when the corrections are applied in an ideal case where the work distribution is necessarily Gaussian, performing simulations under quasi‐equilibrium conditions is still most efficient. Furthermore, it is shown that even for a single methane molecule in water, pulling rates as low as 1 m s^?1 can be problematic. The implications of this finding for studies in which small molecules or even large biomolecules are pulled through inhomogeneous environments at similar pulling rates are discussed.     

Comparisons of classical chemical dynamics simulations of the unimolecular decomposition of classical and quantum microcanonical ensembles

Previous studies have shown that classical trajectory simulations often give accurate results for short-time intramolecular and unimolecular dynamics, particularly for initial non-random energy distributions. To obtain such agreement between experiment and simulation, the appropriate distributions must be sampled to choose initial coordinates and momenta for the ensemble of trajectories. If a molecule's classical phase space is sampled randomly, its initial decomposition will give the classical anharmonic microcanonical (RRKM) unimolecular rate constant for its decomposition. For the work presented here, classical trajectory simulations of the unimolecular decomposition of quantum and classical microcanonical ensembles, at the same fixed total energy, are compared. In contrast to the classical microcanonical ensemble, the quantum microcanonical ensemble does not sample the phase space randomly. The simulations were performed for CH(4), C(2)H(5), and Cl(-)---CH(3)Br using both analytic potential energy surfaces and direct dynamics methods. Previous studies identified intrinsic RRKM dynamics for CH(4) and C(2)H(5), but intrinsic non-RRKM dynamics for Cl(-)---CH(3)Br. Rate constants calculated from trajectories obtained by the time propagation of the classical and quantum microcanonical ensembles are compared with the corresponding harmonic RRKM estimates to obtain anharmonic corrections to the RRKM rate constants. The relevance and accuracy of the classical trajectory simulation of the quantum microcanonical ensemble, for obtaining the quantum anharmonic RRKM rate constant, is discussed.     

Quantum versus classical electron transfer energy as reaction coordinate for the aqueous Ru2+/Ru3+ redox reaction

Applying density functional theory (DFT)-based molecular dynamics simulation methods we investigate the effect of explicit treatment of electronic structure on the solvation free energy of aqueous Ru²⁺ and Ru³⁺.Our approach is based on the Marcus theory of redox half reactions, focussing on the vertical energy gap for reduction or oxidation of a single aqua ion. We compare the fluctuations of the quantum and classical energy gap along the same equilibrium ab initio molecular dynamics trajectory for each oxidation state. The classical gap is evaluated using a standard point charge model for the charge distribution of the solvent molecules (water). The quantum gap is computed from the full DFT electronic ground state energies of reduced and oxidized species, thereby accounting for the delocalization of the electron in the donor orbital and reorganization of the electron cloud after electron transfer (ET). The fluctuations of the quantum ET energy are well approximated by gaussian statistics giving rise to parabolic free energy profiles. The curvature is found to be independent of the oxidation state in agreement with the linear response assumption underlying Marcus theory. By contrast, the diabatic free energy curves evaluated using the classical gap as order parameter, while also quadratic, are asymmetric reflecting the difference in oxidation state. The response of these two order parameters is further analysed by a comparison of the spectral density of the fluctuations and the corresponding reorganization free energies.     

Calculation of solid-liquid interfacial free energy: a classical nucleation theory based approach

We present a simple approach to calculate the solid-liquid interfacial free energy. This new method is based on the classical nucleation theory. Using the molecular dynamics simulation, we employ spherical crystal nuclei embedded in the supercooled liquids to create an ideal model of a homogeneous nucleation. The interfacial free energy is extracted by fitting the relation between the critical nucleus size and the reciprocal of the critical undercooling temperature. The orientationally averaged interfacial free energy is found to be 0.302+/-0.002 (in standard LJ unit). The temperature dependence of the interfacial free energy is also obtained in this work. We find that the interfacial free energy increases slightly with increasing temperature. The positive temperature coefficient of the interfacial free energy is in qualitative agreement with Spaepen's analysis [Solid State Phys. 47, FS181 (1994)] and Turnbull's empirical estimation [J. Appl. Phys. 21, 1022 (1950)].