Determining the density of states for classical statistical models: a random walk algorithm to produce a flat histogram

A water-explicit lattice model of heat-, cold-, and pressure-induced protein unfolding

We investigate the effect of temperature and pressure on polypeptide conformational stability using a two-dimensional square lattice model in which water is represented explicitly. The model captures many aspects of water thermodynamics, including the existence of density anomalies, and we consider here the simplest representation of a protein: a hydrophobic homopolymer. We show that an explicit treatment of hydrophobic hydration is sufficient to produce cold, pressure, and thermal denaturation. We investigate the effects of the enthalpic and entropic components of the water-protein interactions on the overall folding phase diagram, and show that even a schematic model such as the one we consider yields reasonable values for the temperature and pressure ranges within which highly compact homopolymer configurations are thermodynamically stable.

Effective sampling in the configurational space of a small peptide by the multicanonical-multioverlap algorithm

We propose a generalized-ensemble algorithm, which we refer to as the multicanonical-multioverlap algorithm. By utilizing a non-Boltzmann weight factor, this method realizes a random walk in the multidimensional, energy-overlap space and explores widely in the configurational space including specific configurations, where the overlap of a configuration with respect to a reference state is a measure for structural similarity. We apply the multicanonical-multioverlap molecular dynamics method to a penta peptide, Met-enkephalin, in vacuum as a test system. We also apply the multicanonical and multioverlap molecular dynamics methods to this system for the purpose of comparisons. We see that the multicanonical-multioverlap molecular dynamics method realizes effective sampling in the configurational space including specific configurations more than the other two methods. Furthermore, from the results of the multicanonical-multioverlap molecular dynamics simulation, we obtain a local-minimum state of the Met-enkephalin system.

Proposed high-pressure calorimetric experiment to probe theoretical predictions on the liquid-liquid critical point hypothesis

A high-pressure calorimetric experiment to test the validity of a water model's prediction, regarding the behavior of the heat capacity in the liquid-liquid critical point and singularity-free scenarios, is proposed in this paper. The response of the specific heat at high pressure is different depending on the existence, or not, of a second critical point. If the model presents a second critical point, there is a nondivergent maximum in the specific heat at a temperature T approximately T(c) for any pressure P > P(c) (T(c) and P(c) being the temperature and pressure of the model's second critical point). This maximum does not appear for the model with a singularity free scenario.

Accuracy and convergence of the Wang-Landau sampling algorithm

We present estimations of the accuracy and convergence of the Wang-Landau algorithm. Both accuracy and the related length of the Monte Carlo run depend on the modification parameter f and the density of states. The analytical solution obtained for the two-level system was checked numerically on the two-dimensional Ising model. Although the two-level system is a very simple one, it appears that the proposed solution describes the generic features of the Wang-Landau algorithm. The estimations should prove useful in Monte Carlo calculations of protein folding, first-order transitions, and other systems with a rough energy landscape.

Metabasin Approach for Computing the Master Equation Dynamics of Systems with Broken Ergodicity

We propose a technique for computing the master equation dynamics of systems with broken ergodicity. The technique involves a partitioning of the system into components, or metabasins, where the relaxation times within a metabasin are short compared to an observation time scale. In this manner, equilibrium statistical mechanics is assumed within each metabasin, and the intermetabasin dynamics are computed using a reduced set of master equations. The number of metabasins depends upon both the temperature of the system and its derivative with respect to time. With this technique, the integration time step of the master equations is governed by the observation time scale rather than the fastest transition time between basins. We illustrate the technique using a simple model landscape with seven basins and show validation against direct Euler integration. Finally, we demonstrate the use of the technique for a realistic glass-forming system (viz., selenium) where direct Euler integration is not computationally feasible.

Toward absolute density of states calculations for proteins

