A hybrid recursion method to robustly ensure convergence efficiencies in the simulated scaling based free energy simulations

Recently, we developed an efficient free energy simulation technique, the simulated scaling (SS) method [H. Li et al., J. Chem. Phys. 126, 024106 (2007)], in the framework of generalized ensemble simulations. In the SS simulations, random walks in the scaling parameter space are realized so that both phase space overlap sampling and conformational space sampling can be simultaneously enhanced. To flatten the distribution in the scaling parameter space, in the original SS implementation, the Wang-Landau recursion was employed due to its well-known recursion capability. In the Wang-Landau recursion based SS free energy simulation scheme, at the early stage, recursion efficiencies are high and free energy regions are quickly located, although at this stage, the errors of estimated free energy values are large; at the later stage, the errors of estimated free energy values become smaller, however, recursions become increasingly slow and free energy refinements require very long simulation time. In order to robustly resolve this efficiency problem during free energy refinements, a hybrid recursion strategy is presented in this paper. Specifically, we let the Wang-Landau update method take care of the early stage recursion: the location of target free energy regions, and let the adaptive reweighting method take care of the late stage recursion: the refinements of free energy values. As comparably studied in the model systems, among three possible recursion procedures, the adaptive reweighting recursion approach is the least favorable one because of its low recursion efficiency during free energy region locations; and compared to the original Wang-Landau recursion approach, the proposed hybrid recursion technique can be more robust to guarantee free energy simulation efficiencies.     

Free energy calculation using molecular dynamics simulation combined with the three-dimensional reference interaction site model theory. II. Thermodynamic integration along a spatial reaction coordinate

We propose the thermodynamic integration along a spatial reaction coordinate using the molecular dynamics simulation combined with the three-dimensional reference interaction site model theory. This method provides a free energy calculation in solution along the reaction coordinate defined by the Cartesian coordinates of the solute atoms. The proposed method is based on the blue moon algorithm which can, in principle, handle any reaction coordinate as far as it is defined by the solute atom positions. In this article, we apply the present method to the complex formation process of the crown ether 18-Crown-6 (18C6) with the potassium ion in an aqueous solution. The separation between the geometric centers of these two molecules is taken to be the reaction coordinate for this system. The potential of mean force (PMF) becomes the maximum at the separation between the molecular centers being ～4 A?, which can be identified as the free energy barrier in the process of the molecular recognition. In a separation further than the free energy barrier, the PMF is slightly reduced to exhibit a plateau. In the region closer than the free energy barrier, approach of the potassium ion to the center of 18C6 also decreases the PMF. When the potassium ion is accommodated at the center of 18C6, the free energy is lower by -5.7 ± 0.7 kcal/mol than that at the above mentioned plateau or converged state. By comparing the results with those from the free energy calculation along the coupling parameters obtained in our previous paper [T. Miyata, Y. Ikuta, and F. Hirata, J. Chem. Phys. 133, 044114 (2010)], it is found that the effective interaction in water between 18C6 and the potassium ion vanishes beyond the molecular-center-separation of 10 A?. Furthermore, the conformation of 18C6 is found to be significantly changed depending upon the 18C6-K(+) distance. A proper conformational sampling and an accurate solvent treatment are crucial for realizing the accurate PMF, and we believe that the proposed method is useful to evaluate the PMF in a solution. A discussion upon the PMF in terms of the three-dimensional distribution function for the solvent is also presented.     

Solvent interaction energy calculations on molecular dynamics trajectories: increasing the efficiency using systematic frame selection

End-point methods such as linear interaction energy (LIE) analysis, molecular mechanics generalized Born solvent-accessible surface (MM/GBSA), and solvent interaction energy (SIE) analysis have become popular techniques to calculate the free energy associated with protein-ligand binding. Such methods typically use molecular dynamics (MD) simulations to generate an ensemble of protein structures that encompasses the bound and unbound states. The energy evaluation method (LIE, MM/GBSA, or SIE) is subsequently used to calculate the energy of each member of the ensemble, thus providing an estimate of the average free energy difference between the bound and unbound states. The workflow requiring both MD simulation and energy calculation for each frame and each trajectory proves to be computationally expensive. In an attempt to reduce the high computational cost associated with end-point methods, we study several methods by which frames may be intelligently selected from the MD simulation including clustering and address the question of how the number of selected frames influences the accuracy of the SIE calculations.     

