Single-ensemble nonequilibrium path-sampling estimates of free energy differences

We introduce a straightforward, single-ensemble, path sampling approach to calculate free energy differences based on Jarzynski's relation. For a two-dimensional "toy" test system, the new (minimally optimized) method performs roughly one hundred times faster than either optimized "traditional" Jarzynski calculations or conventional thermodynamic integration. The simplicity of the underlying formalism suggests the approach will find broad applicability in molecular systems.     

Density-dependent analysis of nonequilibrium paths improves free energy estimates II. A Feynman-Kac formalism

The nonequilibrium fluctuation theorems have paved the way for estimating equilibrium thermodynamic properties, such as free energy differences, using trajectories from driven nonequilibrium processes. While many statistical estimators may be derived from these identities, some are more efficient than others. It has recently been suggested that trajectories sampled using a particular time-dependent protocol for perturbing the Hamiltonian may be analyzed with another one. Choosing an analysis protocol based on the nonequilibrium density was empirically demonstrated to reduce the variance and bias of free energy estimates. Here, we present an alternate mathematical formalism for protocol postprocessing based on the Feynmac-Kac theorem. The estimator that results from this formalism is demonstrated on a few low-dimensional model systems. It is found to have reduced bias compared to both the standard form of Jarzynski's equality and the previous protocol postprocessing formalism.     

Bayesian estimates of free energies from nonequilibrium work data in the presence of instrument noise

The Jarzynski equality and the fluctuation theorem relate equilibrium free energy differences to nonequilibrium measurements of the work. These relations extend to single-molecule experiments that have probed the finite-time thermodynamics of proteins and nucleic acids. The effects of experimental error and instrument noise have not been considered previously. Here, we present a Bayesian formalism for estimating free energy changes from nonequilibrium work measurements that compensates for instrument noise and combines data from multiple driving protocols. We reanalyze a recent set of experiments in which a single RNA hairpin is unfolded and refolded using optical tweezers at three different rates. Interestingly, the fastest and farthest-from-equilibrium measurements contain the least instrumental noise and, therefore, provide a more accurate estimate of the free energies than a few slow, more noisy, near-equilibrium measurements. The methods we propose here will extend the scope of single-molecule experiments; they can be used in the analysis of data from measurements with atomic force microscopy, optical, and magnetic tweezers.     

Minimum free energy pathways and free energy profiles for conformational transitions based on atomistic molecular dynamics simulations

An efficient method for the calculation of minimum free energy pathways and free energy profiles for conformational transitions is presented. Short restricted perturbation-targeted molecular dynamics trajectories are used to generate an approximate free energy surface. Approximate reaction pathways for the conformational change are constructed from one-dimensional line segments on this surface using a Monte Carlo optimization. Accurate free energy profiles are then determined along the pathways by means of one-dimensional adaptive umbrella sampling simulations. The method is illustrated by its application to the alanine "dipeptide." Due to the low computational cost and memory demands, the method is expected to be useful for the treatment of large biomolecular systems.     

Protein-ligand binding affinity by nonequilibrium free energy methods

Nonequilibrium (NE) free energy methods are embarrassingly parallel and may be very conveniently run on desktop computers using distributed computing software. In recent years there has been a proliferation of NE methods, but these approaches have barely, if at all, been used in the context of calculating protein-ligand binding free energies. In a recent study by these authors, different combinations of NE methods with various test systems were compared and protocols identified which yielded results as accurate as replica exchange thermodynamic integration (RETI). The NE approaches, however, lend themselves to extensive parallelization through the use of distributed computing. Here the best performing of those NE protocols, a replica exchange method using Bennett's acceptance ratio as the free energy estimator (RENE), is applied to two sets of congeneric inhibitors bound to neuraminidase and cyclooxygenase-2. These protein-ligand systems were originally studied with RETI, giving results to which NE and RENE simulations are compared. These NE calculations were carried out on a large, highly distributed group of low-performance desktop computers which are part of a Condor pool. RENE was found to produce results of a predictive quality at least as good as RETI in less than half the wall clock time. However, non-RE NE results were found to be far less predictive. In addition, the RENE method successfully identified a localized region of rapidly changing free energy gradients without the need for prior investigation. These results suggest that the RENE protocol is appropriate for use in the context of predicting protein-ligand binding free energies and that it can offer advantages over conventional, equilibrium approaches.     

Free energy surfaces from nonequilibrium processes without work measurement

Recent developments in statistical mechanics have allowed the estimation of equilibrium free energies from the statistics of work measurements during processes that drive the system out of equilibrium. Here a different class of processes is considered, wherein the system is prepared and released from a nonequilibrium state, and no external work is involved during its observation. For such "clamp-and-release" processes, a simple strategy for the estimation of equilibrium free energies is offered. The method is illustrated with numerical simulations and analyzed in the context of tethered single-molecule experiments.     