The density of states (DOS), which gives the number of conformations with a particular energy E, is a prerequisite in computing most thermodynamic quantities and in elucidating important biological processes such as the mechanism of protein folding. However, current methods for computing DOS of large systems such as proteins generally yield only the ratios of microstate counts for different energies, which could yield absolute conformation counts if the total number of conformations in phase space is known, thus motivating this work. Here, the total number of energy minima of 50-mer polyalanine, whose size corresponds to naturally occurring small proteins, was estimated under an all-atom potential energy function based on the cumulative distribution function (CDF) of conformational differences to be approximately 10(38). This estimate can place any DOS function, such as the Gaussian DOS distribution in the random energy model, on an absolute scale. Comparing the absolute conformational counts from a Gaussian DOS function with those from the CDF derived from quenched molecular dynamics ensembles shows that the former are far greater than the latter, indicating far fewer low-energy minima actually exist. In addition to showing how CDF and relative DOS calculations can yield absolute DOS for a discrete system, we also show how they can yield absolute DOS for continuous variable systems to a specified atomic variance. In the context of protein folding, knowing this phase-space "volume" of conformations in a DOS function, as well as characteristic transition times, constrains the set of possible folding mechanisms.

Structure optimization and folding mechanisms of off-lattice protein models using statistical temperature molecular dynamics simulation: Statistical temperature annealing

The recently proposed statistical temperature molecular dynamics (STMD) algorithm [Kim, Phys. Rev. Lett. 97, 050601 (2006)] is used as the core of an optimization algorithm, statistical temperature annealing (STA), for finding low-lying energy minima of complex potential energy landscapes. Since STMD realizes a random walk in energy, the idea is simply to initiate repeated minimizations from configurations in the low-energy segments of STMD trajectories. STA is tested in structural optimization of various off-lattice AB and extended AB protein models in two and three dimensions with different chain lengths. New putative ground states were found for the two- and three-dimensional AB 55-mer, and for the three-dimensional extended AB 21-mer and 55-mer. The distinct folding features of the models are analyzed in terms of the statistical temperature and other representations of the structure of the potential energy landscape. It is shown that the characteristic behavior of the statistical temperature undergoes a qualitative change with the inclusion of a torsional potential in the extended AB model, as the more rigid backbone makes the potential energy landscape more funnel-like.

Phase behavior of attractive and repulsive ramp fluids: Integral equation and computer simulation studies

Using computer simulations and a thermodynamically self-consistent integral equation we investigate the phase behavior and thermodynamic anomalies of a fluid composed of spherical particles interacting via a two-scale ramp potential (a hard core plus a repulsive and an attractive ramp) and the corresponding purely repulsive model. Both simulation and integral equation results predict a liquid-liquid demixing when attractive forces are present, in addition to a gas-liquid transition. Furthermore, a fluid-solid transition emerges in the neighborhood of the liquid-liquid transition region, leading to a phase diagram with a somewhat complicated topology. This solidification at moderate densities is also present in the repulsive ramp fluid, but in this case inhibits the fluid-fluid separation.

Unexpectedly normal phase behavior of single homopolymer chains

Employing Monte Carlo simulations, we show that the topology of the phase diagram of a single flexible homopolymer chain changes in dependence on the range of an attractive square well interaction between the monomers. For a range of attraction larger than a critical value, the equilibrium phase diagram of the single polymer chain and the corresponding polymer solution phase diagram exhibit vapor (swollen coil, dilute solution), liquid (collapsed globule, dense solution), and solid phases. Otherwise, the liquid-vapor transition vanishes from the equilibrium phase diagram for both the single chain and the polymer solution. This change in topology of the phase diagram resembles the behavior known for colloidal dispersions. The interplay of enthalpy and conformational entropy in the polymer case thus can lead to the same topology of phase diagrams as the interplay of enthalpy and translational entropy in simple liquids.

Critical endpoint behavior in an asymmetric Ising model: application of Wang-Landau sampling to calculate the density of states

Using the Wang-Landau sampling method with a two-dimensional random walk we determine the density of states for an asymmetric Ising model with two- and three-body interactions on a triangular lattice, in the presence of an external field. With an accurate density of states we were able to map out the phase diagram accurately and perform quantitative finite-size analyses at, and away from, the critical endpoint. We observe a clear divergence of the curvature of the spectator phase boundary and of the magnetization coexistence diameter derivative at the critical endpoint, and the exponents for both divergences agree well with previous theoretical predictions.

Wang-Landau molecular dynamics technique to search for low-energy conformational space of proteins

