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

Comparing parallel- and simulated-tempering-enhanced sampling algorithms at phase-transition regimes

Two important enhanced sampling algorithms, simulated (ST) and parallel (PT) tempering, are commonly used when ergodic simulations may be hard to achieve, e.g., due to a phase space separated by large free-energy barriers. This is so for systems around first-order phase transitions, a case still not fully explored with such approaches in the literature. In this contribution we make a comparative study between the PT and ST for the Ising (a lattice gas in the fluid language) and the Blume-Emery-Griffiths (a lattice gas with vacancies) models at phase-transition regimes. We show that although the two methods are equivalent in the limit of sufficiently long simulations, the PT is more advantageous than the ST with respect to all the analysis performed: convergence toward the stationarity; frequency of tunneling between phases at the coexistence; and decay of time-displaced correlation functions of thermodynamic quantities. Qualitative arguments for why one may expect better results from the PT than the ST near phase-transitions conditions are also presented.

Generalized replica exchange method

We present a powerful replica exchange method, particularly suited to first-order phase transitions associated with the backbending in the statistical temperature, by merging an optimally designed generalized ensemble sampling with replica exchanges. The key ingredients of our method are parametrized effective sampling weights, smoothly joining ordered and disordered phases with a succession of unimodal energy distributions by transforming unstable or metastable energy states of canonical ensembles into stable ones. The inverse mapping between the sampling weight and the effective temperature provides a systematic way to design the effective sampling weights and determine a dynamic range of relevant parameters. Illustrative simulations on Potts spins with varying system size and simulation conditions demonstrate a comprehensive sampling for phase-coexistent states with a dramatic acceleration of tunneling transitions. A significant improvement over the power-law slowing down of mean tunneling times with increasing system size is obtained, and the underlying mechanism for accelerated tunneling is discussed.

Metadynamics convergence law in a multidimensional system

Metadynamics is a powerful sampling technique that uses a nonequilibrium history-dependent process to reconstruct the free-energy surface as a function of the relevant collective variables s . In Bussi [Phys. Rev. Lett. 96, 090601 (2006)] it is proved that, in a Langevin process, metadynamics provides an unbiased estimate of the free energy F(s) . We here study the convergence properties of this approach in a multidimensional system, with a Hamiltonian depending on several variables. Specifically, we show that in a Monte Carlo metadynamics simulation of an Ising model the time average of the history-dependent potential converge to F(s) with the same law of an umbrella sampling performed in optimal conditions (i.e., with a bias exactly equal to the negative of the free energy). Remarkably, after a short transient, the error becomes approximately independent on the filling speed, showing that even in out-of-equilibrium conditions metadynamics allows recovering an accurate estimate of F(s) . These results have been obtained introducing a functional form of the history-dependent potential that avoids the onset of systematic errors near the boundaries of the free-energy landscape.

Analytic continuation of quantum Monte Carlo data by stochastic analytical inference

We present an algorithm for the analytic continuation of imaginary-time quantum Monte Carlo data which is strictly based on principles of Bayesian statistical inference. Within this framework we are able to obtain an explicit expression for the calculation of a weighted average over possible energy spectra, which can be evaluated by standard Monte Carlo simulations, yielding as by-product also the distribution function as function of the regularization parameter. Our algorithm thus avoids the usual ad hoc assumptions introduced in similar algorithms to fix the regularization parameter. We apply the algorithm to imaginary-time quantum Monte Carlo data and compare the resulting energy spectra with those from a standard maximum-entropy calculation.

Critical behavior of the pure and random-bond two-dimensional triangular Ising ferromagnet

We investigate the effects of quenched bond randomness on the critical properties of the two-dimensional ferromagnetic Ising model embedded in a triangular lattice. The system is studied in both the pure and disordered versions by the same efficient two-stage Wang-Landau method. In the first part of our study, we present the finite-size scaling behavior of the pure model, for which we calculate the critical amplitude of the specific heat's logarithmic expansion. For the disordered system, the numerical data and the relevant detailed finite-size scaling analysis along the lines of the two well-known scenarios-logarithmic corrections versus weak universality--strongly support the field-theoretically predicted scenario of logarithmic corrections. A particular interest is paid to the sample-to-sample fluctuations of the random model and their scaling behavior that are used as a successful alternative approach to criticality.