流体混合物超额性质的分子模拟研究

提出了逐级取样分子模拟方法。通过逐个改变粒子的类别实现在一次模拟中对混合物整个组成范围的取样,并提出适当的处理方法可以同时获取混合物在整个组成范围内的多种超额热力学性质。用本文方法在零压、115.8K下对氩-氪体系作了模拟计算,与McDonald采用常规NPTMonteCarlo模拟结果的比较表明,本文的逐级取样模拟方法是可靠而有效的。     

Scalable free energy calculation of proteins via multiscale essential sampling

A multiscale simulation method, "multiscale essential sampling (MSES)," is proposed for calculating free energy surface of proteins in a sizable dimensional space with good scalability. In MSES, the configurational sampling of a full-dimensional model is enhanced by coupling with the accelerated dynamics of the essential degrees of freedom. Applying the Hamiltonian exchange method to MSES can remove the biasing potential from the coupling term, deriving the free energy surface of the essential degrees of freedom. The form of the coupling term ensures good scalability in the Hamiltonian exchange. As a test application, the free energy surface of the folding process of a miniprotein, chignolin, was calculated in the continuum solvent model. Results agreed with the free energy surface derived from the multicanonical simulation. Significantly improved scalability with the MSES method was clearly shown in the free energy calculation of chignolin in explicit solvent, which was achieved without increasing the number of replicas in the Hamiltonian exchange.     

Bimolecular electron transfers that follow a Sandros-Boltzmann dependence on free energy

Rate constants (k) for exergonic and endergonic electron-transfer reactions of equilibrating radical cations (A(?+) + B ? A + B(?+)) in acetonitrile could be fit well by a simple Sandros-Boltzmann (SB) function of the reaction free energy (ΔG) having a plateau with a limiting rate constant k(lim) in the exergonic region, followed, near the thermoneutral point, by a steep drop in log k vs ΔG with a slope of 1/RT. Similar behavior was observed for another charge shift reaction, the electron-transfer quenching of excited pyrylium cations (P(+)*) by neutral donors (P(+)* + D → P(?) + D(?+)). In this case, SB dependence was observed when the logarithm of the quenching constant (log k(q)) was plotted vs ΔG + s, where the shift term, s, equals +0.08 eV and ΔG is the free energy change for the net reaction (E(redox) - E(excit)). The shift term is attributed to partial desolvation of the radical cation in the product encounter pair (P(?)/D(?+)), which raises its free energy relative to the free species. Remarkably, electron-transfer quenching of neutral reactants (A* + D → A(?-) + D(?+)) using excited cyanoaromatic acceptors and aromatic hydrocarbon donors was also found to follow an SB dependence of log k(q) on ΔG, with a positive s, +0.06 eV. This positive shift contrasts with the long-accepted prediction of a negative value, -0.06 eV, for the free energy of an A(?-)/D(?+) encounter pair relative to the free radical ions. That prediction incorporated only a Coulombic stabilization of the A(?-)/D(?+) encounter pair relative to the free radical ions. In contrast, the results presented here show that the positive value of s indicates a decrease in solvent stabilization of the A(?-)/D(?+) encounter pair, which outweighs Coulombic stabilization in acetonitrile. These quenching reactions are proposed to proceed via rapidly interconverting encounter pairs with an exciplex as intermediate, A*/D ? exciplex ? A(?-)/D(?+). Weak exciplex fluorescence was observed in each case. For several reactions in the endergonic region, rate constants for the reversible formation and decay of the exciplexes were determined using time-correlated single-photon counting. The quenching constants derived from the transient kinetics agreed well with those from the conventional Stern-Volmer plots. For excited-state electron-transfer processes, caution is required in correlating quenching constants vs reaction free energies when ΔG exceeds ～+0.1 eV. Beyond this point, additional exciplex deactivation pathways-fluorescence, intersystem crossing, and nonradiative decay-are likely to dominate, resulting in a change in mechanism.     

Multiple free energies from a single simulation: extending enveloping distribution sampling to nonoverlapping phase-space distributions