Multicanonical molecular dynamics (MD) is a powerful technique for sampling conformations on rugged potential surfaces such as protein. However, it is notoriously difficult to estimate the multicanonical temperature effectively. Wang and Landau developed a convenient method for estimating the density of states based on a multicanonical Monte Carlo method. In their method, the density of states is calculated autonomously during a simulation. In this paper, we develop a set of techniques to effectively apply the Wang-Landau method to MD simulations. In the multicanonical MD, the estimation of the derivative of the density of states is critical. In order to estimate it accurately, we devise two original improvements. First, the correction for the density of states is made smooth by using the Gaussian distribution obtained by a short canonical simulation. Second, an approximation is applied to the derivative, which is based on the Gaussian distribution and the multiple weighted histogram technique. A test of this method was performed with small polypeptides, Met-enkephalin and Trp-cage, and it is demonstrated that Wang-Landau MD is consistent with replica exchange MD but can sample much larger conformational space.

Monte Carlo method for computing density of states and quench probability of potential energy and enthalpy landscapes

The thermodynamics and kinetics of a many-body system can be described in terms of a potential energy landscape in multidimensional configuration space. The partition function of such a landscape can be written in terms of a density of states, which can be computed using a variety of Monte Carlo techniques. In this paper, a new self-consistent Monte Carlo method for computing density of states is described that uses importance sampling and a multiplicative update factor to achieve rapid convergence. The technique is then applied to compute the equilibrium quench probability of the various inherent structures (minima) in the landscape. The quench probability depends on both the potential energy of the inherent structure and the volume of its corresponding basin in configuration space. Finally, the methodology is extended to the isothermal-isobaric ensemble in order to compute inherent structure quench probabilities in an enthalpy landscape.

Sampling enhancement for the quantum mechanical potential based molecular dynamics simulations: a general algorithm and its extension for free energy calculation on rugged energy surface

An approach is developed in the replica exchange framework to enhance conformational sampling for the quantum mechanical (QM) potential based molecular dynamics simulations. Importantly, with our enhanced sampling treatment, a decent convergence for electronic structure self-consistent-field calculation is robustly guaranteed, which is made possible in our replica exchange design by avoiding direct structure exchanges between the QM-related replicas and the activated (scaled by low scaling parameters or treated with high "effective temperatures") molecular mechanical (MM) replicas. Although the present approach represents one of the early efforts in the enhanced sampling developments specifically for quantum mechanical potentials, the QM-based simulations treated with the present technique can possess the similar sampling efficiency to the MM based simulations treated with the Hamiltonian replica exchange method (HREM). In the present paper, by combining this sampling method with one of our recent developments (the dual-topology alchemical HREM approach), we also introduce a method for the sampling enhanced QM-based free energy calculations.

Ring polymer simulations with global radius of curvature

We simulate three-dimensional flexible off-lattice ring polymers of length L up to L=4000 for various values of the global radius of curvature Rgrc=0.25 , 0.48, and 1.0 and Rgrc=2.0 . We utilize two different ensembles: one with a delta -function constraint on the radius, and the other with a theta -function. For both cases the global radius of curvature provides a valid regularization of polymers with thickness D=2Rgrc . The Flory-type critical exponent nu SAW of self-avoiding rings at D=2 is found to be nu SAW=0.5869(5) from the radii of gyration chain length scaling, while other D values produce consistent results. For our current implementation, the numerical effort of chain thickness calculations is bounded by a number O(LlnL) per single update. We also study low-temperature configurations of spatially dense Lennard-Jones homopolymers on a ring and identify some conformational building blocks.

Computation of free energy profiles with parallel adaptive dynamics

We propose a formulation of an adaptive computation of free energy differences, in the adaptive biasing force or nonequilibrium metadynamics spirit, using conditional distributions of samples of configurations which evolve in time. This allows us to present a truly unifying framework for these methods, and to prove convergence results for certain classes of algorithms. From a numerical viewpoint, a parallel implementation of these methods is very natural, the replicas interacting through the reconstructed free energy. We demonstrate how to improve this parallel implementation by resorting to some selection mechanism on the replicas. This is illustrated by computations on a model system of conformational changes.

Statistical temperature molecular dynamics: Application to coarse-grained β-barrel-forming protein models