Multicritical points and crossover mediating the strong violation of universality: Wang-Landau determinations in the random-bond d=2 Blume-Capel model

The effects of bond randomness on the phase diagram and critical behavior of the square lattice ferromagnetic Blume-Capel model are discussed. The system is studied in both the pure and disordered versions by the same efficient two-stage Wang-Landau method for many values of the crystal field, restricted here in the second-order phase-transition regime of the pure model. For the random-bond version several disorder strengths are considered. We present phase diagram points of both pure and random versions and for a particular disorder strength we locate the emergence of the enhancement of ferromagnetic order observed in an earlier study in the ex-first-order regime. The critical properties of the pure model are contrasted and compared to those of the random model. Accepting, for the weak random version, the assumption of the double-logarithmic scenario for the specific heat we attempt to estimate the range of universality between the pure and random-bond models. The behavior of the strong disorder regime is also discussed and a rather complex and yet not fully understood behavior is observed. It is pointed out that this complexity is related to the ground-state structure of the random-bond version.

The Widom-Rowlinson mixture on a sphere: elimination of exponential slowing down at first-order phase transitions

Computer simulations of first-order phase transitions using 'standard' toroidal boundary conditions are generally hampered by exponential slowing down. This is partly due to interface formation, and partly due to shape transitions. The latter occur when droplets become large such that they self-interact through the periodic boundaries. On a spherical simulation topology, however, shape transitions are absent. We expect that by using an appropriate bias function, exponential slowing down can be largely eliminated. In this work, these ideas are applied to the two-dimensional Widom-Rowlinson mixture confined to the surface of a sphere. Indeed, on the sphere, we find that the number of Monte Carlo steps needed to sample a first-order phase transition does not increase exponentially with system size, but rather as a power law τ α V(α), with α≈2.5, and V the system area. This is remarkably close to a random walk for which α(RW) = 2. The benefit of this improved scaling behavior for biased sampling methods, such as the Wang-Landau algorithm, is investigated in detail.

Finite size scaling and first-order phase transition in a modified XY model

Monte Carlo simulation has been performed in a two-dimensional modified XY -model first proposed by Domany [Phys. Rev. Lett. 52, 1535 (1984)] The cluster algorithm of Wolff has been used and multiple histogram reweighting is performed. The first-order scaling behavior of the quantities such as specific heat and free-energy barrier are found to be obeyed accurately. While the lowest-order correlation function was found to decay to zero at long distance just above the transition, the next-higher-order correlation function shows a nonzero plateau.

Isotropic-to-nematic transition in confined liquid crystals: an essentially nonuniversal phenomenon

Computer simulations are presented of the isotropic-to-nematic transition in a liquid crystal confined between two parallel plates a distance H apart. The plates are neutral and do not impose any anchoring on the particles. Depending on the shape of the pair potential acting between the particles, we find that the transition either changes from first order to continuous at a critical film thickness H=H(x) , or that the transition remains first order irrespective of H . This demonstrates that the isotropic-to-nematic transition in confined geometry is not characterized by any universality class, but rather that its fate is determined by microscopic details. The resulting capillary phase diagrams can thus assume two topologies: one where the isotropic and nematic branches of the binodal meet at H=H(x), and one where they remain separated. For values of H where the transition is strongly first order the shift Deltaepsilon of the transition temperature is in excellent agreement with the Kelvin equation. Not only is the relation Deltaepsilon proportional, variant 1/H recovered but also the prefactor of the shift is in quantitative agreement with the independently measured bulk latent heat and interfacial tension.

Phase transitions of a single polymer chain: A Wang-Landau simulation study

