Multiple histogram reweighting method for the surface tension calculation

The multiple histogram reweighting method takes advantage of calculating ensemble averages over a range of thermodynamic conditions without performing a molecular simulation at each thermodynamic point. We show that this method can easily be extended to the calculation of the surface tension. We develop a new methodology called multiple histogram reweighting with slab decomposition based on the decomposition of the system into slabs along the direction normal to the interface. The surface tension is then calculated from local values of the chemical potential and of the configurational energy using Monte Carlo (MC) simulations. We show that this methodology gives surface tension values in excellent agreement with experiments and with standard NVT MC simulations in the case of the liquid-vapor interface of carbon dioxide.     

Optimized ensemble Monte Carlo simulations of dense Lennard-Jones fluids

We apply the recently developed adaptive ensemble optimization technique to simulate dense Lennard-Jones fluids and a particle-solvent model by broad-histogram Monte Carlo techniques. Equilibration of the simulated fluid is improved by sampling an optimized histogram in radial coordinates that shifts statistical weight towards the entropic barriers between the shells of the liquid. Interstitial states in the vicinity of these barriers are identified with unprecedented accuracy by sharp signatures in the quickly converging histogram and measurements of the local diffusivity. The radial distribution function and potential of mean force are calculated to high precision.     

A statistical analysis of the precision of reweighting-based simulations

Various advanced simulation techniques, which are used to sample the statistical ensemble of systems with complex Hamiltonians, such as those displayed in condensed matters and biomolecular systems, rely heavily on successfully reweighting the sampled configurations. The sampled points of a system from an elevated thermal environment or on a modified Hamiltonian are reused with different statistical weights to evaluate its properties at the initial desired temperature or of the original Hamiltonian. Often, the decrease of accuracy induced by this procedure is ignored and the final results can be far from what is expected. We have addressed the reasons behind such a phenomenon and have provided a quantitative method to estimate the number of sampled points required in the crucial step of reweighting of these advanced simulation methods. We also provided examples from temperature histogram reweighting and accelerated molecular dynamics reweighting to illustrate this idea, which can be generalized to the dynamic reweighting as well. The study shows that this analysis may provide a priori guidance for the strategy of setting up the parameters of advanced simulations before a lengthy one is carried out. The method can therefore provide insights for optimizing the parameters for high accuracy simulations with finite amount of computational resources.     

Free energy calculations,enhanced by a gaussian ansatz,for the “chemical work” distribution

The evaluation of the free energy is essential in molecular simulation because it is intimately related with the existence of multiphase equilibrium. Recently, it was demonstrated that it is possible to evaluate the Helmholtz free energy using a single statistical ensemble along an entire isotherm by accounting for the “chemical work” of transforming each molecule, from an interacting one, to an ideal gas. In this work, we show that it is possible to perform such a free energy perturbation over a liquid vapor phase transition. Furthermore, we investigate the link between a general free energy perturbation scheme and the novel nonequilibrium theories of Crook's and Jarzinsky. We find that for finite systems away from the thermodynamic limit the second law of thermodynamics will always be an inequality for isothermal free energy perturbations, resulting always to a dissipated work that may tend to zero only in the thermodynamic limit. The work, the heat, and the entropy produced during a thermodynamic free energy perturbation can be viewed in the context of the Crooks and Jarzinsky formalism, revealing that for a given value of the ensemble average of the “irreversible” work, the minimum entropy production corresponded to a Gaussian distribution for the histogram of the work. We propose the evaluation of the free energy difference in any free energy perturbation based scheme on the average irreversible “chemical work” minus the dissipated work that can be calculated from the variance of the distribution of the logarithm of the work histogram, within the Gaussian approximation. As a consequence, using the Gaussian ansatz for the distribution of the “chemical work,” accurate estimates for the chemical potential and the free energy of the system can be performed using much shorter simulations and avoiding the necessity of sampling the computational costly tails of the “chemical work.” For a more general free energy perturbation scheme that the Gaussian ansatz may not be valid, the free energy calculation can be expressed in terms of the moment generating function of the “chemical work” distribution. © 2014 Wiley Periodicals, Inc.     

方阱链状分子临界性质的Monte Carlo模拟

在巨正则系综下对阱宽为λ=1.5,链长分别为4、8、16的方阱链状流体实施Monte Carlo模拟,采用建立在完整标度基础上的无偏的Q-参数方法,通过histogram reweighting技术以及有限尺寸标度理论得到了热力学极限下该系列流体的临界温度和临界密度.模拟结果表明,方阱链流体的临界温度随着链长的增加而升高.并且不同链长方阱流体的临界温度均低于已报道的结果.由于本文所采用的完整标度的无偏性,我们估计的临界点更加准确.并且流体的临界温度与链长之间的关系与Flory-Huggins理论相一致.我们还预测了无限链长方阱流体的临界温度,比已有结果略高.     