Recently the authors proposed a novel sampling algorithm, "statistical temperature molecular dynamics" (STMD) [J. Kim et al., Phys. Rev. Lett. 97, 050601 (2006)], which combines ingredients of multicanonical molecular dynamics and Wang-Landau sampling. Exploiting the relation between the statistical temperature and the density of states, STMD generates a flat energy distribution and efficient sampling with a dynamic update of the statistical temperature, transforming an initial constant estimate to the true statistical temperature T(U), with U being the potential energy. Here, the performance of STMD is examined in the Lennard-Jones fluid with diverse simulation conditions, and in the coarse-grained, off-lattice BLN 46-mer and 69-mer protein models, exhibiting rugged potential energy landscapes with a high degree of frustration. STMD simulations combined with inherent structure (IS) analysis allow an accurate determination of protein thermodynamics down to very low temperatures, overcoming quasiergodicity, and illuminate the transitions occurring in folding in terms of the energy landscape. It is found that a thermodynamic signature of folding is significantly suppressed by accurate sampling, due to an incoherent contribution from low-lying non-native IS in multifunneled landscapes. It is also shown that preferred accessibility to such IS during the collapse transition is intimately related to misfolding or poor foldability.

Density of states of a two-dimensional XY model from the Wang-Landau algorithm

Using the Wang-Landau algorithm combined with the analytic method, the density of states of a two-dimensional XY model on square lattices of sizes 16x16, 24x24, and 32x32 is accurately calculated. Thermodynamic quantities, such as internal energy, free energy, entropy, and specific heat are obtained from the resulted density of states by numerical integration. It is found that both the density of states and the extensive thermodynamic quantities obey the scaling law. From the entropy curve, symptoms of phase transition are observed. A general method of calculation of the density of states of continuous models by simulation combined with analytical method is proposed.

Fast algorithm to calculate density of states

An algorithm to calculate the density of states, based on the well-known Wang-Landau method, is introduced. Independent random walks are performed in different restricted ranges of energy, and the resultant density of states is modified by a function of time, F(t) proportional to t-1, for large time. As a consequence, the calculated density of state, gm(E,t) , approaches asymptotically the exact value g(ex)(E) as proportional to t-1/2, avoiding the saturation of the error. It is also shown that the growth of the interface of the energy histogram belongs to the random deposition universality class.

Microcanonical approach to the simulation of first-order phase transitions

A simple microcanonical strategy for the simulation of first-order phase transitions is proposed. At variance with flat-histogram methods, there is no iterative parameters optimization nor long waits for tunneling between the ordered and the disordered phases. We test the method in the standard benchmark: the Q-states Potts model (Q=10 in two dimensions and Q=4 in D=3). We develop a cluster algorithm for this model, obtaining accurate results for systems with more than 10(6) spins.

Phase-space overlap measures. I. Fail-safe bias detection in free energies calculated by molecular simulation

We consider ways to quantify the overlap of the parts of phase space important to two systems, labeled A and B. Of interest is how much of the A-important phase space lies in that important to B, and how much of B lies in A. Two measures are proposed. The first considers four total-energy distributions, formed from all combinations made by tabulating either the A-system or the B-system energy when sampling either the A or B system. Measures for A in B and B in A are given by two overlap integrals defined on pairs of these distributions. The second measure is based on information theory, and defines two relative entropies which are conveniently expressed in terms of the dissipated work for free-energy perturbation (FEP) calculations in the A-->B and B-->A directions, respectively. Phase-space overlap is an important consideration in the performance of free-energy calculations. To demonstrate this connection, we examine bias in FEP calculations applied to a system of independent particles in a harmonic potential. Systems are selected to represent a range of overlap situations, including extreme subset, subset, partial overlap, and nonoverlap. The magnitude and symmetry of the bias (A-->B vs B-->A) are shown to correlate well with the overlap, and consequently with the overlap measures. The relative entropies are used to scale the amount of sampling to obtain a universal bias curve. This result leads to develop a simple heuristic that can be applied to determine whether a work-based free-energy measurement is free of bias. The heuristic is based in part on the measured free energy, but we argue that it is fail-safe inasmuch as any bias in the measurement will not promote a false indication of accuracy.

Free energy landscapes of electron transfer system in dipolar environment below and above the rotational freezing temperature