A single flexible homopolymer chain can assume a variety of conformations which can be broadly classified as expanded coil, collapsed globule, and compact crystallite. Here we study transitions between these conformational states for an interaction-site polymer chain comprised of N=128 square-well-sphere monomers with hard-sphere diameter sigma and square-well diameter lambdasigma. Wang-Landau sampling with bond-rebridging Monte Carlo moves is used to compute the density of states for this chain and both canonical and microcanonical analyses are used to identify and characterize phase transitions in this finite size system. The temperature-interaction range (i.e., T-lambda) phase diagram is constructed for lambda<or=1.30. Chains assume an expanded coil conformation at high temperatures and a crystallite structure at low temperatures. For lambda>1.06 these two states are separated by an intervening collapsed globule phase and thus, with decreasing temperature a chain undergoes a continuous coil-globule (collapse) transition followed by a discontinuous globule-crystal (freezing) transition. For well diameters lambda<1.06 the collapse transition is pre-empted by the freezing transition and thus there is a direct first-order coil-crystal phase transition. These results confirm the recent prediction, based on a lattice polymer model, that a collapsed globule state is unstable with respect to a solid phase for flexible polymers with sufficiently short-range monomer-monomer interactions.

Introducing sampling entropy in repository based adaptive umbrella sampling

Determining free energy surfaces along chosen reaction coordinates is a common and important task in simulating complex systems. Due to the complexity of energy landscapes and the existence of high barriers, one widely pursued objective to develop efficient simulation methods is to achieve uniform sampling among thermodynamic states of interest. In this work, we have demonstrated sampling entropy (SE) as an excellent indicator for uniform sampling as well as for the convergence of free energy simulations. By introducing SE and the concentration theorem into the biasing-potential-updating scheme, we have further improved the adaptivity, robustness, and applicability of our recently developed repository based adaptive umbrella sampling (RBAUS) approach [H. Zheng and Y. Zhang, J. Chem. Phys. 128, 204106 (2008)]. Besides simulations of one dimensional free energy profiles for various systems, the generality and efficiency of this new RBAUS-SE approach have been further demonstrated by determining two dimensional free energy surfaces for the alanine dipeptide in gas phase as well as in water.

Phase diagram and critical behavior of the square-lattice Ising model with competing nearest-neighbor and next-nearest-neighbor interactions

Using the parallel tempering algorithm and graphics processing unit accelerated techniques, we have performed large-scale Monte Carlo simulations of the Ising model on a square lattice with antiferromagnetic (repulsive) nearest-neighbor and next-nearest-neighbor interactions of the same strength and subject to a uniform magnetic field. Both transitions from the (2x1) and row-shifted (2x2) ordered phases to the paramagnetic phase are continuous. From our data analysis, re-entrance behavior of the (2x1) critical line and a bicritical point which separates the two ordered phases at T=0 are confirmed. Based on the critical exponents we obtained along the phase boundary, Suzuki's weak universality seems to hold.

Six-state clock model on the square lattice: fisher zero approach with Wang-Landau sampling

We investigate the six-state clock model with nearest-neighbor interactions on the square lattice. We obtain the density of states of the finite systems up to L=28 using the Wang-Landau sampling. With the density of states and the Fisher zero approach, we successfully find two different critical temperatures 0.632(2) and 0.997(2) for the clock model. It seems that this study supports the recent results by [Lapilli Phys. Rev. Lett. 96, 140603 (2006)] that the transitions are not of Kosterlitz and Thouless type.

Some recent techniques for free energy calculations

A few recent techniques to calculate free energies in the context of molecular dynamics simulations are discussed: temperature-accelerated molecular dynamics, which is a method to explore fast the important regions in the free energy landscape associated with a set of continuous collective variables without having to know where these regions are beforehand; the single sweep method, which is a variational method to interpolate the free energy globally given a set of mean forces (i.e., a set of gradients of the free energy) calculated at specific points, or centers, on the free energy landscape; and a Voronoi-based free energy method for the calculation of the free energy of the Voronoi tessellation associated with a set of centers. We also discuss how this last technique can be used in conjunction with the string method, and how kinetic information such as reaction rates can be calculated by milestoning using the edges of a Voronoi tessellation as milestones.

Temperature and anharmonic effects on the infrared absorption spectrum from a quantum statistical approach: application to naphthalene