Calculations of crystal-melt interfacial free energies by nonequilibrium work measurements

We developed a multistep thermodynamic perturbation method to compute the interfacial free energies by nonequilibrium work measurements with cleaving potential procedure. Using this method, we calculated the interfacial free energies of different crystal orientations for the Lennard-Jones system. Our results are in good agreement with the results by thermodynamic integration method. Compared with thermodynamic integration method, the multistep thermodynamic perturbation method is more efficient. For each stage of the cleaving process, only a few thermodynamic perturbation steps are needed, and there is no requirement on the reversibility of the path.     

A new algorithm for extended nonequilibrium molecular dynamics simulations of mixed flow

In this work, we develop a new algorithm for nonequilibrium molecular dynamics of fluids under planar mixed flow, a linear combination of planar elongational flow and planar Couette flow. To date, the only way of simulating mixed flow using nonequilibrium molecular dynamics techniques was to impose onto the simulation box irreversible transformations. This would bring the simulation to an end as soon as the minimum lattice space requirements were violated. In practical terms, this meant repeating the short simulations to improve statistics and extending the box dimensions to increase the total simulation time. Our method, similar to what has already been done for pure elongational flow, allows a cuboid box to deform in time following the streamlines of the mixed flow and, after a period of time determined by the elongational field, to be mapped back and recover its initial shape. No discontinuity in physical properties is present during the mapping and the simulation can, in this way, be extended indefinitely. We also show that the most general form of mixed flow, in which the angle between the expanding (or contracting) direction and the velocity gradient axis varies, can be cast in a so-called canonical form, in which the angle assumes values that are multiples of π (when a mixed flow exists), by an appropriate choice of the field parameters.     

Equilibrium and nonequilibrium molecular dynamics simulations of the thermal conductivity of molten alkali halides

The thermal conductivity of molten NaCl and KCl was calculated through the Evans-Gillan nonequilibrium molecular dynamics (NEMD) algorithm and Green-Kubo equilibrium molecular dynamics (EMD) simulations. The EMD simulations were performed for a "binary" ionic mixture and the NEMD simulations assumed a pure system for reasons discussed in this work. The cross thermoelectric coefficient obtained from Green-Kubo EMD simulations is discussed in terms of the homogeneous thermoelectric power or Seebeck coefficient of these materials. The thermal conductivity obtained from NEMD simulations is found to be in very good agreement with that obtained through Green-Kubo EMD simulations for a binary ionic mixture. This result points to a possible cancellation between the neglected "partial enthalpy" contribution to the heat flux associated with the interdiffusion of one species through the other and that part of the thermal conductivity related to the coupled fluxes of charge and heat in "binary" ionic mixtures.     

Determination of free energy profiles by repository based adaptive umbrella sampling: bridging nonequilibrium and quasiequilibrium simulations

We propose a new adaptive sampling approach to determine free energy profiles with molecular dynamics simulations, which is called as "repository based adaptive umbrella sampling" (RBAUS). Its main idea is that a sampling repository is continuously updated based on the latest simulation data, and the accumulated knowledge and sampling history are then employed to determine whether and how to update the biasing umbrella potential for subsequent simulations. In comparison with other adaptive methods, a unique and attractive feature of the RBAUS approach is that the frequency for updating the biasing potential depends on the sampling history and is adaptively determined on the fly, which makes it possible to smoothly bridge nonequilibrium and quasiequilibrium simulations. The RBAUS method is first tested by simulations on two simple systems: a double well model system with a variety of barriers and the dissociation of a NaCl molecule in water. Its efficiency and applicability are further illustrated in ab initio quantum mechanics/molecular mechanics molecular dynamics simulations of a methyl-transfer reaction in aqueous solution.     

Accuracy and convergence of free energy differences calculated from nonequilibrium switching processes

The molecular-dynamics-based calculation of accurate free energy differences for biomolecular systems is a challenging task. Accordingly, convergence and accuracy of established equilibrium methods has been subject of many studies, often focusing at small test systems. In contrast, the potential of more recently proposed nonequilibrium methods, derived from the Jarzynski and Crooks equalities, has not yet fully been explored. Here, we compare the performance of these methods by calculating free energy differences for test systems at different levels of complexity and varying extent of the involved perturbations. We consider the interconversion of ethane into methanol, the switching of a tryptophane-sidechain in a tripeptide, and the binding of two different ligands to the globular protein snurportin 1. On the basis of our results, we suggest and assess a new nonequilibrium free energy method, Crooks Gaussian Intersection (CGI), which combines the advantages of existing methods. CGI is highly parallelizable and, for the test systems considered here, is shown to outperform the other studied equilibrium and nonequilibrium methods.     