Dynamical reweighting: improved estimates of dynamical properties from simulations at multiple temperatures

Dynamical averages based on functionals of dynamical trajectories, such as time-correlation functions, play an important role in determining kinetic or transport properties of matter. At temperatures of interest, the expectations of these quantities are often dominated by contributions from rare events, making the precise calculation of these quantities by molecular dynamics simulation difficult. Here, we present a reweighting method for combining simulations from multiple temperatures (or from simulated or parallel tempering simulations) to compute an optimal estimate of the dynamical properties at the temperature of interest without the need to invoke an approximate kinetic model (such as the Arrhenius law). Continuous and differentiable estimates of these expectations at any temperature in the sampled range can also be computed, along with an assessment of the associated statistical uncertainty. For rare events, aggregating data from multiple temperatures can produce an estimate with the desired precision at greatly reduced computational cost compared with simulations conducted at a single temperature. Here, we describe use of the method for the canonical (NVT) ensemble using four common models of dynamics (canonical distribution of Hamiltonian trajectories, Andersen thermostatting, Langevin, and overdamped Langevin or Brownian dynamics), but it can be applied to any thermodynamic ensemble provided the ratio of path probabilities at different temperatures can be computed. To illustrate the method, we compute a time-correlation function for solvated terminally-blocked alanine peptide across a range of temperatures using trajectories harvested using a modified parallel tempering protocol.     

A differential fluctuation theorem

We derive a nonequilibrium thermodynamics identity (the "differential fluctuation theorem") that connects forward and reverse joint probabilities of nonequilibrium work and of arbitrary generalized coordinates corresponding to states of interest. This identity allows us to estimate the free energy difference between domains of these states. Our results follow from a general symmetry relation between averages over nonequilibrium forward and backward path functions derived by Crooks [Crooks, G. E. Phys. Rev. E 2000, 61, 2361-2366]. We show how several existing nonequilibrium thermodynamic identities can be obtained directly from the differential fluctuation theorem. We devise an approach for measuring conformational free energy differences, and we demonstrate its applicability to the analysis of molecular dynamics simulations by estimating the free energy difference between two conformers of the alanine dipeptide model system. We anticipate that these developments can be applied to the analysis of laboratory experiments.     

Multiple Markov transition matrix method: Obtaining the stationary probability distribution from multiple simulations

We herein propose the multiple Markov transition matrix method (MMMM), an algorithm by which to estimate the stationary probability distribution from independent multiple molecular dynamics simulations with different Hamiltonians. Applications to the potential of mean force calculation in combination with the umbrella sampling method are presented. First, the performance of the MMMM is examined in the case of butane. Compared with the weighted histogram analysis method (WHAM), the MMMM has an advantage with respect to the reasonable evaluation of the stationary probability distribution even from nonequilibrium trajectories. This method is then applied to Met‐enkephalin nonequilibrium simulation. © 2008 Wiley Periodicals, Inc. J Comput Chem, 2009     

Probability distribution analysis of single-molecule fluorescence anisotropy and resonance energy transfer

Analysis of anisotropy in single-molecule fluorescence experiments using the probability distribution analysis (PDA) method is presented. The theory of anisotropy-PDA is an extension of the PDA theory recently developed for the analysis of F?rster resonance energy transfer (FRET) signals [Antonik, M.; et al. J. Phys. Chem. B 2006, 110, 6970]. The PDA method predicts the shape of anisotropy histograms for any given expected ensemble anisotropy, signal intensity distribution, and background. Further improvements of the PDA theory allow one to work with very low photon numbers, i.e., starting from the level of background signal. Analysis of experimental and simulated data shows that PDA has the major advantage to unambiguously distinguish between shot noise broadening and broadening caused by heterogeneities in the sample. Fitting of experimental histograms yields anisotropy values of individual species, which can be directly compared with those measured in ensemble experiments. Excellent agreement between the ensemble data and the results of PDA demonstrates a good absolute accuracy of the PDA method. The precision in determination of mean values depends mainly on the total number of photons, whereas the ability of PDA to detect the presence of heterogeneities strongly depends on the time window length. In its present form PDA can be also applied to computed fluorescence parameters such as FRET efficiency and scatter-corrected fluorescence anisotropy. Extension of the PDA theory to low photon numbers makes it possible to apply PDA to dynamic systems, for which high time resolution is required. In this way PDA is developed as a sensitive tool to detect biomolecular heterogeneities in space and time.     

The nature of the calculation of the pressure in molecular simulations of continuous models from volume perturbations

