On the Convergence of Parallel Tempering Monte Carlo Simulations of LJ38

100 Years of the Lennard-Jones Potential

It is now 100 years since Lennard-Jones published his first paper introducing the now famous potential that bears his name. It is therefore timely to reflect on the many achievements, as well as the limitations, of this potential in the theory of atomic and molecular interactions, where applications range from descriptions of intermolecular forces to molecules, clusters, and condensed matter.

Analysis and interpretation of first passage time distributions featuring rare events

In this contribution we consider theory and associated computational tools to treat the kinetics associated with competing pathways on multifunnel energy landscapes. Multifunnel landscapes are associated with molecular switches and multifunctional materials, and are expected to exhibit multiple relaxation time scales and associated thermodynamic signatures in the heat capacity. Our focus here is on the first passage time distribution, which is encoded in a kinetic transition network containing all the locally stable states and the pathways between them. This network can be renormalised to reduce the dimensionality, while exactly conserving the mean first passage time and approximately conserving the full distribution. The structure of the reduced network can be visualised using disconnectivity graphs. We show how features in the first passage time distribution can be associated with specific kinetic traps, and how the appearance of competing relaxation time scales depends on the starting conditions. The theory is tested for two model landscapes and applied to an atomic cluster and a disordered peptide. Our most important contribution is probably the reconstruction of the full distribution for long time scales, where numerical problems prevent direct calculations. Here we combine accurate treatment of the mean first passage time with the reliable part of the distribution corresponding to faster time scales. Hence we now have a fundamental understanding of both thermodynamic and kinetic signatures of multifunnel landscapes.

Combining experiment and energy landscapes to explore anaerobic heme breakdown in multifunctional hemoproteins

To survive, many pathogens extract heme from their host organism and break down the porphyrin scaffold to sequester the Fe2+ ion via a heme oxygenase. Recent studies have revealed that certain pathogens can anaerobically degrade heme. Our own research has shown that one such pathway proceeds via NADH-dependent heme degradation, which has been identified in a family of hemoproteins from a range of bacteria. HemS, from Yersinia enterocolitica, is the main focus of this work, along with HmuS (Yersinia pestis), ChuS (Escherichia coli) and ShuS (Shigella dysenteriae). We combine experiments, Energy Landscape Theory, and a bioinformatic investigation to place these homologues within a wider phylogenetic context. A subset of these hemoproteins are known to bind certain DNA promoter regions, suggesting not only that they can catalytically degrade heme, but that they are also involved in transcriptional modulation responding to heme flux. Many of the bacterial species responsible for these hemoproteins (including those that produce HemS, ChuS and ShuS) are known to specifically target oxygen-depleted regions of the gastrointestinal tract. A deeper understanding of anaerobic heme breakdown processes exploited by these pathogens could therefore prove useful in the development of future strategies for disease prevention.

Rare events and first passage time statistics from the energy landscape

We analyze the probability distribution of rare first passage times corresponding to transitions between product and reactant states in a kinetic transition network. The mean first passage times and the corresponding rate constants are analyzed in detail for two model landscapes and the double funnel landscape corresponding to an atomic cluster. Evaluation schemes based on eigendecomposition and kinetic path sampling, which both allow access to the first passage time distribution, are benchmarked against mean first passage times calculated using graph transformation. Numerical precision issues severely limit the useful temperature range for eigendecomposition, but kinetic path sampling is capable of extending the first passage time analysis to lower temperatures, where the kinetics of interest constitute rare events. We then investigate the influence of free energy based state regrouping schemes for the underlying network. Alternative formulations of the effective transition rates for a given regrouping are compared in detail to determine their numerical stability and capability to reproduce the true kinetics, including recent coarse-graining approaches that preserve occupancy cross correlation functions. We find that appropriate regrouping of states under the simplest local equilibrium approximation can provide reduced transition networks with useful accuracy at somewhat lower temperatures. Finally, a method is provided to systematically interpolate between the local equilibrium approximation and exact intergroup dynamics. Spectral analysis is applied to each grouping of states, employing a moment-based mode selection criterion to produce a reduced state space, which does not require any spectral gap to exist, but reduces to gap-based coarse graining as a special case. Implementations of the developed methods are freely available online.

Exploring Energy Landscapes

Recent advances in the potential energy landscapes approach are highlighted, including both theoretical and computational contributions. Treating the high dimensionality of molecular and condensed matter systems of contemporary interest is important for understanding how emergent properties are encoded in the landscape and for calculating these properties while faithfully representing barriers between different morphologies. The pathways characterized in full dimensionality, which are used to construct kinetic transition networks, may prove useful in guiding such calculations. The energy landscape perspective has also produced new procedures for structure prediction and analysis of thermodynamic properties. Basin-hopping global optimization, with alternative acceptance criteria and generalizations to multiple metric spaces, has been used to treat systems ranging from biomolecules to nanoalloy clusters and condensed matter. This review also illustrates how all this methodology, developed in the context of chemical physics, can be transferred to landscapes defined by cost functions associated with machine learning.