A method is developed to calculate the finite-temperature infrared absorption spectrum of polyatomic molecules with energy levels described by a second-order Dunham expansion. The anharmonic couplings are fully incorporated in the calculation of the quantum density of states, achieved using a Wang-Landau Monte Carlo procedure, as well as in the determination of transition energies. Additional multicanonical simulations provide the microcanonical absorption intensity as a function of both the absorption wavelength and the internal energy of the molecule. The finite-temperature spectrum is finally obtained by Laplace transformation of this microcanonical histogram. The present scheme is applied to the infrared spectrum of naphthalene, for which we quantify the shifting, broadening, and third-order effects as a continuous function of temperature. The influence of anharmonicity and couplings is manifested on the nontrivial variations of these features with increasing temperature.

Replica exchange statistical temperature Monte Carlo

The replica exchange statistical temperature Monte Carlo algorithm (RESTMC) is presented, extending the single-replica STMC algorithm [J. Kim, J. E. Straub, and T. Keyes, Phys. Rev. Lett. 97, 050601 (2006)] to alleviate the slow convergence of the conventional temperature replica exchange method (t-REM) with increasing system size. In contrast to the Gibbs-Boltzmann sampling at a specific temperature characteristic of the standard t-REM, RESTMC samples a range of temperatures in each replica and achieves a flat energy sampling employing the generalized sampling weight, which is automatically determined via the dynamic modification of the replica-dependent statistical temperature. Faster weight determination, through the dynamic update of the statistical temperature, and the flat energy sampling, maximizing energy overlaps between neighboring replicas, lead to a considerable acceleration in the convergence of simulations even while employing significantly fewer replicas. The performance of RESTMC is demonstrated and quantitatively compared with that of the conventional t-REM under varying simulation conditions for Lennard-Jones 19, 31, and 55 atomic clusters, exhibiting single- and double-funneled energy landscapes.

A one-dimensional dipole lattice model for water in narrow nanopores

We present a recently developed one-dimensional dipole lattice model that accurately captures the key properties of water in narrow nanopores. For this model, we derive three equivalent representations of the Hamiltonian that together yield a transparent physical picture of the energetics of the water chain and permit efficient computer simulations. In the charge representation, the Hamiltonian consists of nearest-neighbor interactions and Coulomb-like interactions of effective charges at the ends of dipole ordered segments. Approximations based on the charge picture shed light on the influence of the Coulomb-like interactions on the structure of nanopore water. We use Monte Carlo simulations to study the system behavior of the full Hamiltonian and its approximations as a function of chemical potential and system size and investigate the bimodal character of the density distribution occurring at small system sizes.

Calculating thermodynamics properties of quantum systems by a non-Markovian Monte Carlo procedure

We present a history-dependent Monte Carlo scheme for the efficient calculation of the free energy of quantum systems inspired by Wang-Landau and metadynamics. In the two-dimensional quantum Ising model, chosen here for illustration, the accuracy of free energy, critical temperature, and specific heat is demonstrated as a function of simulation time and successfully compared with the best available approaches. The approach is based on a path integral formulation of the quantum problem and can be applied without modifications to quantum Hamiltonians of any level of complexity. The combination of high accuracy and performance with a much broader applicability is a major advance with respect to other available methods.

Effect of platykurtic and leptokurtic distributions in the random-field Ising model: mean-field approach

The influence of the tail features of the local magnetic field probability density function (PDF) on the ferromagnetic Ising model is studied in the limit of infinite range interactions. Specifically, we assign a quenched random field whose value is in accordance with a generic distribution that bears platykurtic and leptokurtic distributions depending on a single parameter tau<3 to each site. For tau<5/3, such distributions, which are basically Student-t and r distribution extended for all plausible real degrees of freedom, present a finite standard deviation, if not the distribution has got the same asymptotic power-law behavior as a alpha-stable Lévy distribution with alpha=(3-tau)/(tau-1). For every value of tau, at specific temperature and width of the distribution, the system undergoes a continuous phase transition. Strikingly, we impart the emergence of an inflexion point in the temperature-PDF width phase diagrams for distributions broader than the Cauchy-Lorentz (tau=2) which is accompanied with a divergent free energy per spin (at zero temperature).