Electron transfer reaction in a polar solvent is modeled by a solute dipole surrounded by dipolar molecules with simple rotational dynamics posted on the three-dimensional distorted lattice sites. The interaction energy between the solute and solvent dipoles as a reaction coordinate is adopted and free energy landscapes are calculated by generating all possible states for a 26 dipolar system and by employing Wang-Landau sampling algorithm for a 92 dipolar system. For temperatures higher than the energy scale of dipole-dipole interactions, the free energy landscapes for the small reaction coordinate region have quadratic shape as predicted by Marcus [Rev. Mod. Phys. 65, 599 (1993)] whereas for the large reaction coordinate region, the landscapes exhibit a nonquadratic shape. When the temperature drops, small notched structures appear on the free energy profiles because of the frustrated interactions among dipoles. The formation of notched structure is analyzed with statistical approach and it is shown that the amplitude of notched structure depend upon the segment size of the reaction coordinate and is characterized by the interaction energy among the dipoles. Using simulated free energy landscapes, the authors calculate the reaction rates as a function of the energy gap for various temperatures. At high temperature, the reactions rates follow a bell shaped (inverted parabolic) energy gap law in the small energy gap regions, while it becomes steeper than the parabolic shape in a large energy gap regions due to the nonquadratic shape of the free energy landscape. The peak position of parabola also changes as the function of temperature. At low temperature, the profile of the reaction rates is no longer smooth because of the many local minima of the free energy landscape.

Dense orientationally ordered states of a single semiflexible macromolecule: an expanded ensemble Monte Carlo simulation

Using a coarse-grained model we perform a Monte Carlo simulation of the state behavior of an individual semiflexible macromolecule. Chains consisting of N = 256 and 512 monomer units have been investigated. A recently proposed enhanced sampling Monte Carlo technique for the bond fluctuation model in an expanded ensemble in four-dimensional coordinate space was applied. The algorithm allows one to accelerate the sampling of statistically independent three-dimensional conformations in a dense globular state. We found that the temperature of the intraglobular liquid-solid transition decreases with increasing chain stiffness. We have investigated the possible intraglobular orientationally ordered (i.e., liquid-crystalline) structures and obtained a diagram of states for chains consisting of N = 256 monomer units. This diagram contains regions of stability of coil, two spherical globules (liquid and solid), and rod-like globule conformations. Transitions between the globular states are rounded first-order ones since the states of liquid, solid, and cylinder-like globules do have different internal symmetry.

Serial replica exchange

Parallel tempering (or the replica exchange method (REM)) is a powerful method for speeding up the sampling of conformational states of systems with rough energy landscapes, like proteins, where stable conformational states can be separated by large energy barriers. The usual implementation of the REM is performed on local computer clusters (or parallel processors) where the different replicas must be run synchronously. Here, we present serial replica exchange (SREM), a method that is equivalent to the standard REM in terms of efficiency yet runs asynchronously on a distributed network of computers. A second advantage is the method's greatly enhanced fault tolerance, which enables the study of biological systems on worldwide distributed computing environments, such as Folding@Home. For proof of concept, we apply the SREM to a single alanine dipeptide molecule in explicit water. We show that the SREM reproduces the thermodynamic and structural properties determined by the REM.

Entropic sampling of flexible polyelectrolytes within the Wang-Landau algorithm

We extend the Monte Carlo methods developed in our previous papers [J. Phys. A 37, 1573 (2004); Macromol. Theory-Simul. 14, 491 (2005)] and based on entropic sampling within the Wang-Landau algorithm to simulation of a lattice model of flexible polyelectrolytes. We consider a strongly charged polyelectrolyte chain accompanied by neutralizing counterions on a simple cubic lattice with periodic boundary conditions. The Coulomb potential and the excluded volume condition between different ions or beads are taken into account. The obtained energy distributions make possible the calculation of canonical properties such as conformational energy, heat capacity, entropy, free energy, and mean-square end-to-end distance over a wide temperature range in a single simulation.

Extended multicanonical method combined with thermodynamically optimized potential: Application to the liquid-crystal transition of silicon

A novel method is proposed to study first-order phase transition in real materials. It is applied to the liquid-crystal transition of silicon successfully. It consists of two parts: a direct simulation of the transition by an extended multicanonical ensemble with an order parameter defined with structure factors that characterize the transition, and optimization of a model interatomic potential in terms of the ensemble from an accurate one. These provide a principle to project a first-principles approach on a model-based approach conserving thermodynamic properties of multiple phases.

Accurate estimation of the density of states from Monte Carlo transition probability data