A molecular dynamics computer simulation study of room-temperature ionic liquids. II. Equilibrium and nonequilibrium solvation dynamics

The molecular dynamics (MD) simulation study of solvation structure and free energetics in 1-ethyl-3-methylimidazolium chloride and 1-ethyl-3-methylimidazolium hexafluorophosphate using a probe solute in the preceding article [Y. Shim, M. Y. Choi and H. J. Kim, J. Chem. Phys. 122, 044510 (2005)] is extended to investigate dynamic properties of these liquids. Solvent fluctuation dynamics near equilibrium are studied via MD and associated time-dependent friction is analyzed via the generalized Langevin equation. Nonequilibrium solvent relaxation following an instantaneous change in the solute charge distribution and accompanying solvent structure reorganization are also investigated. Both equilibrium and nonequilibrium solvation dynamics are characterized by at least two vastly different time scales--a subpicosecond inertial regime followed by a slow diffusive regime. Solvent regions contributing to the subpicosecond nonequilibrium relaxation are found to vary significantly with initial solvation configurations, especially near the solute. If the solvent density near the solute is sufficiently high at the outset of the relaxation, subpicosecond dynamics are mainly governed by the motions of a few ions close to the solute. By contrast, in the case of a low local density, solvent ions located not only close to but also relatively far from the solute participate in the subpicosecond relaxation. Despite this difference, linear response holds reasonably well in both ionic liquids.     

Diffusive dynamics on multidimensional rough free energy surfaces

The dynamics of processes relevant to chemistry and biophysics on rough free energy landscapes is investigated using a recently developed algorithm to solve the Smoluchowski equation. Two different processes are considered: ligand rebinding in MbCO and protein folding. For the rebinding dynamics of carbon monoxide (CO) to native myoglobin (Mb) from locations around the active site, the two-dimensional free energy surface (FES) is constructed using extensive molecular dynamics simulations. The surface describes the minima in the A state (bound MbCO), CO in the distal pocket and in the Xe4 pocket, and the transitions between these states and allows to study the diffusion of CO in detail. For the folding dynamics of protein G, a previously determined two-dimensional FES was available. To follow the diffusive dynamics on these rough free energy surfaces, the Smoluchowski equation is solved using the recently developed hierarchical discrete approximation method. From the relaxation of the initial nonequilibrium distribution, experimentally accessible quantities such as the rebinding time for CO or the folding time for protein G can be calculated. It is found that the free energy barrier for CO in the Xe4 pocket and in the distal pocket (B state) closer to the heme iron is approximately 6 kcal/mol which is considerably larger than the inner barrier which separates the bound state and the B state. For the folding of protein G, a barrier of approximately 10 kcal/mol between the unfolded and the folded state is consistent with folding times of the order of milliseconds.     

Rosenbluth-sampled nonequilibrium work method for calculation of free energies in molecular simulation

We present methods that introduce concepts from Rosenbluth sampling [M. N. Rosenbluth and A. W. Rosenbluth, J. Chem. Phys. 23, 356 (1955)] into the Jarzynski nonequilibrium work (NEW) free-energy calculation technique [C. Jarzynski, Phys. Rev. Lett. 78, 2690 (1997)]. The proposed hybrid modifies the way steps are taken in the NEW process. With it, each step is selected from a range of alternatives, with bias given to steps that contribute the least work. The definition of the work average is modified to account for the bias. We introduce two variants of this method, lambda-bias sampling and configuration-bias sampling, respectively; a combined lambda- and configuration-bias method is also considered. By reducing the likelihood that large nonequilibrated work values enter the ensemble average, the Rosenbluth sampling aids in remedying problems of inaccuracy of the calculation. We demonstrate the performance of the proposed methods through a model system of N independent harmonic oscillators. This model captures the difficulties involved in calculating free energies in real systems while retaining many tractable features that are helpful to the study. We examine four variants of this model that differ qualitatively in the nature of their phase-space overlap. Results indicate that the lambda-bias sampling method is most useful for systems with entropic sampling barriers, while the configuration-bias methods are best for systems with energetic sampling barriers. The Rosenbluth-sampling schemes yield much more accurate results than the unbiased nonequilibrium work method. Typically the accuracy can be improved by about an order of magnitude for a given amount of sampling; this improvement translates into two or more orders of magnitude less sampling required to obtain a given level of accuracy, owing to the generally slow convergence of the NEW calculation when the inaccuracy is large.     

Continuous medium theory for nonequilibrium solvation: I. How to correctly evaluate solvation free energy of nonequilibrium