A recently proposed method to obtain free energy differences for multiple end states from a single simulation of a reference state which was called enveloping distribution sampling (EDS) [J. Chem. Phys. 126, 184110 (2007)] is expanded to situations where the end state configuration space densities do not show overlap. It uses a reference state Hamiltonian suggested by Han in 1992 [Phys. Lett. A 165, 28 (1992)] in a molecular dynamics implementation. The method allows us to calculate multiple free energy differences "on the fly" from a single molecular dynamics simulation. The influence of the parameters on the accuracy and precision of the obtained free energy differences is investigated. A connection is established between the presented method and the Bennett acceptance ratio method. The method is applied to four two-state test systems (dipole inversion, van der Waals perturbation, charge inversion, and water to methanol conversion) and two multiple-state test systems [dipole inversion with five charging states and five (dis-)appearing water molecules]. Accurate results could be obtained for all test applications if the parameters of the reference state Hamiltonian were optimized according to a given algorithm. The deviations from the exact result or from an independent calculation were at most 0.6 kJ/mol. An accurate estimation of the free energy difference is always possible, independent of how different the end states are. However, the convergence times of the free energy differences are longer in cases where the end state configuration space densities do not show overlap [charge inversion, water to methanol conversion, (dis-)appearing water molecules] than in cases where the configuration space densities do show some overlap [(multiple) dipole inversion and van der Waals perturbation].     

Computational methods for the study of enzymic reaction mechanisms III: a perturbation plus QM/MM approach for calculating relative free energies of protonation

We describe a coupling parameter, that is, perturbation, approach to effectively create and annihilate atoms in the quantum mechanical Hamiltonian within the closed shell restricted Hartree-Fock formalism. This perturbed quantum mechanical atom (PQA) method is combined with molecular mechanics (MM) methods (PQA/MM) within a molecular dynamics simulation, to model the protein environment (MM region) effects that also make a contribution to the overall free energy change. Using the semiempirical PM3 method to model the QM region, the application of this PQA/MM method is illustrated by calculation of the relative protonation free energy of the conserved OD2 (Asp27) and the N5 (dihydrofolate) proton acceptor sites in the active site of Escherichia coli dihydrofolate reductase (DHFR) with the bound nicotinamide adenine dinucleotide phosphate (NADPH) cofactor. For a number of choices for the QM region, the relative protonation free energy was calculated as the sum of contributions from the QM region and the interaction between the QM and MM regions via the thermodynamic integration (TI) method. The results demonstrate the importance of including the whole substrate molecule in the QM region, and the overall protein (MM) environment in determining the relative stabilities of protonation sites in the enzyme active site. The PQA/MM free energies obtained by TI were also compared with those estimated by a less computationally demanding nonperturbative method based on the linear response approximation (LRA). For some choices of QM region, the total free energies calculated using the LRA method were in very close agreement with the PQA/MM values. However, the QM and QM/MM component free energies were found to differ significantly between the two methods.     

A theoretical investigation on the honeycomb potential fluid

A local self-consistent Ornstein-Zernike (OZ) integral equation theory (IET) is proposed to provide a rapid route for obtaining thermodynamic and structural information for any thermodynamically stable or metastable state points in the bulk phase diagram without recourse to traditional thermodynamic integration, and extensive NVT-Monte Carlo simulations are performed on a recently proposed honeycomb potential in three dimensions to test the theory's reliability. The simulated quantities include radial distribution function (rdf) and excess internal energy, pressure, excess chemical potential, and excess Helmholtz free energy. It is demonstrated that (i) the theory reproduces the rdf very satisfactorily only if the bulk state does not enter deep into a two phases coexistence region; (ii) the excess internal energy is the only one of the four thermodynamic quantities investigated amenable to the most accurate prediction by the present theory, and the simulated pressure is somewhat overestimated by the theoretical calculations, but the deviation tends to vanish along with rising of the temperature; (iii) using the structural functions from the present local self-consistent OZ IET, a previously derived local expression, due to the present author, achieves even a higher accuracy in calculating for the excess chemical potential than the exact virial pressure formula for the pressure, and the resulting excess Helmholtz free energy is in surprisingly same with the simulation results due to offset of the errors. Based on the above observations, it is suggested that it may be a good procedure to integrate the theoretical excess internal energy along the isochors to get the excess Helmholtz free energy, which is then fitted to a polynomial to be used for calculation of all of other thermodynamic quantities in the framework of the OZ IET.     