This study develops an efficient approach for calculating the density of states from energy transition probability matrices generated from extended sampling Monte Carlo simulations. Direct and iterative variants of the method are shown to achieve high accuracy when applied to the two-dimensional Ising model for which the density of states function can be determined exactly. They are also used to calculate the density of states of lattice protein and Lennard-Jones models which generate more complex nonzero matrix structures. Whereas the protein simulations test the method on a system exhibiting a rugged free energy landscape, the Lennard-Jones calculations highlight implementation details that arise in applications to continuous energy systems. Density of states results for these two systems agree with estimates from multiple histogram reweighting, demonstrating that the new method provides an alternative approach for computing the thermodynamic properties of complex systems.

Conformational study of Met-enkephalin based on the ECEPP force fields

We report a computational study of the small peptide Met-enkephalin based on the ECEPP/2 and ECEPP/3 force fields using the basin paving method. We have located a new global minimum when using the ECEPP/3 force field with peptide angles omega fixed at 180 degrees. With this new result, we can conclude that the lowest energy configurations of Met-enkephalin predicted based on all four versions of ECEPP have a classic gamma-turn centered at residue Gly3 and a beta-turn at residues Gly3-Phe4. However, minor differences between the structures also exist.

Free-energy calculations with multiple Gaussian modified ensembles

We present a Monte Carlo algorithm, which samples free energies of complex systems. Less probable configurations are populated with the help of a multitude of additional Gaussian weights and parallel tempering is used for efficient Monte Carlo moves within phase space. The algorithm is easily parallelized and can be applied to a wide class of problems. We discuss algorithmic performance for the case of low-temperature phase separation in two-dimensional and three-dimensional Ising models, where we determine the magnetic interface tension. Multiple Gaussian modified ensemble simulations, unlike multicanonical ensemble simulations do not require a priori knowledge of the free energy and are of similar efficiency as multicanonical ensemble and Wang-Landau simulations.

Efficient and precise solvation free energies via alchemical adiabatic molecular dynamics

A new molecular dynamics method for calculating free energies associated with transformations of the thermodynamic state or chemical composition of a system (also known as alchemical transformations) is presented. The new method extends the adiabatic dynamics approach recently introduced by Rosso et al. [J. Chem. Phys. 116, 4389 (2002)] and is based on the use of an additional degree of freedom, lambda, that is used as a switching parameter between the potential energy functions that characterize the two states. In the new method, the coupling parameter lambda is introduced as a fictitious dynamical variable in the Hamiltonian, and a system of switching functions is employed that leads to a barrier in the lambda free energy profile between the relevant thermodynamic end points. The presence of such a barrier, therefore, enhances sampling in the end point (lambda = 0 and lambda = 1) regions which are most important for computing relevant free energy differences. In order to ensure efficient barrier crossing, a high temperature T(lambda) is assigned to lambda and a fictitious mass m(lambda) is introduced as a means of creating an adiabatic separation between lambda and the rest of the system. Under these conditions, it is shown that the lambda free energy profile can be directly computed from the adiabatic probability distribution function of lambda without any postprocessing or unbiasing of the output data. The new method is illustrated on two model problems and in the calculation of the solvation free energy of amino acid side-chain analogs in TIP3P water. Comparisons to previous work using thermodynamic integration and free energy perturbation show that the new lambda adiabatic free energy dynamics method results in very precise free energy calculations using significantly shorter trajectories.

Statistical-temperature Monte Carlo and molecular dynamics algorithms

A simulation method is presented that achieves a flat energy distribution by updating the statistical temperature instead of the density of states in Wang-Landau sampling. A novel molecular dynamics algorithm (STMD) applicable to complex systems and a Monte Carlo algorithm are developed from this point of view. Accelerated convergence for large energy bins, essential for large systems, is demonstrated in tests on the Ising model, the Lennard-Jones fluid, and bead models of proteins. STMD shows a superior ability to find local minima in proteins and new global minima are found for the 55 bead AB model in two and three dimensions. Calculations of the occupation probabilities of individual protein inherent structures provide new insights into folding and misfolding.

Annealing contour Monte Carlo algorithm for structure optimization in an off-lattice protein model

We present a space annealing version for a contour Monte Carlo algorithm and show that it can be applied successfully to finding the ground states for an off-lattice protein model. The comparison shows that the algorithm has made a significant improvement over the pruned-enriched-Rosenbluth method and the Metropolis Monte Carlo method in finding the ground states for AB models. For all sequences, the algorithm has renewed the putative ground energy values in the two-dimensional AB model and set the putative ground energy values in the three-dimensional AB model.

