Metropolis Monte Carlo Sampling Research Articles

Metropolis Monte Carlo sampling: convergence, localization transition and optimality

Among random sampling methods, Markov chain Monte Carlo (MC) algorithms are foremost. Using a combination of analytical and numerical approaches, we study their convergence properties toward the steady state, within a random walk Metropolis scheme. Analyzing the relaxation properties of some model algorithms sufficiently simple to enable analytic progress, we show that the deviations from the target steady-state distribution can feature a localization transition as a function of the characteristic length of the attempted jumps defining the random walk. While the iteration of the MC algorithm converges to equilibrium for all choices of jump parameters, the localization transition changes drastically the asymptotic shape of the difference between the probability distribution reached after a finite number of steps of the algorithm and the target equilibrium distribution. We argue that the relaxation before and after the localization transition is respectively limited by diffusion and rejection rates.

Iterative cluster expansion approach for predicting the structure evolution of mixed Ruddelsden-Popper oxides La2−xSrxNi1−yFeyO4±δ

Layered perovskites viz. Ruddlesden-Popper (RP) compounds have attracted a significant research interest lately due to their superior electrocatalytic activities, on top of their rich physics. Although the versatility of their compositional space is the key to tune their stability and applicability, its complexity also makes the understanding and design of new functional RP-based materials difficult. Herein by parameterizing a cluster expansion model Hamiltonian with density-functional inputs, we have successfully scanned over the complex configurational space of an RP phase, taking ${\mathrm{La}}_{2\ensuremath{-}x}{\mathrm{Sr}}_{x}{\mathrm{Ni}}_{1\ensuremath{-}y}{\mathrm{Fe}}_{y}{\mathrm{O}}_{4\ifmmode\pm\else\textpm\fi{}\ensuremath{\delta}}$ as an example. This surrogate model is then used to predict the energies of unknown compositions and performing Metropolis Monte Carlo sampling for finding the ground-state structures. By analyzing these stable configurations at different $x, y$'s, and $\ensuremath{\delta}$=0, 0.125, we found that apical vacancies in the rock-salt layer are more favorable due to lesser steric forces but changes in $A$-site oxidation state due to Sr substitution can facilitate equatorial O-vacancy formation even in the perovskite layers. Our results uncover the structural evolution of RP oxides with increasing $A$-site (Sr) and $B$-site (Fe) substitutions. We find that a critical Sr concentration is required to stabilize Fe in these systems. With increasing Sr concentration, the $A$-site charge state changes gradually from ${\mathrm{La}}^{+3}$ to ${\mathrm{Sr}}^{+2}$, which facilitates both O-vacancy formation and stabilization of the $B$-site dopants. The findings further highlight the effects of configurational changes for the same composition on electronic structure. Thus, our work clearly demonstrates the importance of exploring configurational space of these complex oxides when using electronic-structure features as descriptors of their catalytic and other functional properties.

Efficient low temperature Monte Carlo sampling using quantum annealing

Quantum annealing is an efficient technology to determine ground state configurations of discrete binary optimization problems, described through Ising Hamiltonians. Here we show that—at very low computational cost—finite temperature properties can be calculated. The approach is most efficient at low temperatures, where conventional approaches like Metropolis Monte Carlo sampling suffer from high rejection rates and therefore large statistical noise. To demonstrate the general approach, we apply it to spin glasses and Ising chains.

A Metropolis Monte Carlo algorithm for merging single-particle diffraction intensities.

Single-particle imaging with X-ray free-electron lasers depends crucially on algorithms that merge large numbers of weak diffraction patterns despite missing measurements of parameters such as particle orientations. The expand-maximize-compress (EMC) algorithm is highly effective at merging single-particle diffraction patterns with missing orientation values, but most implementations exhaustively sample the space of missing parameters and may become computationally prohibitive as the number of degrees of freedom extends beyond orientation angles. This paper describes how the EMC algorithm can be modified to employ Metropolis Monte Carlo sampling rather than grid sampling, which may be favorable for reconstruction problems with more than three missing parameters. Using simulated data, this variant is compared with the standard EMC algorithm.

Q-Canonical Monte Carlo Sampling for Modeling the Linkage between Charge Regulation and Conformational Equilibria of Peptides.