Application of the Wang-Landau algorithm to the dimerization of glycophorin A

A two-step Monte Carlo procedure is developed to investigate the dimerization process of the homodimer glycophorin A. In the first step, the energy density of states of the system is estimated by the Wang-Landau algorithm. In the second step, a production run is performed during which various energetical and structural observables are sampled to provide insight into the thermodynamics of the system. All seven residues LIxxGVxxGVxxT constituting the contact interface play a dominating role in the dimerization, however at different stages of the process. The leucine motif and to some extent the GxxxG motif are involved at the very beginning of the dimerization when the two helices come into contact, ensuring an interface already similar to the native one. At a lower temperature, the threonine motif stabilizes by hydrogen bonding the dimer, which finally converges toward its native state at around 300 K. The power and flexibility of the procedure employed here makes it an interesting alternative to other Monte Carlo methods for the study of similar protein systems.

Crossover from first- to second-order transition in frustrated Ising antiferromagnetic films

In the bulk state, the Ising face-centered-cubic (fcc) antiferromagnet is fully frustrated and is known to have a very strong first-order transition. In this paper, we study the nature of this phase transition in the case of a thin film as a function of the film thickness. Using Monte Carlo simulations, we show that the transition remains first order down to a thickness of four fcc cells (eight atomic layers). It becomes clearly second order at a thickness of two fcc cells, i.e., four atomic layers. It is also interesting to note that the presence of the surface reduces the ground-state degeneracy found in the bulk. For the two-cell thickness, the surface magnetization is larger than the interior one. It undergoes a second-order phase transition at a temperature TC while interior spins become disordered at a lower temperature TD. This loss of order is characterized by a peak of the interior spins susceptibility and a peak of the specific heat which do not depend on the lattice size suggesting that either it is not a real transition or it is a Kosterlitz-Thouless nature. The surface transition, on the other hand, is shown to be of second order with critical exponents deviated from those of pure two-dimensional Ising universality class. We also show results obtained from the Green's function method. A discussion is given.

Enhanced sampling in generalized ensemble with large gap of sampling parameter: case study in temperature space random walk

We present an efficient sampling method for computing a partition function and accelerating configuration sampling. The method performs a random walk in the lambda space, with lambda being any thermodynamic variable that characterizes a canonical ensemble such as the reciprocal temperature beta or any variable that the Hamiltonian depends on. The partition function is determined by minimizing the difference of the thermal conjugates of lambda (the energy in the case of lambda = beta), defined as the difference between the value from the dynamically updated derivatives of the partition function and the value directly measured from simulation. Higher-order derivatives of the partition function are included to enhance the Brownian motion in the lambda space. The method is much less sensitive to the system size, and to the size of lambda window than other methods. On the two dimensional Ising model, it is shown that the method asymptotically converges the partition function, and the error of the logarithm of the partition function is much smaller than the algorithm using the Wang-Landau recursive scheme. The method is also applied to off-lattice model proteins, the AB models, in which cases many low energy states are found in different models.

Versatile approach to access the low temperature thermodynamics of lattice polymers and proteins

We show that Wang-Landau sampling, combined with suitable Monte Carlo trial moves, provides a powerful method for both the ground state search and the determination of the density of states for the hydrophobic-polar (HP) protein model and the interacting self-avoiding walk (ISAW) model for homopolymers. We obtain accurate estimates of thermodynamic quantities for HP sequences with >100 monomers and for ISAWs up to >500 monomers. Our procedure possesses an intrinsic simplicity and overcomes the limitations inherent in more tailored approaches making it interesting for a broad range of protein and polymer models.

All-or-none proteinlike folding transition of a flexible homopolymer chain

Here we report a first-order all-or-none transition from an expanded coil to a compact crystallite for a flexible polymer chain. Wang-Landau sampling is used to construct the complete density of states for square-well chains up to length 256. Analysis within both the microcanonical and canonical ensembles shows a direct freezing transition for finite length chains with sufficiently short-range interactions. This type of transition is a distinctive feature of "one-step" protein folding and our findings demonstrate that a simple homopolymer model can exhibit protein-folding thermodynamics.

Energy landscape of a spin-glass model: exploration and characterization