On the use of Bennett's acceptance ratio method in multi-canonical-type simulations

A common strategy for mapping coexistence curves is to employ multi-canonical (MUCA) sampling to simulate along a macrostate path connecting two phases. Central to this approach is the task of accurately calculating the importance weights used in the MUCA procedure, which are needed for both effective sampling and accurate determination of phase boundaries. The purpose of this study is to develop a strategy for determining the importance weights that is built upon Bennett's optimized acceptance ratio method. This approach is shown to be closely related to transition matrix schemes, and is used to compute the vapor-liquid equilibrium of a Lennard-Jones fluid and the liquid-liquid equilibrium of a n-hexane/n-perfluorohexane mixture. For the Lennard-Jones system, the importance weights as a function of the number of particles "N" (at fixed temperature and volume) are obtained by using Bennett's method to estimate free energy differences between N and N+1 particle systems over the desired range of N values. In this application, the method is found to perform slightly better than a related transition matrix scheme. For the n-hexane/n-perfluorohexane liquid mixture, the method is designed to obtain weights as a function of composition (for fixed temperature, pressure, and total number of particles); in this case, the method is found to outperform the Gibbs ensemble approach.

Molecular simulation of the reversible mechanical unfolding of proteins

In this work we have combined a Wang-Landau sampling scheme [F. Wang and D. Landau, Phys. Rev. Lett. 86, 2050 (2001)] with an expanded ensemble formalism to yield a simple and powerful method for computing potentials of mean force. The new method is implemented to investigate the mechanical deformation of proteins. Comparisons are made with analytical results for simple model systems such as harmonic springs and Rouse chains. The method is then illustrated on a model 15-residue alanine molecule in an implicit solvent. Results for mechanical unfolding of this oligopeptide are compared to those of steered molecular dynamics calculations.

Structure of molecular liquids: cavity and bridge functions of the hard spheroid fluid

We present methodologies for calculating the direct correlation function c(1,2), the cavity function y(1,2), and the bridge function b(1,2), for molecular liquids, from Monte Carlo simulations. As an example we present results for the isotropic hard spheroid fluid with elongation e = 3. The simulation data are compared with the results from integral equation theory. In particular, we solve the Percus-Yevick and hypernetted chain equations. In addition, we calculate the first two terms in the virial expansion of the bridge function and incorporate this into the closure. At low densities, the bridge functions calculated by theory and from simulation are in good agreement, lending support to the correctness of our numerical procedures. At higher densities, the hypernetted chain results are brought into closer agreement with simulation by incorporating the approximate bridge function, but significant discrepancies remain.

Universal features and tail analysis of the order-parameter distribution of the two-dimensional Ising model: an entropic sampling Monte Carlo study

We present a numerical study of the order-parameter probability density function (PDF) of the square Ising model for lattices with linear sizes L = 80 - 140. A recent efficient entropic sampling scheme, combining the Wang-Landau and broad histogram methods and based on the high levels of the Wang-Landau process in dominant energy subspaces is employed. We find that for large lattices there exists a stable window of the scaled order-parameter in which the full ansatz including the pre-exponential factor for the tail regime of the universal PDF is well obeyed. This window is used to estimate the equation of state exponent and to observe the behavior of the universal constants implicit in the functional form of the universal PDF. The probability densities are used to estimate the universal Privman-Fisher coefficient and to investigate whether one could obtain reliable estimates of the universal constants controlling the asymptotic behavior of the tail regime.

Quantum Monte Carlo simulation in the canonical ensemble at finite temperature

A quantum Monte Carlo method with a nonlocal update scheme is presented. The method is based on a path-integral decomposition and a worm operator which is local in imaginary time. It generates states with a fixed number of particles and respects other exact symmetries. Observables like the equal-time Green's function can be evaluated in an efficient way. To demonstrate the versatility of the method, results for the one-dimensional Bose-Hubbard model and a nuclear pairing model are presented. Within the context of the Bose-Hubbard model the efficiency of the algorithm is discussed.

Simulation of the density of states in isothermal and adiabatic ensembles