The overall charge content and the patterning of charged residues have a profound impact on the conformational ensembles adopted by intrinsically disordered proteins. These parameters can be altered by charge regulation, which refers to the effects of post-translational modifications, pH-dependent changes to charge, and conformational fluctuations that modify the pKa values of ionizable residues. Although atomistic simulations have played a prominent role in uncovering the major sequence-ensemble relationships of IDPs, most simulations assume fixed charge states for ionizable residues. This may lead to erroneous estimates for conformational equilibria if they are linked to charge regulation. Here, we report the development of a new method we term q-canonical Monte Carlo sampling for modeling the linkage between charge regulation and conformational equilibria. The method, which is designed to be interoperable with the ABSINTH implicit solvation model, operates as follows: For a protein sequence with n ionizable residues, we start with all 2n charge microstates and use a criterion based on model compound pKa values to prune down to a subset of thermodynamically relevant charge microstates. This subset is then grouped into mesostates, where all microstates that belong to a mesostate have the same net charge. Conformational distributions, drawn from a canonical ensemble, are generated for each of the charge microstates that make up a mesostate using a method we designate as proton walk sampling. This method combines Metropolis Monte Carlo sampling in conformational space with an auxiliary Markov process that enables interconversions between charge microstates along a mesostate. Proton walk sampling helps identify the most likely charge microstate per mesostate. We then use thermodynamic integration aided by the multistate Bennett acceptance ratio method to estimate the free energies for converting between mesostates. These free energies are then combined with the per-microstate weights along each mesostate to estimate standard state free energies and pH-dependent free energies for all thermodynamically relevant charge microstates. The results provide quantitative estimates of the probabilities and preferred conformations associated with every thermodynamically accessible charge microstate. We showcase the application of q-canonical sampling using two model systems. The results establish the soundness of the method and the importance of charge regulation in systems characterized by conformational heterogeneity.

Hamiltonian Switch Metropolis Monte Carlo Simulations for Improved Conformational Sampling of Intrinsically Disordered Regions Tethered to Ordered Domains of Proteins

There is growing interest in the topic of intrinsically disordered proteins (IDPs). Atomistic Metropolis Monte Carlo (MMC) simulations based on novel implicit solvation models have yielded useful insights regarding sequence-ensemble relationships for IDPs modeled as autonomous units. However, a majority of naturally occurring IDPs are tethered to ordered domains. Tethering introduces additional energy scales and this creates the challenge of broken ergodicity for standard MMC sampling or molecular dynamics that cannot be readily alleviated by using generalized tempering methods. We have designed, deployed, and tested our adaptation of the Nested Markov Chain Monte Carlo sampling algorithm. We refer to our adaptation as Hamiltonian Switch Metropolis Monte Carlo (HS-MMC) sampling. In this method, transitions out of energetic traps are enabled by the introduction of an auxiliary Markov chain that draws conformations for the disordered region from a Boltzmann distribution that is governed by an alternative potential function that only includes short-range steric repulsions and conformational restraints on the ordered domain. We show using multiple, independent runs that the HS-MMC method yields conformational distributions that have similar and reproducible statistical properties, which is in direct contrast to standard MMC for equivalent amounts of sampling. The method is efficient and can be deployed for simulations of a range of biologically relevant disordered regions that are tethered to ordered domains.

MagiC: Software Package for Multiscale Modeling.

We present software package MagiC, which is designed to perform systematic structure-based coarse graining of molecular models. The effective pairwise potentials between coarse-grained sites of low-resolution molecular models are constructed to reproduce structural distribution functions obtained from the modeling of the system in a high resolution (atomistic) description. The software supports coarse-grained tabulated intramolecular bond and angle interactions, as well as tabulated nonbonded interactions between different site types in the coarse-grained system, with the treatment of long-range electrostatic forces by the Ewald summation. Two methods of effective potential refinement are implemented: iterative Boltzmann inversion and inverse Monte Carlo, the latter accounting for cross-correlations between pair interactions. MagiC uses its own Metropolis Monte Carlo sampling engine, allowing parallel simulation of many copies of the system with subsequent averaging of the properties, which provides fast convergence of the method with nearly linear scaling at parallel execution.

Ordering and segregation in isolated Au–Pd icosahedral nanoclusters and nanowires and the consequences of their encapsulation inside carbon nanotubes

Metropolis Monte Carlo sampling in the semi-grand canonical ensemble with empirical potentials is used to predict equilibrium ordered structures and segregation properties of small icosahedral Au–Pd nanoclusters and helical nanowire segments over the whole range of compositions at low temperatures. The cases of free-standing clusters and wires are compared with the same systems encapsulated in carbon nanotubes.A number of chemically ordered structures and segregation states are identified and found to be consistent with the same interplay between size mismatch, mixing enthalpy and surface energies of elemental metals which determines the thermodynamic equilibrium of binary metal alloys. Encapsulation has the effect of modifying the surface energies of nanoclusters and wires, with considerable consequences on their thermodynamic states, although the metal–graphite interaction strength is low as compared with the metal cohesive energy and the carbon–carbon binding energy.