Decoding heat capacity features from the energy landscape

A general scheme is derived to connect transitions in configuration space with features in the heat capacity. A formulation in terms of occupation probabilities for local minima that define the potential energy landscape provides a quantitative description of how contributions arise from competition between different states. The theory does not rely on a structural interpretation for the local minima, so it is equally applicable to molecular energy landscapes and the landscapes defined by abstract functions. Applications are presented for low-temperature solid-solid transitions in atomic clusters, which involve just a few local minima with different morphologies, and for cluster melting, which is driven by the landscape entropy associated with the more numerous high-energy minima. Analyzing these features in terms of the balance between states with increasing and decreasing occupation probabilities provides a direct interpretation of the underlying transitions. This approach enables us to identify a qualitatively different transition that is caused by a single local minimum associated with an exceptionally large catchment volume in configuration space for a machine learning landscape.

Grand and Semigrand Canonical Basin-Hopping

We introduce grand and semigrand canonical global optimization approaches using basin-hopping with an acceptance criterion based on the local contribution of each potential energy minimum to the (semi)grand potential. The method is tested using local harmonic vibrational densities of states for atomic clusters as a function of temperature and chemical potential. The predicted global minima switch from dissociated states to clusters for larger values of the chemical potential and lower temperatures, in agreement with the predictions of a model fitted to heat capacity data for selected clusters. Semigrand canonical optimization allows us to identify particularly stable compositions in multicomponent nanoalloys as a function of increasing temperature, whereas the grand canonical potential can produce a useful survey of favorable structures as a byproduct of the global optimization search.

Communication: optimal parameters for basin-hopping global optimization based on Tsallis statistics

A fundamental problem associated with global optimization is the large free energy barrier for the corresponding solid-solid phase transitions for systems with multi-funnel energy landscapes. To address this issue we consider the Tsallis weight instead of the Boltzmann weight to define the acceptance ratio for basin-hopping global optimization. Benchmarks for atomic clusters show that using the optimal Tsallis weight can improve the efficiency by roughly a factor of two. We present a theory that connects the optimal parameters for the Tsallis weighting, and demonstrate that the predictions are verified for each of the test cases.

Replica exchange statistical temperature molecular dynamics algorithm

The replica exchange statistical temperature molecular dynamics (RESTMD) algorithm is presented, designed to alleviate an extensive increase of the number of replicas required as system size increases in the conventional temperature replica exchange method (tREM), and to obtain improved sampling in individual replicas. RESTMD optimally integrates multiple STMD (Phys. Rev. Lett. 2006, 97, 050601) runs with replica exchanges, giving rise to a flat energy sampling in each replica with a self-adjusting weight determination. The expanded flat energy dynamic sampling range allows the use of significantly fewer STMD replicas while maintaining the desired acceptance probability for replica exchanges. The computational advantages of RESTMD over conventional REM and single-replica STMD are explicitly demonstrated with an application to a coarse-grained protein model. The effect of two different kinetic temperature control schemes on the sampling efficiency is explored for diverse simulation conditions.

An infinite swapping approach to the rare-event sampling problem

We describe a new approach to the rare-event Monte Carlo sampling problem. This technique utilizes a symmetrization strategy to create probability distributions that are more highly connected and, thus, more easily sampled than their original, potentially sparse counterparts. After discussing the formal outline of the approach and devising techniques for its practical implementation, we illustrate the utility of the technique with a series of numerical applications to Lennard-Jones clusters of varying complexity and rare-event character.

Spatially local parallel tempering for thermal-equilibrium sampling

Parallel tempering (PT) is a set of techniques for accelerating thermal-equilibrium sampling in systems where the exploration of configuration space is hindered by energy barriers. With standard PT algorithms, the computational effort scales unfavorably with system size, so that it is difficult to apply them to large systems. We propose local PT algorithms, for which the computational effort is proportional to the number of degrees of freedom. We demonstrate the effectiveness of the new algorithms on two one-dimensional model systems, showing that results for selected observables are correctly reproduced, and that practical linear scaling is achieved. We show also that the algorithms are readily applied to systems in higher dimensions. We note the prospects for studying large extended systems, including surfaces and interfaces.

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.

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.

Global optimization method for cluster structures

A time-going-backward quasidynamics method is developed for global optimization of cluster structures, and its merits are examined by a simple classical mechanics model, indicating that the probability for the system to jump over high potential barriers by this method is much higher than that by common annealing methods. The method is then used to investigate the isomers of a Lennard-Jones cluster containing 38 atoms and the C60 cluster with the Brenner potential, and can easily give the most stable structures, which are difficult to obtain by common annealing methods. In addition, for small carbon clusters C_{n} (n=21-30) , most of the potential energies optimized by this method are much lower than those obtained by a genetic algorithm.

Equilibrium density of states and thermodynamic properties of a model glass former