We consider some fundamental aspects of the calculation of the pressure from simulations by performing volume perturbations. The method, initially proposed for hard-core potentials by Eppenga and Frenkel [Mol. Phys.52, 1303 (1984)] and then extended to continuous potentials by Harismiadis et al. [J. Chem. Phys. 105, 8469 (1996)], is based on the numerical estimate of the change in Helmholtz free energy associated with the perturbation which, in turn, can be expressed as an ensemble average of the corresponding Boltzmann factor. The approach can be easily generalized to the calculation of components of the pressure tensor and also to ensembles other than the canonical ensemble. The accuracy of the method is assessed by comparing simulation results obtained from the volume-perturbation route with those obtained from the usual virial expression for several prototype fluid models. Monte Carlo simulation data are reported for bulk fluids and for inhomogeneous systems containing a vapor-liquid interface.     

Single molecule photobehavior of a chromophore interacting with silica-based nanomaterials

Single molecule studies of the free DY-630-MI and interacting with MCM-41 and (Al)MCM-41, show the conformational diversity of the molecule. The free dye is characterized by a single broad (fwhm = 0.7 ns) lifetime distribution histogram centered on 1.47 ns, which is also reflected in the broadness of the polarization value distribution histogram, covering almost the full range of values from -1 to 1. The fluorescence intensity traces of the free DY-630-MI show strong blinking behavior and weak photostability. Upon interaction with the mesoporous silica nanomaterials, MCM-41 and (Al)MCM-41, the dye molecule becomes more stable, with less blinking present in the fluorescence traces. The lifetime distribution histogram in the case of DY-630-MI/MCM-41 complexes is fitted by 3 Gaussians, indicating 3 distinct interaction sites. The Gaussian with the largest amplitude is centered on 2.19 ns, consistent with the confinement effect of MCM-41 and in agreement with the ensemble average studies. The polarization value distribution histogram becomes narrower in comparison with the free molecule and is more biased towards the positive limit. Replacing few Si(4+) ions with Al(3+) ones in the regular MCM-41 changes the local electrostatic field within the nanotube. This atomic substitution in the nanohosts results in a more selective orientation of the dye molecules, giving two populations with time constants 1.56 and 2.10 ns.     

Development of the theory of a self-consistent field for polyelectrolyte solutions

A new theoretical approach to calculating the thermodynamic and structure functions of polyelectrolyte solutions is proposed, based on the method of Gaussian equivalent representation for calculating the functional integrals. Formulas for the mean-force potential, osmotic pressure, and complete monomer-monomer pair distribution functions are presented. A sodium polystyrene sulfonate solution with NaCl additives is considered as an example.     

Geometrical cluster ensemble analysis of random sphere packings

We introduce a geometric analysis of random sphere packings based on the ensemble averaging of hard-sphere clusters generated via local rules including a nonoverlap constraint for hard spheres. Our cluster ensemble analysis matches well with computer simulations and experimental data on random hard-sphere packing with respect to volume fractions and radial distribution functions. To model loose as well as dense sphere packings various ensemble averages are investigated, obtained by varying the generation rules for clusters. Essential findings are a lower bound on volume fraction for random loose packing that is surprisingly close to the freezing volume fraction for hard spheres and, for random close packing, the observation of an unexpected split peak in the distribution of volume fractions for the local configurations. Our ensemble analysis highlights the importance of collective and global effects in random sphere packings by comparing clusters generated via local rules to random sphere packings and clusters that include collective effects.     

Multiple Histogram Method and Static Monte Carlo Sampling

Summary: We describe an approach to use multiple‐histogram methods in combination with static, biased Monte Carlo simulations. To illustrate this, we computed the force‐extension curve of an athermal polymer from multiple histograms constructed in a series of static Rosenbluth Monte Carlo simulations. From the complete histogram of the distribution function of the end‐to‐end vectors of the polymer chain, we can efficiently compute the polymer force‐extension curve.

Comparison of the stress‐strain curves for the stress ensemble (symbols) and the strain ensemble (lines). Results obtained for N = 100, 200, 400, and 600. For small x, f(x) = −F′(x) was computed by aproximating F(x) by a second degree polynomial and then taking the derivative. For large x, f(x) = −F′(x) was computed numerically.     

Bayesian analysis of heterogeneity in the distribution of binding properties of immobilized surface sites