Thermodynamic properties of Au–Pd nanostructured surfaces studied by atomic scale modelling

AbstractWe address the question of the sensitivity of the segregation and ordering properties to the surface energy of Au‐Pd alloys and nanoalloys. To this purpose, the parameterization of an EAM potential fitted to the vacancy formation energy was modified in order to predict surface energies closer to the experimentally measured ones. Metropolis Monte Carlo sampling was used to predict segregation and ordering. We show in the detail that modifying the surface energies of Au and Pd by up to 50 percent and, their surface energy difference by about 35 percent is of no major importance for the ordering and segregation properties of nanostructured Au‐Pd alloy surfaces (© 2010 WILEY‐VCH Verlag GmbH & Co. KGaA, Weinheim)

Methods for calculating the absolute entropy and free energy of biological systems based on ideas from polymer physics.

The commonly used simulation techniques, Metropolis Monte Carlo (MC) and molecular dynamics (MD) are of a dynamical type which enables one to sample system configurations i correctly with the Boltzmann probability, P(i)(B), while the value of P(i)(B) is not provided directly; therefore, it is difficult to obtain the absolute entropy, S approximately -ln P(i)(B), and the Helmholtz free energy, F. With a different simulation approach developed in polymer physics, a chain is grown step-by-step with transition probabilities (TPs), and thus their product is the value of the construction probability; therefore, the entropy is known. Because all exact simulation methods are equivalent, i.e. they lead to the same averages and fluctuations of physical properties, one can treat an MC or MD sample as if its members have rather been generated step-by-step. Thus, each configuration i of the sample can be reconstructed (from nothing) by calculating the TPs with which it could have been constructed. This idea applies also to bulk systems such as fluids or magnets. This approach has led earlier to the "local states" (LS) and the "hypothetical scanning" (HS) methods, which are approximate in nature. A recent development is the hypothetical scanning Monte Carlo (HSMC) (or molecular dynamics, HSMD) method which is based on stochastic TPs where all interactions are taken into account. In this respect, HSMC(D) can be viewed as exact and the only approximation involved is due to insufficient MC(MD) sampling for calculating the TPs. The validity of HSMC has been established by applying it first to liquid argon, TIP3P water, self-avoiding walks (SAW), and polyglycine models, where the results for F were found to agree with those obtained by other methods. Subsequently, HSMD was applied to mobile loops of the enzymes porcine pancreatic alpha-amylase and acetylcholinesterase in explicit water, where the difference in F between the bound and free states of the loop was calculated. Currently, HSMD is being extended for calculating the absolute and relative free energies of ligand-enzyme binding. We describe the whole approach and discuss future directions.

A multi-range order parameter for binary alloy bulk materials and nanoparticles

A multi-range order parameter for binary alloys is considered, valid for bulk materials and nanoparticles. The difference between two common definitions is discussed, arising when applied to small systems. Its properties are illustrated on the example of the order/disorder transitions in Au-Pd fcc bulk materials and nanoparticles. The thermodynamic equilibrium states of these systems are predicted by Metropolis Monte Carlo sampling in the canonical (NPT) ensemble and an embedded-atom model poten- tial. Short range order is found partially preserved by the transition while long range order vanishes. The order of intermediate ranges is found to be particularly sensitive to the presence of anti-phase boundaries (APBs) in both, bulk alloys and nanoparticles.

DIADORIM: a Monte Carlo Program for liquid simulations including quantum mechanics and molecular mechanics (QM/MM) facilities: applications to liquid ethanol

This paper presents a Metropolis Monte Carlo (MMC) computer program developed to simulate liquids and solutions including QM/MM facilities: the energy from intermolecular interactions is calculated with classical force field functions and the internal molecular energies are calculated using Quantum Chemistry methods. The following facilities were implemented: (i) the semiempirical MOPAC 6 quantum chemistry package was included as a subroutine of the main MMC simulation program; (ii) alternatively, an interface with an external data bank was developed to allow the use of energies previously obtained; (iii) calculation in NpT and NVT ensembles are available; (iv) the trajectory generated along the MMC sampling can be saved using standard xtc file format allowing trajectory visualization and data analysis using external programs. The program was used to calculate thermodynamical and structural properties of liquids ethanol (ET) in the NpT ensemble at 1.0 atm and 298 K. Bond angles and bond distances of ethanol molecules were kept constant but torsions along the dihedral angle defined by the H-O-C-C atoms were sampled in the simulation. QM/MM calculations were performed using the MOPAC AM1, PM3 and MNDO Hamiltonians to calculate the internal molecular energy. Calculations using internal energies obtained with OPLS-AA force field, ab initio HF/6-31g*, MP2/6-31g* and B3LYP/6-31g* calculations were performed. The dihedral angles distributions obtained with different methodologies reveal extraordinary qualitative and quantitative differences. However, the pair energies and radial distributions functions obtained with the different methodologies are in good agreement. The values calculated for the density and enthalpy of vaporization of liquid ethanol are in good agreement with experimental data.