This paper provides a unified treatment of the fundamental methods used to obtain the density of states via molecular simulations with isothermal ensembles (IEs) and adiabatic ensembles (AEs). Our analysis and results show that provides a natural bridge to go back and forth between IE and AE simulation data. They also underline the difference between the density of states of potential energy macrostates and that of total energy macrostates Omega, even though both provide access to the thermodynamic properties of the system. Visited-states approaches and transition matrix methods are described and applied to the Lennard-Jones fluid to target omega and Omega as functions of energy and volume macrostates. It is shown that one can obtain omega via a generalized acceptance-ratio formula that is applicable regardless of the conditions at which the ensemble is simulated. In this way, one can obtain while performing conventional IE or AE simulations, and do it at no extra cost and with a higher accuracy than is achievable with histogram methods.

Performances of Wang-Landau algorithms for continuous systems

The relative performances of different implementations of the Wang-Landau method are assessed on two classes of systems with continuous degrees of freedom, namely, two polypeptides and two atomic Lennard-Jones clusters. Parallel tempering Monte Carlo simulations serve as a reference, and we pay particular attention to the variations of the multiplicative factor f during the course of the simulation. For the systems studied, the Wang-Landau method is found to be of comparable accuracy as parallel tempering, but has significant difficulties in reproducing low-temperature transitions exhibited by the Lennard-Jones clusters at low temperature. Using a complementary order parameter and calculating a two-dimensional joint density of states significantly improves the situation, especially for the notoriously difficult LJ(38) system. However, while parallel tempering easily converges for LJ(31), we have not been able to get data of comparable accuracy with Wang-Landau multicanonical sampling.

Improved Wang-Landau sampling through the use of smoothed potential-energy surfaces

A method is presented to improve the speed of convergence of Wang-Landau simulations as used to calculate the density of states of continuous systems. The density of states is first crudely estimated with calculations employing a smoothed potential-energy surface. This estimate is then used as a seed for subsequent Wang-Landau simulations using the original potential. The performance of the method is demonstrated by employing several simple models, including an analytically solvable harmonic system as well as a Go model of a protein. For all systems considered, the seeded simulations were found to converge significantly faster and with higher accuracy than the standard Wang-Landau simulations.

Wang-Landau algorithm for continuous models and joint density of states

We present a modified Wang-Landau algorithm for models with continuous degrees of freedom. We demonstrate this algorithm with the calculation of the joint density of states of ferromagnet Heisenberg models and a model polymer chain. The joint density of states contains more information than the density of states of a single variable-energy, but is also much more time consuming to calculate. We present strategies to significantly speed up this calculation for large systems over a large range of energy and order parameter.

Free energy landscape analysis of two-dimensional dipolar solvent model at temperatures below and above the rotational freezing point

Ionic solvation in a polar solvent is modeled by a central charge surrounded by dipolar molecules posted on two-dimensional distorted lattice sites with simple rotational dynamics. Density of states is calculated by applying the Wang-Landau algorithm to both the energy and polarization states. The free energy landscapes of solvent molecules as a function of polarization are depicted to explore the competition between the thermal fluctuation and solvation energy. Without a central charge, for temperatures higher than the energy scale of the dipole-dipole interactions, the energy landscape for the small polarization region exhibits a parabolic shape as predicted by Marcus [Rev. Mod. Phys. 65, 599 (1993)] for electron transfer reaction, while there is an additional quartic contribution to the landscape for the large polarization region. When the temperature drops, the simulated free energy landscapes are no longer smooth due to the presence of multiple local minima arising from the frustrated interaction among the dipoles. The parabolic contribution becomes negligible and the energy landscape becomes quartic in shape. For a strong central charge, the energy landscape exhibits an asymmetric profile due to the contributions of linear and cubic terms that arise from the charge-dipole interactions.

Wang-Landau Monte Carlo simulation of the Blume-Capel model

We carry out a study of the two-dimensional Blume-Capel model using the Wang-Landau Monte Carlo method which estimates the density of states g(E) directly. This work validates the applicability of this method to multiparametric systems, since only one computer run is needed for all range of macroscopic parameters (temperature, anisotropy, etc.). The location of the tricritical point is determined as kBTt/J=0.609(3), Dt/J=1.966(2), and is in excellent agreement with previous estimates. The free energy and the entropy, which are not directly accessible by conventional Monte Carlo simulations, are obtained simply using g(E).