This paper presents an analysis of the thermodynamics of a model glass former. We have performed equilibrium sampling of a popular binary Lennard-Jones model, employing parallel tempering Monte Carlo to cover the crystalline, amorphous, and liquid regions of configuration space. Disconnectivity graphs are used to visualize the potential energy landscape in the vicinity of a crystalline geometry and in an amorphous region of configuration space. The crystalline global minimum is separated from the bulk of the minima by a large potential energy gap, leading to broken ergodicity in conventional simulations. Our sampling reveals crystalline global minima that are lower in potential energy than some of the previous candidates. We present equilibrium thermodynamic properties based on parallel tempering simulations, including heat capacities and free energy profiles, which depend explicitly on the crystal structure. We also report equilibrium melting temperatures.

Multiple structural transformations in Lennard-Jones clusters: Generic versus size-specific behavior

The size-temperature "phase diagram" for Lennard-Jones clusters LJn with sizes up to n=147 is constructed based on the analysis of the heat capacities and orientational bond order parameter distributions computed by the exchange Monte Carlo method. Two distinct types of "phase transitions" accompanied by peaks in the heat capacities are proven to be generic. Clusters with Mackay atom packing in the overlayer undergo a lower-temperature melting (or Mackay-anti-Mackay) transition that occurs within the overlayer. All clusters undergo a higher-temperature transition, which for the three-layer clusters is proven to be the 55-atom-core-melting transition. For the two-layer clusters, the core/overlayer subdivision is ambiguous, so the higher-temperature transition is better characterized as the breaking of the local icosahedral coordination symmetry. A pronounced size-specific behavior can typically be observed at low temperatures and often occurs in clusters with highly symmetric global minima. An example of such behavior is LJ135, which undergoes a low-temperature solid-solid transition, besides the two generic transitions, i.e., the overlayer reconstruction and the core melting.

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.

Structural transitions in the 309-atom magic number Lennard-Jones cluster

The thermal behavior of the 309-atom Lennard-Jones cluster, whose structure is a complete Mackay icosahedron, has been studied by parallel tempering Monte Carlo simulations. Surprisingly for a magic number cluster, the heat capacity shows a very pronounced peak before melting, which is attributed to several coincident structural transformation processes. The main transformation is somewhat akin to surface roughening and involves a cooperative condensation of vacancies and adatoms that leads to the formation of pits and islands one or two layers thick on the Mackay icosahedron. The second transition in order of importance involves a whole scale transformation of the cluster structure and leads to a diverse set of twinned structures that are assemblies of face-centered-cubic tetrahedra with six atoms along their edges, i.e., one atom more than the edges of the 20 tetrahedra that make up the 309-atom Mackay icosahedron. A surface reconstruction of the icosahedron from a Mackay to an anti-Mackay overlayer is also observed, but with a lower probability.

Equilibrium thermodynamics from basin-sampling

We present a "basin-sampling" approach for calculation of the potential energy density of states for classical statistical models. It combines a Wang-Landau-type uniform sampling of local minima and a novel approach for approximating the relative contributions from local minima in terms of the volumes of basins of attraction. We have employed basin-sampling to study phase changes in atomic clusters modeled by the Lennard-Jones potential and for ionic clusters. The approach proves to be efficient for systems involving broken ergodicity and has allowed us to calculate converged heat capacity curves for systems that could previously only be treated using the harmonic superposition approximation. Benchmarks are also provided by comparison with parallel tempering and Wang-Landau simulations, where these proved feasible.

Structural Transitions and Melting in LJ74-78 Lennard-Jones Clusters from Adaptive Exchange Monte Carlo Simulations

Phase changes in Lennard-Jones (LJ) clusters containing between 74 and 78 atoms are investigated by means of exchange Monte Carlo simulations in the canonical ensemble. The replica temperatures are self-adapted to facilitate the convergence. Although the 74- and 78-atom clusters have icosahedral global minima, the clusters with 75-77 atoms have decahedral ground-state structures and they undergo a structural transition to icosahedral minima before melting. The structural transitions are characterized by quenching and by looking at the Q4 and Q6 orientational bond order parameters. The transition temperatures are estimated to be 0.114, 0.065, and 0.074 reduced units for LJ75, LJ76, and LJ77, respectively. These values, their ordering and the associated latent heats are compared with other estimates based on the harmonic superposition approach.

All-exchanges parallel tempering

An alternative exchange strategy for parallel tempering simulations is introduced. Instead of attempting to swap configurations between two randomly chosen but adjacent replicas, the acceptance probabilities of all possible swap moves are calculated a priori. One specific swap move is then selected according to its probability and enforced. The efficiency of the method is illustrated first on the case of two Lennard-Jones (LJ) clusters containing 13 and 31 atoms, respectively. The convergence of the caloric curve is seen to be at least twice as fast as in conventional parallel tempering simulations, especially for the difficult case of LJ31. Further evidence for an improved efficiency is reported on the ergodic measure introduced by Mountain and Thirumalai [J. Phys. Chem. 93, 6975 (1989)], calculated here for LJ13 close to the melting point. Finally, tests on two simple spin systems indicate that the method should be particularly useful when a limited number of replicas are available.