Correlations between designability and various structural characteristics of protein lattice models

Using six kinds of lattice types (4 x 4, 5 x 5, and 6 x 6 square lattices; 3 x 3 x 3 cubic lattice; and 2+3+4+3+2 and 4+5+6+5+4 triangular lattices), three different size alphabets (HP, HNUP, and 20 letters), and two energy functions, the designability of protein structures is calculated based on random samplings of structures and common biased sampling (CBS) of protein sequence space. Then three quantities stability (average energy gap), foldability, and partnum of the structure, which are defined to elucidate the designability, are calculated. The authors find that whatever the type of lattice, alphabet size, and energy function used, there will be an emergence of highly designable (preferred) structure. For all cases considered, the local interactions reduce degeneracy and make the designability higher. The designability is sensitive to the lattice type, alphabet size, energy function, and sampling method of the sequence space. Compared with the random sampling method, both the CBS and the Metropolis Monte Carlo sampling methods make the designability higher. The correlation coefficients between the designability, stability, and foldability are mostly larger than 0.5, which demonstrate that they have strong correlation relationship. But the correlation relationship between the designability and the partnum is not so strong because the partnum is independent of the energy. The results are useful in practical use of the designability principle, such as to predict the protein tertiary structure.

Structure and Stability of Small Self-interstitial Clusters in Si: Prediction by First Principles-based Monte Carlo Simulations

We present an efficient computational method for determining the ground-state configurations of self-interstitial clusters in Si, based on a combination of continuous random network model, Metropolis Monte Carlo sampling, and Keating-like valence bond model. Using the new approach, we have determined the structure and stability of small interstitial clusters (In, 3 < n < 12). The results are in good agreement with first principles results.

Exhaustive Metropolis Monte Carlo sampling and analysis of polyalanine conformations adopted under the influence of hydrogen bonds

We propose a novel Metropolis Monte Carlo procedure for protein modeling and analyze the influence of hydrogen bonding on the distribution of polyalanine conformations. We use an atomistic model of the polyalanine chain with rigid and planar polypeptide bonds, and elastic alpha carbon valence geometry. We adopt a simplified energy function in which only hard-sphere repulsion and hydrogen bonding interactions between the atoms are considered. Our Metropolis Monte Carlo procedure utilizes local crankshaft moves and is combined with parallel tempering to exhaustively sample the conformations of 16-mer polyalanine. We confirm that Flory's isolated-pair hypothesis (the steric independence between the dihedral angles of individual amino acids) does not hold true in long polypeptide chains. In addition to 3(10)- and alpha-helices, we identify a kink stabilized by 2 hydrogen bonds with a shared acceptor as a common structural motif. Varying the strength of hydrogen bonds, we induce the helix-coil transition in the model polypeptide chain. We compare the propensities for various hydrogen bonding patterns and determine the degree of cooperativity of hydrogen bond formation in terms of the Hill coefficient. The observed helix-coil transition is also quantified according to Zimm-Bragg theory.

A Metropolis Monte Carlo Implementation of Bayesian Time-Domain Parameter Estimation: Application to Coupling Constant Estimation from Antiphase Multiplets

The Bayesian perspective on statistics asserts that it makes sense to speak of a probability of an unknown parameter having a particular value. Given a model for an observed, noise-corrupted signal, we may use Bayesian methods to estimate not only the most probable value for each parameter but also their distributions. We present an implementation of the Bayesian parameter estimation formalism developed by G. L. Bretthorst (1990,J. Magn. Reson.88, 533) using the Metropolis Monte Carlo sampling algorithm to perform the parameter and error estimation. This allows us to make very few assumptions about the shape of the posterior distribution, and allows the easy introduction of prior knowledge about constraints among the model parameters. We present evidence that the error estimates obtained in this manner are realistic, and that the Monte Carlo approach can be used to accurately estimate coupling constants from antiphase doublets in synthetic and experimental data.