The disconnectivity graph (DG) is widely used to represent energy landscapes. Although powerful numerical methods have been developed to construct DGs for continuous potential-energy surfaces, they have difficulties in applications to discrete Hamiltonians as the case of spin-glass models. When the configuration space is large, brute force enumeration of all configurations to build a DG is not practical. We propose an alternative approach to construct DGs based on recursive partition of Monte Carlo samples from microcanonical ensembles. To characterize energy landscapes, we define the local density of states (LDOS) on a DG, with which one can compute many thermodynamic properties over local energy basins for any temperature. Estimation of LDOS is developed with DG construction. We further propose the concepts of tree entropy and local escape probability, both of which are functions of local density of states, to capture the symmetry and the roughness of a Boltzmann distribution, respectively. Our approach is applied to a study of the Sherrington-Kirkpatrick spin-glass model with N varying between 20 and 100 spins. We observe that the energy landscape is extremely asymmetric and there exists a sharp increase in local escape probability preceding the transition from spin glass to paramagnetic phase.

An optimized replica exchange molecular dynamics method

We introduce a new way to perform swaps between replicas in replica exchange molecular dynamics simulations. The method is based on a generalized canonical probability distribution function and flattens the potential of mean force along the temperature coordinate, such that a random walk in temperature space is achieved. Application to a Go model of protein A showed that the method is more efficient than conventional replica exchange. The method results in a constant probability distribution of the replicas over the thermostats, yields a minimum round-trip time between extremum temperatures, and leads to faster ergodic convergence.

Thermodynamics of a conformational change using a random walk in energy-reaction coordinate space: Application to methane dimer hydrophobic interactions

A random walk sampling algorithm allows the extraction of the density of states distribution in energy-reaction coordinate space. As a result, the temperature dependences of thermodynamic quantities such as relative energy, entropy, and heat capacity can be calculated using first-principles statistical mechanics. The strategies for optimal convergence of the algorithm and control of its accuracy are proposed. We show that the saturation of the error [Q. Yan and J. J. de Pablo, Phys. Rev. Lett. 90, 035701 (2003); E. Belardinelli and V. D. Pereyra, J. Chem. Phys. 127, 184105 (2007)] is due to the use of histogram flatness as a criterion of convergence. An application of the algorithm to methane dimer hydrophobic interactions is presented. We obtained a quantitatively accurate energy-entropy decomposition of the methane dimer cavity potential. The presented results confirm the previous results, and they provide new information regarding the thermodynamics of hydrophobic interactions. We show that the finite-difference approximation, which is widely used in molecular dynamic simulations for the energy-entropy decomposition of a free energy potential, can lead to a significant error.

Optimal replica exchange method combined with Tsallis weight sampling

A unified framework integrating the generalized ensemble sampling associated with the Tsallis weight [C. Tsallis, J. Stat. Phys. 52, 479 (1988)] and the replica exchange method (REM) has been proposed to accelerate the convergence of the conventional temperature REM (t-REM). Using the effective temperature formulation of the Tsallis weight sampling, it is shown that the average acceptance probability for configurational swaps between neighboring replicas in the combination of Tsallis weight sampling and REM (Tsallis-REM) is directly proportional to an overlap integral of the energy distributions of neighboring replicas as in the t-REM. Based on this observation, we suggest a robust method to select optimal Tsallis parameters in the conventional parametrization scheme and present new parametrization schemes for the Tsallis-REM, which significantly improves the acceptance of configurational swaps by systematically modulating energy overlaps between neighboring replicas. The distinguished feature of our method is that all relevant parameters in the Tsallis-REM are automatically determined from the equilibrium phase simulation using the t-REM. The overall performance of our method is explicitly demonstrated for various simulation conditions for the Lennard-Jones 31 atom clusters, exhibiting a double-funneled energy landscape.

From multidimensional replica-exchange method to multidimensional multicanonical algorithm and simulated tempering