Improving the efficiency and reliability of free energy perturbation calculations using overlap sampling methods

A challenge in free energy calculation for complex molecular systems by computer simulation is to obtain a reliable estimate within feasible computational time. In this study, we suggest an answer to this challenge by exploring a simple method, overlap sampling (OS), for producing reliable free-energy results in an efficient way. The formalism of the OS method is based on ensuring sampling of important overlapping phase space during perturbation calculations. This technique samples both forward and reverse free energy perturbation (FEP) to improve the free-energy calculation. It considers the asymmetry of the FEP calculation and features an ability to optimize both the precision and the accuracy of the measurement without affecting the simulation process itself. The OS method is tested at two optimization levels: no optimization (simple OS), and full optimization (equivalent to Bennett's method), and compared to conventional FEP techniques, including the widely used direct FEP averaging method, on three alchemical mutation systems: (a) an anion transformation in water solution, (b) mutation between methanol and ethane, and (c) alchemical change of an adenosine molecule. It is consistently shown that the reliability of free-energy estimates can be greatly improved using the OS techniques at both optimization levels, while the performance of Bennett's method is particularly striking. In addition, the efficiency of a calculation can be significantly improved because the method is able to (a) converge to the right answer quickly, and (b) work for large perturbations. The basic two-stage OS method can be extended to admit additional stages, if needed. We suggest that the OS method can be used as a general perturbation technique for computing free energy differences in molecular simulations.     

Application of the constrained fluid lambda-integration path to the calculation of high temperature Au(110) surface free energies

Recently a method termed constrained fluid lambda-integration was proposed for calculating the free energy difference between bulk solid and liquid reference states via the construction of a reversible thermodynamic integration path; coupling the two states in question. The present work shows how the application of the constrained fluid lambda-integration concept to solid/liquid slab simulation cells makes possible a generally applicable computer simulation methodology for calculating the free energy of any surface and/or surface defect structure, including surfaces requiring variations in surface atom or density number, such as the (1 x 5) Au(100) or (1 x 2) missing row Au(110) reconstructed surfaces or excess adatom/vacancy/step populated surfaces. We evaluate the methodology by calculating the free energy of various disordered high temperature Au(110) embedded atom method surfaces constrained to differing excess surface atom numbers [including those corresponding to the (1 x 2) missing row reconstructed surface] and obtained the interesting result that at 1000 K (as distinct from lower temperatures) the free energy difference between these surfaces is reduced to zero; a result which is consistent with an expected order-disorder phase transition for the Au(110) surface at such high temperatures.     

Molecular dynamics study of micellization thermodynamics in AOT/hexane system

The thermodynamics of reverse micelle formation from an ionic surfactant, sodium bis(2-ethylhexyl) sulfosuccinate (Aerosol OT, AOT), in hexane is studied by molecular dynamics simulation. A change in the Gibbs free energy upon the addition of one AOT molecule to a reverse micelle is calculated as depending on aggregation number N by the thermodynamic integration method. This dependence has a minimum at N ≈ 20 and maximum at N ≈ 35 and predetermines the monotonically decreasing character of the standard chemical potential of AOT in a micelle with the increase of the aggregation number. The simulation results predict the formation of reverse AOT micelles with an average aggregation number of ≈30, which is in good agreement with experimental data.     

An application of coupled reference interaction site model/molecular dynamics to the conformational analysis of the alanine dipeptide

We present an application of our recently proposed coupled reference interaction site model (RISM) molecular dynamics (MD) solvation free energy methodology [Freedman and Truong, Chem. Phys. Lett. 381, 362 (2003); J. Chem. Phys. 121, 2187 (2004)] to study the conformational stability of alanine dipeptide in aqueous solution. In this methodology, radial distribution functions obtained from a single MD simulation are substituted into a RISM expression for solvation free energy. Consequently, iterative solution of the RISM equation is not needed. The relative solvation free energies of seven different conformations of the alanine dipeptide in aqueous solution are calculated. Results from the coupled RISM/MD methodology are in good agreement with those from earlier simulations using the accurate free energy perturbation approach, showing that the alphaR conformation is most stabilized by solution. This study establishes a framework for applying this coupled RISM/MD method to larger biological systems.     

Kinetics of Extremely Fast Ligand Exchange Reactions with Methylmercury(II)-complexes of 1-Methylpyridine-2-thione and 1-Methyl-quinaldine-4-thione: Rate-Equilibria Correlations

Methylmercury(II) transfer from two heterocyclic thiones to a wide variety of other ligands and the reverse reaction have been investigated by the temperature-jump method. The reactions are almost diffusion controlled, even when the free energy difference is negligible. The only exceptions are reactions to and from hydroxide. They are about hundred times slower than those with other ligands of comparable stability at zero free energy difference. This is in agreement with the behaviour of hydroxide in other methylmercury(II) exchange reactions. They follow an associative reaction mechanism. An empirical correlation of rate-and equilibrium-constants for atom transfer reactions is proposed. This equation shows some similarities to that of the Marcus type for atom transfer. However, the proposed equation not only describes reactions governed by an associative mechanism, but also those of dissociative type and cases between the two extremes.     

On the convergence improvement in the metadynamics simulations: a Wang-Landau recursion approach

As a popular tool in exploring free energy landscapes, the metadynamics method has been widely applied to elucidate various chemical or biochemical processes. As deeply discussed by Laio et al. [J. Phys. Chem. B 109, 6714 (2005)], the size of the updating Gaussian function is pivotal to the free energy convergence toward the target free energy surface. For instance, a greater Gaussian height can facilitate the quick visit of a conformation region of interest; however, it may lead to a larger error of the calculated free energy surface. In contrast, a lower Gaussian height can guarantee a better resolution of the calculated free energy surface; however, it will take longer time for such a simulation to navigate through the defined conformational region. In order to reconcile such confliction, the authors present a method by implementing the Wang-Landau recursion scheme in the metadynamics simulations to adaptively update the height of the unit Gaussian function. As demonstrated in their model studies on both a toy system, and a realistic molecular system treated with the hybrid quantum mechanical and molecular mechanical (QMMM) potential, the present approach can quickly result in more decently converged free energy surfaces, compared with the classical metadynamics simulations employing the fixed Gaussian heights.     

Application of a renormalization-group treatment to the statistical associating fluid theory for potentials of variable range (SAFT-VR)

An accurate prediction of phase behavior at conditions far and close to criticality cannot be accomplished by mean-field based theories that do not incorporate long-range density fluctuations. A treatment based on renormalization-group (RG) theory as developed by White and co-workers has proven to be very successful in improving the predictions of the critical region with different equations of state. The basis of the method is an iterative procedure to account for contributions to the free energy of density fluctuations of increasing wavelengths. The RG method has been combined with a number of versions of the statistical associating fluid theory (SAFT), by implementing White's earliest ideas with the improvements of Prausnitz and co-workers. Typically, this treatment involves two adjustable parameters: a cutoff wavelength L for density fluctuations and an average gradient of the wavelet function Φ. In this work, the SAFT-VR (variable range) equation of state is extended with a similar crossover treatment which, however, follows closely the most recent improvements introduced by White. The interpretation of White's latter developments allows us to establish a straightforward method which enables Φ to be evaluated; only the cutoff wavelength L then needs to be adjusted. The approach used here begins with an initial free energy incorporating only contributions from short-wavelength fluctuations, which are treated locally. The contribution from long-wavelength fluctuations is incorporated through an iterative procedure based on attractive interactions which incorporate the structure of the fluid following the ideas of perturbation theories and using a mapping that allows integration of the radial distribution function. Good agreement close and far from the critical region is obtained using a unique fitted parameter L that can be easily related to the range of the potential. In this way the thermodynamic properties of a square-well (SW) fluid are given by the same number of independent intermolecular model parameters as in the classical equation. Far from the critical region the approach provides the correct limiting behavior reducing to the classical equation (SAFT-VR). In the critical region the β critical exponent is calculated and is found to take values close to the universal value. In SAFT-VR the free energy of an associating chain fluid is obtained following the thermodynamic perturbation theory of Wertheim from the knowledge of the free energy and radial distribution function of a reference monomer fluid. By determining L for SW fluids of varying well width a unique equation of state is obtained for chain and associating systems without further adjustment of critical parameters. We use computer simulation data of the phase behavior of chain and associating SW fluids to test the accuracy of the new equation.     

Single-sweep methods for free energy calculations

A simple, efficient, and accurate method is proposed to map multidimensional free energy landscapes. The method combines the temperature-accelerated molecular dynamics (TAMD) proposed in [L. Maragliano and E. Vanden-Eijnden, Chem. Phys. Lett. 426, 168 (2006)] with a variational reconstruction method using radial-basis functions for the representation of the free energy. TAMD is used to rapidly sweep through the important regions of the free energy landscape and to compute the gradient of the free energy locally at points in these regions. The variational method is then used to reconstruct the free energy globally from the mean force at these points. The algorithmic aspects of the single-sweep method are explained in detail, and the method is tested on simple examples and used to compute the free energy of the solvated alanine dipeptide in two and four dihedral angles.     

Accurate and efficient generalized born model based on solvent accessibility: derivation and application for LogP octanol/water prediction and flexible peptide docking

A novel method for fast and accurate evaluation of the generalized Born radii in macromolecular solvation electrostatics calculations is proposed, based on the solvent accessibility of the first two solvation layers around an atom. The reverse generalized Born radii calculated by the method have correlation coefficient of 98.7% and RMSD of 0.031 A(-1) with the values obtained using a precise but significantly slower numerical boundary element solution. The method is applied to derive an estimate of the free solvation energy difference between octanol and water and to predict LogP octanol-water. A nine-parameter model is optimized on an 81 compound training set and applied to predict LogP(ow) for an external evaluation set of 19 drug molecules with RMSD of 0.9. The new GB approximation is also tested in Monte Carlo docking simulations of the fully flexible p53 peptide fragment to MDM2. The best energy solution found in the simulations has RMSD of 2.8 A to the X-ray structure.     

Construction of the free energy landscape of biomolecules via dihedral angle principal component analysis

A systematic approach to construct a low-dimensional free energy landscape from a classical molecular dynamics (MD) simulation is presented. The approach is based on the recently proposed dihedral angle principal component analysis (dPCA), which avoids artifacts due to the mixing of internal and overall motions in Cartesian coordinates and circumvents problems associated with the circularity of angular variables. Requiring that the energy landscape reproduces the correct number, energy, and location of the system's metastable states and barriers, the dimensionality of the free energy landscape (i.e., the number of essential components) is obtained. This dimensionality can be determined from the distribution and autocorrelation of the principal components. By performing an 800 ns MD simulation of the folding of hepta-alanine in explicit water and using geometric and kinetic clustering techniques, it is shown that a five-dimensional dPCA energy landscape is a suitable and accurate representation of the full-dimensional landscape. In the second step, the dPCA energy landscape can be employed (e.g., in a Langevin simulation) to facilitate a detailed investigation of biomolecular dynamics in low dimensions. Finally, several ways to visualize the multidimensional energy landscape are discussed.     

γ辐照 聚全氟乙丙烯(F46)的电子自旋共振研究 I.真空中俘获自由基的定量研究

A method for determining respectively the concentration of free radicials...CF2-CF2.(RI), ...CF2-CF-CF2...(RII)and...CF2- C(CF3)-CF2...(RIII) trapped in γ-irradiated F46 in vacuum at room temperature by analysing the ESR overlapping spectra is proposed. A quantitative relation between the concentration of free readicals and the total dose of γ-irradiation has been studied. The generating rate o free radicals caused by irradiation, i.e. G-value, has been calculated; it is 2.0 for RIII, 1.2 for RI, and 0.59 for RII. the magnitude of G value of free radicals is in the reverse order of their stability. Hence, we assume that the G value of free radicals in F46 mainly depends on the escaping rate of F2 and other fluorine-containing small molecules from surface of polymer into space. Therefore the looser the molecular chain, the larger becomes the G value of free radical. From that, we consider that both unstable free radicals RIII and RI with larger G-value mainly locate in the amorphous region, and the stable free radical RII with lower G-value is in the crystalline region. Recombinations between free radicals RII and RIII, or between RII and RI are forbidden at the temperature lower than 80C. An explanation for the saturation phenomenon of concentration of free radicals is given.