Once a homogeneous ensemble of a protein ligand is taken from solution and immobilized to a surface, for many reasons the resulting ensemble of surface binding sites to soluble analytes may be heterogeneous. For example, this can be due to the intrinsic surface roughness causing variations in the local microenvironment, nonuniform density distribution of polymeric linkers, or nonuniform chemical attachment producing different protein orientations and conformations. We previously described a computational method for determining the distribution of affinity and rate constants of surface sites from analysis of experimental surface binding data. It fully exploits the high signal/noise ratio and reproducibility provided by optical biosensor technology, such as surface plasmon resonance. Since the computational analysis is ill conditioned, the previous approach used a regularization strategy assuming a priori all binding parameters to be equally likely, resulting in the broadest possible parameter distribution consistent with the experimental data. We now extended this method in a Bayesian approach to incorporate the opposite assumption, i.e., that the surface sites a priori are expected to be uniform (as one would expect in free solution). This results in a distribution of binding parameters as close to monodispersity as possible given the experimental data. Using several model protein systems immobilized on a carboxymethyl dextran surface and probed with surface plasmon resonance, we show microheterogeneity of the surface sites in addition to broad populations of significantly altered affinity. The distributions obtained are highly reproducible. Immobilization conditions and the total surface density of immobilized sites can have a substantial impact on the functional distribution of the binding sites.     

Theory of binless multi-state free energy estimation with applications to protein-ligand binding

The weighted histogram analysis method (WHAM) is routinely used for computing free energies and expectations from multiple ensembles. Existing derivations of WHAM require observations to be discretized into a finite number of bins. Yet, WHAM formulas seem to hold even if the bin sizes are made arbitrarily small. The purpose of this article is to demonstrate both the validity and value of the multi-state Bennet acceptance ratio (MBAR) method seen as a binless extension of WHAM. We discuss two statistical arguments to derive the MBAR equations, in parallel to the self-consistency and maximum likelihood derivations already known for WHAM. We show that the binless method, like WHAM, can be used not only to estimate free energies and equilibrium expectations, but also to estimate equilibrium distributions. We also provide a number of useful results from the statistical literature, including the determination of MBAR estimators by minimization of a convex function. This leads to an approach to the computation of MBAR free energies by optimization algorithms, which can be more effective than existing algorithms. The advantages of MBAR are illustrated numerically for the calculation of absolute protein-ligand binding free energies by alchemical transformations with and without soft-core potentials. We show that binless statistical analysis can accurately treat sparsely distributed interaction energy samples as obtained from unmodified interaction potentials that cannot be properly analyzed using standard binning methods. This suggests that binless multi-state analysis of binding free energy simulations with unmodified potentials offers a straightforward alternative to the use of soft-core potentials for these alchemical transformations.     

Thermodynamics and partitioning of homopolymers into a slit-A grand canonical Monte Carlo simulation study

Grand canonical ensemble Monte Carlo simulation (GCMC) combined with the histogram reweighting technique was used to study the thermodynamic equilibrium of a homopolymer solution between a bulk and a slit pore. GCMC gives the partition coefficients that agree with those from canonical ensemble Monte Carlo simulations in a twin box, and it also gives results that are not accessible through the regular canonical ensemble simulation such as the osmotic pressure of the solution. In a bulk polymer solution, the calculated osmotic pressure agrees very well with the scaling theory predictions both for the athermal polymer solution and the theta solution. However, one cannot obtain the osmotic pressure of the confined solution in the same way since the osmotic pressure of the confined solution is anisotropic. The chemical potentials in GCMC simulations were found to differ by a translational term from the chemical potentials obtained from canonical ensemble Monte Carlo simulations with the chain insertion method. This confirms the equilibrium condition of a polymer solution partition between the bulk and a slit pore: the chemical potentials of the polymer chain including the translational term are equal at equilibrium. The histogram reweighting method enables us to obtain the partition coefficients in the whole range of concentrations based on a limited set of simulations. Those predicted bulk-pore partition coefficient data enable us to perform further theoretical analysis. Scaling predictions of the partition coefficient at different regimes were given and were confirmed by the simulation data.     

A distribution-based method to resolve single-molecule Förster resonance energy transfer observations

We introduce a new approach to analyze single-molecule Fo?rster resonance energy transfer (FRET) data. The method recognizes that FRET efficiencies assumed by traditional ensemble methods are unobservable for single molecules. We propose instead a method to predict distributions of FRET parameters obtained directly from the data. Distributions of FRET rates, given the data, are precisely defined using Bayesian methods and increase the information derived from the data. Benchmark comparisons find that the response time of the new method outperforms traditional methods of averaging. Our approach makes no assumption about the number or distribution of underlying FRET states. The new method also yields information about joint parameter distributions going beyond the standard framework of FRET analysis. For example, the running distribution of FRET means contains more information than any conceivable single measure of FRET efficiency. The method is tested against simulated data and then applied to a pilot-study sample of calmodulin molecules immobilized in lipid vesicles, revealing evidence for multiple dynamical states.