We discuss multidimensional generalizations of multicanonical algorithm, simulated tempering, and replica-exchange method. We generalize the original potential-energy function E0 by adding any physical quantity V of interest as a new energy term with a coupling constant lambda. We then perform a multidimensional multicanonical simulation where a random walk in E0 and V spaces is realized. We can alternately perform a multidimensional simulated-tempering simulation where a random walk in temperature T and parameter lambda is realized. The results of the multidimensional replica-exchange simulations can be used to determine the weight factors for these multidimensional multicanonical and simulated-tempering simulations.

QM/MM methods for biomolecular systems

Combined quantum-mechanics/molecular-mechanics (QM/MM) approaches have become the method of choice for modeling reactions in biomolecular systems. Quantum-mechanical (QM) methods are required for describing chemical reactions and other electronic processes, such as charge transfer or electronic excitation. However, QM methods are restricted to systems of up to a few hundred atoms. However, the size and conformational complexity of biopolymers calls for methods capable of treating up to several 100,000 atoms and allowing for simulations over time scales of tens of nanoseconds. This is achieved by highly efficient, force-field-based molecular mechanics (MM) methods. Thus to model large biomolecules the logical approach is to combine the two techniques and to use a QM method for the chemically active region (e.g., substrates and co-factors in an enzymatic reaction) and an MM treatment for the surroundings (e.g., protein and solvent). The resulting schemes are commonly referred to as combined or hybrid QM/MM methods. They enable the modeling of reactive biomolecular systems at a reasonable computational effort while providing the necessary accuracy.

Momentum-shell renormalization-group flow from simulation

Our recently developed Fourier Monte Carlo algorithm permits a nonperturbative calculation of momentum-shell renormalization-group flows by simulation which despite its apparent simplicity is illustrative both numerically as well as conceptually interesting. We study the example of a varphi(4) model with long-range lattice interaction. For this model we show that the topology of the renormalization flow is globally accessible in a particularly convenient way. The nontrivial fixed point of Wilson-Fisher type is observed its accompanying critical exponents are numerically determined from fitting its surrounding flow pattern to a linearized renormalization-group transformation. The results are compared to those obtained from perturbation theory, -expansion and earlier Monte Carlo simulations. Application of our method is also expected to be rewarding in other models with long-range interactions.

Identification of amyloidogenic peptide sequences using a coarse-grained physicochemical model

Cross-beta amyloid is implicated in over 20 human diseases. Experiments suggest that specific sequence elements within amyloidogenic proteins play a major role in seeding amyloid formation. Identifying these seeding sequences is important for rationalizing the molecular mechanisms of amyloid formation and for elaborating therapeutic strategies that target amyloid. Theoretical techniques play an important role in facilitating the identification and structural characterization of putative seeding sequences; most amyloid species are not amenable to high resolution experimental structure techniques. In this study we have combined a coarse-grained physicochemical protein model with a highly efficient Monte Carlo sampling technique to identify amyloidogenic sequences in four proteins for which respective experimental peptide fragmentation data exist. Peptide sequences were defined as amyloidogenic if the ensemble structure predicted for three interacting peptides described a stable and regular three-stranded beta-sheet. For such peptides, free energies were calculated to provide a measure of amyloid propensity. The overall agreement between the experimental and predicted data is good, and we correctly identify several self-recognition motifs proposed to define the cross-beta amyloid fibril architectures of two of the proteins. Our results compare very favorably with those obtained using atomistic molecular dynamics methods, though our simulations are 30-40 times faster.

WHAT MAKES MOLECULAR DYNAMICS WORK?

The equations of motion for deterministic molecular dynamics (MD) are chaotic, creating problems for their numerical treatment due to the exponential growth of error with time. Indeed, modeling and computational errors overwhelm numerical trajectories in typical simulations. Consequently, accuracy is expected only in a statistical sense, based on random initial conditions. Of great interest then is the relationship between errors in the dynamics and their effects on the accuracy of statistical quantities, specifically, expectations. This article provides a formula for the effect of a perturbation on an ensemble average, which explains the accuracy of such calculations. It also provides a formula for the effect of a perturbation on a time correlation function, which, however, fails to explain accuracy for these calculations. Additionally, this article clarifies the relationships among various dynamical properties of MD and provides an extension to a theory of non-Hamiltonian MD.