“Mining Minima”: Direct Computation of Conformational Free Energy

We describe a novel, two-step method for directly computing the conformational free energy of a molecule. In the first step, a finite set of low-energy conformations is identified, and its contribution to the configuration integral is evaluated by a straightforward Monte Carlo technique. The method of finding energy minima incorporates certain features of the global-underestimator method and of a genetic algorithm. In the second step, the contribution to the configuration integral due to conformations not included in the initial integration is determined by Metropolis Monte Carlo sampling. Applications to alanine oligopeptides and to three cyclic urea inhibitors of HIV protease are presented.

DMC: A multifunctional hybrid dynamics/Monte Carlo simulation algorithm for the evaluation of conformational space

AbstractA hybrid conformational search algorithm (DMC) is described that combines a modified form of molecular dynamics with Metropolis Monte Carlo sampling, using the COSMIC(90) force field. Trial configurations are generated by short bursts of high‐temperature dynamics in which the initial kinetic energy is focused into single bond rotations or alternatively into “corner‐flapping” motions in ring systems. Constant temperature and simulated annealing search protocols have been applied to the conformational analysis of several model hydrocarbons (cyclopentane, cyclohexane, cycloheptane, cyclooctane, cycloheptadecane, decane, and tetradecane), and the performance compared with conventional molecular dynamics and Monte Carlo sampling methods. Optimum Metropolis sampling temperatures have been determined and range from 1000–2000 K for acyclic molecules to 3000 K for cyclic systems. Simulated annealing runs are most successful at locating the global minimum when cooling slowing from these optimum temperatures.

Supercoiled DNA energetics and dynamics by computer simulation

A new formulation is presented for investigating supercoiled DNA configurations by deterministic techniques. Thus far, the computational difficulties involved in applying deterministic methods to supercoiled DNA studies have generally limited computer simulations to stochastic approaches. While stochastic methods, such as simulated annealing and Metropolis-Monte Carlo sampling, are successful at generating a large number of configurations and estimating thermodynamic properties of topoisomer ensembles, deterministic methods offer an accurate characterization of the minima and a systematic following of their dynamics. To make this feasible, we model circular duplex DNA compactly by a B-spline ribbon-like model in terms of a small number of control vertices. We associate an elastic deformation energy composed of bending and twisting integrals and represent intrachain contact by a 6–12 Lennard Jones potential. The latter is parameterized to yield an energy minimum at the observed DNA-helix diameter inclusive of a hydration shell. A penalty term to ensure fixed contour length is also included. First and second partial derivatives of the energy function have been derived by using various mathematical simplifications. First derivatives are essential for Newton-type minimization as well as molecular dynamics, and partial second-derivative information can significantly accelerate minimization convergence through preconditioning. Here we apply a new large-scale truncated-Newton algorithm for minimization and a Langevin/implicit-Euler scheme for molecular dynamics. Our truncated-Newton method exploits the separability of potential energy functions into terms of differing complexity. It relies on a preconditioned conjugate gradient method that is efficient for large-scale problems to solve approximately for the search direction at every step. Our dynamics algorithm is numerically stable over large time steps. It also introduces a frequency-discriminating mechanism so that vibrational modes with frequencies greater than a chosen cutoff frequency are essentially frozen by the method. With these tools, we rapidly identify corresponding circular and interwound energy minima for small DNA rings for a series of imposed linking-number differences. These structures are consistent with available electron microscopy data. The energetic exchange of stability between the circle and the figure-8, in very good agreement with analytical results, is also detailed. Molecular dynamics trajectories at 100 femtosecond time steps then reveal the rapid folding of the unstable circular state into supercoiled forms. Significant bending and twisting motions of the interwound structures are also observed. Such information may be useful for understanding transition states along the folding pathway and the role of enzymes that regulate supercoiling. More generally, new quantitative data obtained by such deterministic approaches may help in interpreting the effects of supercoiling on key biological processes.

Monte Carlo-minimization approach to the multiple-minima problem in protein folding.

A Monte Carlo-minimization method has been developed to overcome the multiple-minima problem. The Metropolis Monte Carlo sampling, assisted by energy minimization, surmounts intervening barriers in moving through successive discrete local minima in the multidimensional energy surface. The method has located the lowest-energy minimum thus far reported for the brain pentapeptide [Met5]enkephalin in the absence of water. Presumably it is the global minimum-energy structure. This supports the concept that protein folding may be a Markov process. In the presence of water, the molecules appear to exist as an ensemble of different conformations.