Considering the influences of electrostatic potential Phi upon the change of solute charge distribution deltarho and rho upon the change deltaPhi at the same time, a more reasonable integral formula of dG = (1/2) integral (V) (rhodeltaPhi + Phideltarho)dV is used to calculate the change of the electrostatic free energy in charging the solute-solvent system to a nonequilibrium state, instead of the one of dG = integral (V) PhideltarhodV used before. This modification improves the expressions of electrostatic free energy and solvation free energy, in which no quantity of the intermediate equilibrium state is explicitly involved. Detailed investigation reveals that the solvation free energy of nonequilibrium only contains the interaction energy between the field due to the solute charge in vacuum, and the dielectric polarization at the nonequilibrium state. The solvent reorganization energies of forward and backward electron transfer reactions have been redefined because the derivations lead to a remarkable feature that these quantities are direction-dependent, unlike the theoretical models developed before. The deductions are given in the electric field-displacement form. Relevant discussions on the reliability of theoretical models suggested in this work have also been presented.     

Linear free energy relations in solid state reactions II

Kinetics of the reaction between solid cobalt(II) acetate and solid substituted aniline hydrochlorides have been studied by thermal analysis. All the reactions go to completion. The crystalline products were characterized by physico-chemical methods. The energy of activation is 83.4,44.2,64.2,100.9,79.8,99.9, 135.9 and 168.7 kJ mol^–1for the reactions of aniline, 3- and 4-chloro-, 4-iodo-, 4-bromo-, 4-fluoro-, 4-methyl- and 4-methoxy aniline hydrochlorides with cobalt acetate. A plot of energy of activation against Hammett's constant is linear with a slope of –169.9 kJ mol^–1. A comparison between the reactivity of these reactions by the capillary and the thermal methods has been done.

---

Zusammenfassung Mittels Thermoanalyse wurde die Kinetik der Festkörperreaktion zwischen Kobalt(II)azetat und substituierten Anilinhydrochloriden untersucht. Alle Reaktionen laufen vollständig ab. Die kristallinen Produkte wurden mittels physikalisch-chemischer Methoden charakterisiert. Die Aktivierungsenergie der Reaktion von Kobalt(II)azetat mit Anilin, 3- und 4-Chlor-, 4-Jod-,4-Brom-, 4-Fluor-, 4-Methoxy- und 4-Methylanilin beträgt 83.4,44.2,64.2,100.9,79.8,99.9,168.7 bzw. 135.9 kJ/mol. Bei der Darstellung der Aktivierungsenergie in Abhängigkeit von der Hammettkonstante konnte Linearität mit einem Anstieg von – 169.9 kJ/mol festgestellt werden. Die Reaktivität dieser Reaktionen wurde mittels Kapillar- und thermischen Verfahrens verglichen.

---

. . - . , 3- 4--, 4--, 4--, 4--, 4- - 4- , , 83,4; 44,2; 64,2; 100,9; 79,8; 99,9; 135,9 168,7 ·^–1. - –169,9 ·^–1. , .

     

Variance minimization of free energy estimates from optimized expanded ensembles

We explored the possibility of improving the accuracy and precision of free-energy differences estimated via expanded ensembles by manipulation of the biasing weights. Three different weighing approaches were compared: the flat histogram (FH) method, the optimized ensemble (OE) method, and a method introduced in this work, denoted MinVar, which aims to explicitly minimize the expected variance. The performance of these three methods was tested for the simulation of chemical potentials in systems of symmetric diblock copolymers with chain lengths of either 10 or 4 beads, and a system of one large hard sphere of diameter 10 d immersed in a fluid of hard spheres of diameter d. In addition, the effect of the weighing method on the observed accuracy was investigated for different choices of macrostate staging and for both optimized and nonoptimized acceptance ratio methods for calculating free-energy differences. In the diblock copolymer systems, we found that the maximum attainable accuracy can be limited by correlations between the samples, causing the "real" observed variances to be much larger than the expected "ideal" ones. Hence, if the formal minimization of the variance, as aimed by the MinVar method, occurs at the expense of increasing the correlations in the data, the accuracy may actually decrease. Although maximizing the number of round trips between initial and final macrostates (as aimed by the OE method) was found to be directly related to data decorrelation, this only translates into increased accuracy if the correlations are the major source of errors in the free energy estimates. Finally, for the hard sphere system, we found that the MinVar method performs better than both the OE and FH methods even though the MinVar method in this case never completes a round trip, illustrating that maximizing the number of round trips for fixed computational cost does not necessarily lead to increased precision.     

Gamma-point lattice free energy estimates from O1 force calculations

We present a new method for estimating the vibrational free energy of crystal (and molecular) structures employing only a single force calculation, for a particularly displaced configuration, in addition to the calculation of the ground state configuration. This displacement vector is the sum of the phonon eigenvectors obtained from a fast-relative to, e.g., density-functional theory (DFT)-Hessian calculation using interatomic potentials. These potentials are based here on effective charges obtained from a DFT calculation of the ground state electronic charge density but could also be based on other, e.g., empiric approaches.