Generalized simulated annealing algorithms using Tsallis statistics: Application to conformational optimization of a tetrapeptide

Entropy defect in thermodynamics

This paper describes the physical foundations of the newly discovered "entropy defect" as a basic concept of thermodynamics. The entropy defect quantifies the change in entropy caused by the order induced in a system through the additional correlations among its constituents when two or more subsystems are assembled. This defect is closely analogous to the mass defect that arises when nuclear particle systems are assembled. The entropy defect determines how the entropy of the system compares to its constituent's entropies and stands on three fundamental properties: each constituent's entropy must be (i) separable, (ii) symmetric, and (iii) bounded. We show that these properties provide a solid foundation for the entropy defect and for generalizing thermodynamics to describe systems residing out of the classical thermal equilibrium, both in stationary and nonstationary states. In stationary states, the consequent thermodynamics generalizes the classical framework, which was based on the Boltzmann-Gibbs entropy and Maxwell-Boltzmann canonical distribution of particle velocities, into the respective entropy and canonical distribution associated with kappa distributions. In nonstationary states, the entropy defect similarly acts as a negative feedback, or reduction of the increase of entropy, preventing its unbounded growth toward infinity.

Generation of Quantum Configurational Ensembles Using Approximate Potentials

Conformational analysis is of paramount importance in drug design: it is crucial to determine pharmacological properties, understand molecular recognition processes, and characterize the conformations of ligands when unbound. Molecular Mechanics (MM) simulation methods, such as Monte Carlo (MC) and molecular dynamics (MD), are usually employed to generate ensembles of structures due to their ability to extensively sample the conformational space of molecules. The accuracy of these MM-based schemes strongly depends on the functional form of the force field (FF) and its parametrization, components that often hinder their performance. High-level methods, such as ab initio MD, provide reliable structural information but are still too computationally expensive to allow for extensive sampling. Therefore, to overcome these limitations, we present a multilevel MC method that is capable of generating quantum configurational ensembles while keeping the computational cost at a minimum. We show that FF reparametrization is an efficient route to generate FFs that reproduce QM results more closely, which, in turn, can be used as low-cost models to achieve the gold standard QM accuracy. We demonstrate that the MC acceptance rate is strongly correlated with various phase space overlap measurements and that it constitutes a robust metric to evaluate the similarity between the MM and QM levels of theory. As a more advanced application, we present a self-parametrizing version of the algorithm, which combines sampling and FF parametrization in one scheme, and apply the methodology to generate the QM/MM distribution of a ligand in aqueous solution.

Quantum Chemistry Calculations for Metabolomics

A primary goal of metabolomics studies is to fully characterize the small-molecule composition of complex biological and environmental samples. However, despite advances in analytical technologies over the past two decades, the majority of small molecules in complex samples are not readily identifiable due to the immense structural and chemical diversity present within the metabolome. Current gold-standard identification methods rely on reference libraries built using authentic chemical materials ("standards"), which are not available for most molecules. Computational quantum chemistry methods, which can be used to calculate chemical properties that are then measured by analytical platforms, offer an alternative route for building reference libraries, i.e., in silico libraries for "standards-free" identification. In this review, we cover the major roadblocks currently facing metabolomics and discuss applications where quantum chemistry calculations offer a solution. Several successful examples for nuclear magnetic resonance spectroscopy, ion mobility spectrometry, infrared spectroscopy, and mass spectrometry methods are reviewed. Finally, we consider current best practices, sources of error, and provide an outlook for quantum chemistry calculations in metabolomics studies. We expect this review will inspire researchers in the field of small-molecule identification to accelerate adoption of in silico methods for generation of reference libraries and to add quantum chemistry calculations as another tool at their disposal to characterize complex samples.

Biomolecular QM/MM Simulations: What Are Some of the "Burning Issues"?

QM/MM simulations have become an indispensable tool in many chemical and biochemical investigations. Considering the tremendous degree of success, including recognition by a 2013 Nobel Prize in Chemistry, are there still "burning challenges" in QM/MM methods, especially for biomolecular systems? In this short Perspective, we discuss several issues that we believe greatly impact the robustness and quantitative applicability of QM/MM simulations to many, if not all, biomolecules. We highlight these issues with observations and relevant advances from recent studies in our group and others in the field. Despite such limited scope, we hope the discussions are of general interest and will stimulate additional developments that help push the field forward in meaningful directions.

Multi-level free energy simulation with a staged transformation approach

Combining multiple levels of theory in free energy simulations to balance computational accuracy and efficiency is a promising approach for studying processes in the condensed phase. While the basic idea has been proposed and explored for quite some time, it remains challenging to achieve convergence for such multi-level free energy simulations as it requires a favorable distribution overlap between different levels of theory. Previous efforts focused on improving the distribution overlap by either altering the low-level of theory for the specific system of interest or ignoring certain degrees of freedom. Here, we propose an alternative strategy that first identifies the degrees of freedom that lead to gaps in the distributions of different levels of theory and then treats them separately with either constraints or restraints or by introducing an intermediate model that better connects the low and high levels of theory. As a result, the conversion from the low level to the high level model is done in a staged fashion that ensures a favorable distribution overlap along the way. Free energy components associated with different steps are mostly evaluated explicitly, and thus, the final result can be meaningfully compared to the rigorous free energy difference between the two levels of theory with limited and well-defined approximations. The additional free energy component calculations involve simulations at the low level of theory and therefore do not incur high computational costs. The approach is illustrated with two simple but non-trivial solution examples, and factors that dictate the reliability of the result are discussed.

Nonextensive Statistical Mechanics: Equivalence Between Dual Entropy and Dual Probabilities

The concept of duality of probability distributions constitutes a fundamental “brick” in the solid framework of nonextensive statistical mechanics—the generalization of Boltzmann–Gibbs statistical mechanics under the consideration of the q-entropy. The probability duality is solving old-standing issues of the theory, e.g., it ascertains the additivity for the internal energy given the additivity in the energy of microstates. However, it is a rather complex part of the theory, and certainly, it cannot be trivially explained along the Gibb’s path of entropy maximization. Recently, it was shown that an alternative picture exists, considering a dual entropy, instead of a dual probability. In particular, the framework of nonextensive statistical mechanics can be equivalently developed using q- and 1/q- entropies. The canonical probability distribution coincides again with the known q-exponential distribution, but without the necessity of the duality of ordinary-escort probabilities. Furthermore, it is shown that the dual entropies, q-entropy and 1/q-entropy, as well as, the 1-entropy, are involved in an identity, useful in theoretical development and applications.

The Structure of Adamantane Clusters: Atomistic vs. Coarse-Grained Predictions From Global Optimization

Candidate structures for the global minima of adamantane clusters, (C₁₀H₁₆)_N, are presented. Based on a rigid model for individual molecules with atom-atom pairwise interactions that include Lennard-Jones and Coulomb contributions, low-energy structures were obtained up to N = 42 using the basin-hopping method. The results indicate that adamantane clusters initially grow accordingly with an icosahedral packing scheme, followed above N = 14 by a structural transition toward face-centered cubic structures. The special stabilities obtained at N = 13, 19, and 38 are consistent with these two structural families, and agree with recent mass spectrometry measurements on cationic adamantane clusters. Coarse-graining the intermolecular potential by averaging over all possible orientations only partially confirm the all-atom results, the magic numbers at 13 and 38 being preserved. However, the details near the structural transition are not captured well, because despite their high symmetry the adamantane molecules are still rather anisotropic.

Enhanced sampling method in molecular simulations using genetic algorithm for biomolecular systems

We propose a molecular simulation method using genetic algorithm (GA) for biomolecular systems to obtain ensemble averages efficiently. In this method, we incorporate the genetic crossover, which is one of the operations of GA, to any simulation method such as conventional molecular dynamics (MD), Monte Carlo, and other simulation methods. The genetic crossover proposes candidate conformations by exchanging parts of conformations of a target molecule between a pair of conformations during the simulation. If the candidate conformations are accepted, the simulation resumes from the accepted ones. While conventional simulations are based on local update of conformations, the genetic crossover introduces global update of conformations. As an example of the present approach, we incorporated genetic crossover to MD simulations. We tested the validity of the method by calculating ensemble averages and the sampling efficiency by using two kinds of peptides, ALA3 and (AAQAA)₃ . The results show that for ALA3 system, the distribution probabilities of backbone dihedral angles are in good agreement with those of the conventional MD and replica-exchange MD simulations. In the case of (AAQAA)₃ system, our method showed lower structural correlation of α-helix structures than the other two methods and more flexibility in the backbone ψ angles than the conventional MD simulation. These results suggest that our method gives more efficient conformational sampling than conventional simulation methods based on local update of conformations. © 2018 Wiley Periodicals, Inc.

Analysis of Temperature andq-Parameter Dependency of FCM with Tsallis Entropy Maximization

The Tsallis entropy is aq-parameter extension of the Shannon entropy. By maximizing it within the framework of fuzzyc-means, statistical mechanical membership functions can be derived. We propose a clustering algorithm that includes the membership function and deterministic annealing. One of the major issues for this method is the determination of an appropriate values forqand an initial annealing temperature for a given data distribution. Accordingly, in our previous study, we investigated the relationship betweenqand the annealing temperature. We quantitatively compared the area of the membership function for various values ofqand for various temperatures. The results showed that the effect ofqon the area was nearly the inverse of that of the temperature. In this paper, we analytically investigate this relationship by directly integrating the membership function, and the inversely proportional relationship betweenqand the temperature is approximately confirmed. Based on this relationship, aq-incrementation deterministic annealing fuzzyc-means (FCM) algorithm is developed. Experiments are performed, and it is confirmed that the algorithm works properly. However, it is also confirmed that differences in the shape of the membership function of the annealing method and that of theq-incrementation method are remained.

Conformational sampling of metastable states: Tq-REM as a novel replica exchange method

Although the replica exchange methods (REMs) were developed as efficient conformational sampling methods for bio-molecular simulations, their application to very large bio-systems is somewhat limited. We propose a new replica exchange scheme (Tq-REM) created by combining the conventional temperature-REM (T-REM) and one of the Hamiltonian-REMs, q-REM, using the effective potential with reduced barriers. In the proposed Tq-REM scheme, high temperature replicas in T-REM are substituted with q-replicas. This combined scheme is expected to exploit advantages of the T-REM and q-REM resulting in improved sampling efficiency while minimizing the drawbacks of both approaches. We investigated the performance of Tq-REM compared with T-REM by performing all-atom MD simulations on Met-enkephalin, (AAQAA)₃, and Trpzip2. It was found that convergence of the free energy surfaces was improved by Tq-REM over the conventional T-REM. In particular, the trajectories of Tq-REM were able to sample the relevant conformations for all of the metastable folding intermediates, while some of the local minimum structures are poorly represented by T-REM. The results of the present study suggest that Tq-REM can provide useful tools to investigate systems where metastable states play important roles.

Enhancing the branching ratios in the dissociation channels for O(16)O(16)O(18) molecule by designing optimum laser pulses: A study using stochastic optimization

We propose a strategy of using a stochastic optimization technique, namely, simulated annealing to design optimum laser pulses (both IR and UV) to achieve greater fluxes along the two dissociating channels (O(18) + O(16)O(16) and O(16) + O(16)O(18)) in O(16)O(16)O(18) molecule. We show that the integrated fluxes obtained along the targeted dissociating channel is larger with the optimized pulse than with the unoptimized one. The flux ratios are also more impressive with the optimized pulse than with the unoptimized one. We also look at the evolution contours of the wavefunctions along the two channels with time after the actions of both the IR and UV pulses and compare the profiles for unoptimized (initial) and optimized fields for better understanding the results that we achieve. We also report the pulse parameters obtained as well as the final shapes they take.

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.

Generalized essential energy space random walks to more effectively accelerate solute sampling in aqueous environment

Molecular dynamics sampling can be enhanced via the promoting of potential energy fluctuations, for instance, based on a Hamiltonian modified with the addition of a potential-energy-dependent biasing term. To overcome the diffusion sampling issue, which reveals the fact that enlargement of event-irrelevant energy fluctuations may abolish sampling efficiency, the essential energy space random walk (EESRW) approach was proposed earlier. To more effectively accelerate the sampling of solute conformations in aqueous environment, in the current work, we generalized the EESRW method to a two-dimension-EESRW (2D-EESRW) strategy. Specifically, the essential internal energy component of a focused region and the essential interaction energy component between the focused region and the environmental region are employed to define the two-dimensional essential energy space. This proposal is motivated by the general observation that in different conformational events, the two essential energy components have distinctive interplays. Model studies on the alanine dipeptide and the aspartate-arginine peptide demonstrate sampling improvement over the original one-dimension-EESRW strategy; with the same biasing level, the present generalization allows more effective acceleration of the sampling of conformational transitions in aqueous solution. The 2D-EESRW generalization is readily extended to higher dimension schemes and employed in more advanced enhanced-sampling schemes, such as the recent orthogonal space random walk method.

TUNING STRATEGIES IN CONSTRAINED SIMULATED ANNEALING FOR NONLINEAR GLOBAL OPTIMIZATION

This paper studies various strategies in constrained simulated annealing (CSA), a global optimization algorithm that achieves asymptotic convergence to constrained global minima (CGM) with probability one for solving discrete constrained nonlinear programming problems (NLPs). The algorithm is based on the necessary and sufficient condition for discrete constrained local minima (CLM) in the theory of discrete Lagrange multipliers and its extensions to continuous and mixed-integer constrained NLPs. The strategies studied include adaptive neighborhoods, distributions to control sampling, acceptance probabilities, and cooling schedules. We report much better solutions than the best-known solutions in the literature on two sets of continuous benchmarks and their discretized versions.

First-step mutations for adaptation at elevated temperature increase capsid stability in a virus

The relationship between mutation, protein stability and protein function plays a central role in molecular evolution. Mutations tend to be destabilizing, including those that would confer novel functions such as host-switching or antibiotic resistance. Elevated temperature may play an important role in preadapting a protein for such novel functions by selecting for stabilizing mutations. In this study, we test the stability change conferred by single mutations that arise in a G4-like bacteriophage adapting to elevated temperature. The vast majority of these mutations map to interfaces between viral coat proteins, suggesting they affect protein-protein interactions. We assess their effects by estimating thermodynamic stability using molecular dynamic simulations and measuring kinetic stability using experimental decay assays. The results indicate that most, though not all, of the observed mutations are stabilizing.

Simulated Q-annealing: conformational search with an effective potential

We have tested a version of the generalized simulated annealing algorithm based on molecular dynamics simulations with effective potential suggested by Tsallis statistics. The generalized annealing method, termed "simulated Q-annealing (SQ)," is applied to the simulations of a synthetic 11-residue peptide segment (1AQG). In SQ, the energy barriers between local minima change as the parameter q is varied and specific degrees of freedom can be selectively heated up and annealed. Conformational dynamics obtained by ordinary simulated annealing (SA) and SQ simulations are compared in order to illustrate the effectiveness of the SQ approach in conformational searching. We show that SQ can navigate the potential energy surface efficiently with a simple annealing protocol and demonstrate that conformations sampled by SQ can represent the funnel-like free energy surface.

Generalized simulated tempering for exploring strong phase transitions

An extension of the simulation tempering algorithm is proposed. It is shown to be particularly suited to the exploration of first-order phase transition systems characterized by the backbending or S-loop in the statistical temperature or a microcanonical caloric curve. A guided Markov process in an auxiliary parameter space systematically combines a set of parametrized Tsallis-weight ensemble simulations, which are targeted to transform unstable or metastable energy states of canonical ensembles into stable ones and smoothly join ordered and disordered phases across phase transition regions via a succession of unimodal energy distributions. The inverse mapping between the sampling weight and the effective temperature enables an optimal selection of relevant Tsallis-weight parameters. A semianalytic expression for the biasing weight in parameter space is adaptively updated "on the fly" during the simulation to achieve rapid convergence. Accelerated tunneling transitions with a comprehensive sampling for phase-coexistent states are explicitly demonstrated in systems subject to strong hysteresis including Potts and Ising spin models and a 147 atom Lennard-Jones cluster.

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.

CHARMM: the biomolecular simulation program

CHARMM (Chemistry at HARvard Molecular Mechanics) is a highly versatile and widely used molecular simulation program. It has been developed over the last three decades with a primary focus on molecules of biological interest, including proteins, peptides, lipids, nucleic acids, carbohydrates, and small molecule ligands, as they occur in solution, crystals, and membrane environments. For the study of such systems, the program provides a large suite of computational tools that include numerous conformational and path sampling methods, free energy estimators, molecular minimization, dynamics, and analysis techniques, and model-building capabilities. The CHARMM program is applicable to problems involving a much broader class of many-particle systems. Calculations with CHARMM can be performed using a number of different energy functions and models, from mixed quantum mechanical-molecular mechanical force fields, to all-atom classical potential energy functions with explicit solvent and various boundary conditions, to implicit solvent and membrane models. The program has been ported to numerous platforms in both serial and parallel architectures. This article provides an overview of the program as it exists today with an emphasis on developments since the publication of the original CHARMM article in 1983.

Coalescence phenomena in 1D silver nanostructures

Different coalescence processes on 1D silver nanostructures synthesized by a PVP assisted reaction in ethylene glycol at 160 °C were studied experimentally and theoretically. Analysis by TEM and HRTEM shows different defects found on the body of these materials, suggesting that they were induced by previous coalescence processes in the synthesis stage. TEM observations showed that irradiation with the electron beam eliminates the boundaries formed near the edges of the structures, suggesting that this process can be carried out by the application of other means of energy (i.e. thermal). These results were also confirmed by theoretical calculations by Monte Carlo simulations using a Sutton-Chen potential. A theoretical study by molecular dynamics simulation of the different coalescence processes on 1D silver nanostructures is presented, showing a surface energy driven sequence followed to form the final coalesced structure. Calculations were made at 1000-1300 K, which is near the melting temperature of silver (1234 K). Based on these results, it is proposed that 1D nanostructures can grow through a secondary mechanism based on coalescence, without losing their dimensionality.

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.

Two-time Green's functions and spectral density method in nonextensive quantum statistical mechanics

We extend the formalism of the thermodynamic two-time Green's functions to nonextensive quantum statistical mechanics. Working in the optimal Lagrangian multiplier representation, the q -spectral properties and the methods for a direct calculation of the two-time q Green's functions and the related q -spectral density ( q measures the nonextensivity degree) for two generic operators are presented in strict analogy with the extensive (q=1) counterpart. Some emphasis is devoted to the nonextensive version of the less known spectral density method whose effectiveness in exploring equilibrium and transport properties of a wide variety of systems has been well established in conventional classical and quantum many-body physics. To check how both the equations of motion and the spectral density methods work to study the q -induced nonextensivity effects in nontrivial many-body problems, we focus on the equilibrium properties of a second-quantized model for a high-density Bose gas with strong attraction between particles for which exact results exist in extensive conditions. Remarkably, the contributions to several thermodynamic quantities of the q -induced nonextensivity close to the extensive regime are explicitly calculated in the low-temperature regime by overcoming the calculation of the q grand-partition function.

Myosin dynamics on the millisecond time scale

Myosin is a motor protein associating with actin and ATP. It translates along actin filaments against a force by transduction of free energy liberated with ATP hydrolysis. Various myosin crystal structures define time points during ATPase showing the protein undergoes large conformation change during transduction over a cycle with approximately 10 ms periodicity. The protein conformation trajectory between two intermediates in the cycle is surmised by non-equilibrium Monte Carlo simulation utilizing free-energy minimization. The trajectory shows myosin transduction of free energy to mechanical work giving evidence for: (i) a causal relationship between product release and work production in the native isoform that is correctly disrupted in a chemically modified protein, (ii) the molecular basis of ATP-sensitive tryptophan fluorescence enhancement and acrylamide quenching, (iii) an actin-binding site peptide containing the free-energy barrier to ATPase product release defining the rate limiting step and, (iv) a scenario for actin-activation of myosin ATPase.

Enhanced sampling via strong coupling to a heat bath: relationship between Tsallis and multicanonical algorithms

We show that Tsallis and multicanonical statistical mechanics are equivalent under specific conditions and that they describe a system strongly coupled to a heat bath. The concept of the strong coupling to a heat bath, in which energy fluctuation is larger than that in the canonical ensemble [J. Chem. Phys. 119, 7075 (2003)], plays a key role in relating Tsallis formalism to multicanonical formalism. The equivalence between these formalisms allows us to obtain an appropriate q parameter in the Tsallis algorithm to enhance the sampling in the phase space in a manner similar to the multicanonical algorithm. An enhanced sampling in the configurational space by use of the strong coupling formalism is demonstrated in a Lennard-Jones fluid.

Generalized simulated annealing applied to protein folding studies

During the last few years, computational simulations based on the atomic description of biological molecules have resulted in significant advances in the comprehension of biological processes. It is well known, however, that a molecular system may have a great number of conformations due to the large number of rotation of degrees of freedom around chemical bonds, leading to several local minima on the energy hypersurface. It has been proposed though, that proteins express their biological function when their structure is close to a conformation with energy global minimum. To help solve the protein-folding problem, we use a new strategy based on Simulated Annealing methods. These methods have been well suited for a large extent of optimization problems, especially those containing many local minima. In fact, this work applies the Generalized Simulated Annealing method (GSA) coupled to the GROMOS96 Molecular Force Field to research the minimum energy conformation of 18-alanine. We show that the q(T) GSA parameter can be used to control the freezing process during the annealing procedure, and to avoid polypeptide chains to be trapped in energy local minima. We scanned the q-values for visiting (q(V)) and accepting (q(A)) functions for q(T) values ranging from 1 to 3, and found the best values to obtain an alpha-helix conformation for the polyalanine peptide, which is the conformation with energy global minimum. Global optimization methods also exemplify a class of applications that requires a large amount of computational resources, being suitable for Grid computing. To implement a Grid computing platform, we developed and tested a Grid environment based on MYGRID middleware, which is a technology that can employ all machines accessed by the user to run the application.

Folding simulations with novel conformational search method

A novel scheme for fast conformational search has been developed by combining the replica exchange method (REM) with the generalized effective potential concept. The new method, referred to Q-REM [S. Jang et al. Phys. Rev. Lett. 91, 058305 (2003)], is expected to provide a useful alternative to the conventional REM for effective conformational sampling of complex systems. The authors have performed folding simulations of the Trp-cage miniprotein using Q-REM. All atom level simulations with generalized Born solvent access-area solvation model show that successful folding can be observed with much smaller number of replicas in Q-REM compared to the conventional REM. It can be concluded that the new method has potential to significantly improve sampling efficiency, allowing simulations of more challenging systems.

Conformational sampling via a self-regulating effective energy surface

The difficulty of efficiently sampling the phase space of complex systems with rough energy surfaces is well known. Typical solutions to the problem involve accelerating the crossing of barriers, but such methods often have the secondary problem that the low-energy states of interest are inadequately sampled, unless the parameters of the search algorithm are modified as the system evolves. A method is presented to improve the sampling with particular emphasis on the low-energy conformations, which make the most important contributions to the thermodynamics of the system. The algorithm proposed here samples the details of the minima, while easily surmounting barriers. This is achieved by introducing a self-regulating sampling variable which depends on the current state of the system. Two replicas of the system are introduced and the sampling variable is treated as a particle coupled to the physical system. The method is illustrated with a simple model system and is applied to the realistic example of barrier crossing in a protein-ligand complex.

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.

Monte Carlo simulations of biomolecules: The MC module in CHARMM

We describe the implementation of a general and flexible Monte Carlo (MC) module for the program CHARMM, which is used widely for modeling biomolecular systems with empirical energy functions. Construction and use of an almost arbitrary move set with only a few commands is made possible by providing several predefined types of moves that can be combined. Sampling can be enhanced by noncanonical acceptance criteria, automatic optimization of step sizes, and energy minimization. A systematic procedure for improving MC move sets is introduced and applied to simulations of two peptides. The resulting move sets allow MC to sample the configuration spaces of these systems much more rapidly than Langevin dynamics. The rate of convergence of the difference in free energy between ethane and methanol in explicit solvent is also examined, and comparable performances are observed for MC and the Nosé-Hoover algorithm. Its ease of use combined with its sampling efficiency make the MC module in CHARMM an attractive alternative for exploring the behavior of biomolecular systems.

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.

On the Convergence of Parallel Tempering Monte Carlo Simulations of LJ38

The convergence of parallel tempering Monte Carlo simulations of the 38-atom Lennard-Jones cluster starting from the Oh global minimum and from the C(5v) second-lowest-energy minimum is investigated. It is found that achieving convergence is appreciably more difficult, particularly at temperatures in the vicinity of the Oh --> C(5v) transformation when starting from the C(5v) structure. A strategy combining the Tsallis generalized ensemble and the parallel tempering algorithm is implemented and used to improve the convergence of the simulations in the vicinity of the Oh --> C(5v) transformation.

Generalized simulated tempering realized on expanded ensembles of non-Boltzmann weights

A generalized version of the simulated tempering operated in the expanded ensembles of non-Boltzmann weights has been proposed to mitigate a quasiergodicity problem occurring in simulations of rough energy landscapes. In contrast to conventional simulated tempering employing the Boltzmann weight, our method utilizes a parametrized, generalized distribution as a workhorse for stochastic exchanges of configurations and subensembles transitions, which allows a considerable enhancement for the rate of convergence of Monte Carlo and molecular dynamics simulations using delocalized weights. A feature of our method is that the exploration of the parameter space encouraging subensembles transitions is greatly accelerated using the dynamic update scheme for the weight via the average guide specific to the energy distribution. The performance and characteristic feature of our method have been validated in the liquid-solid transition of Lennard-Jones clusters and the conformational sampling of alanine dipeptide by taking two types of Tsallis [C. Tsallis, J. Stat. Phys. 52, 479 (1988)] expanded ensembles associated with different parametrization schemes.

Dynamical origin of enhanced conformational searches of Tsallis statistics sampling

The characteristic sampling dynamics of importance samplings driven by the Tsallis weight [C. Tsallis, J. Stat. Phys. 52, 479 (1988)] has been analyzed in terms of recently developed Langevin stochastic model by considering the effects of the density of states and the potential smoothing of the Tsallis transformation. Our study reveals that the fixed points, which are determined by the crossing points of the statistical temperature and the Tsallis effective temperature, play a critical role in overall dynamics of the Tsallis statistics sampling. The dynamical origin of enhanced conformational searches of the Tsallis weight has been investigated by unveiling the intimate relationship between the sampling dynamics and the stability change of corresponding fixed points. Based on this stochastic analysis, we propose one effective method to realize a broad energy distribution in the Tsallis statistics sampling by determining optimal Tsallis parameters systematically based on preliminary canonical samplings. The effectiveness of our method has been validated in the folding simulation of Met-Enkephalin and liquid-solid transition simulation of Lennard-Jones cluster systems.

Stochastic formulation of sampling dynamics in generalized ensemble methods

The fundamental relations of the sampling process in the multicanonical ensemble and simulated tempering have been studied in the expanded ensembles formalism. The simulated tempering is identified as the multicanonical sampling with the generalized weight determined by a Laplace transform of the temperature weight. The characteristic dynamics of both sampling methods has been verified in the stochastic formulation of the sampling process. Our study gives a necessary and sufficient condition for the weights to realize a uniform sampling in the energy and temperature spaces. Based on the stochastic model, robust force biased iteration schemes have been proposed to allow automatic determinations of uniform sampling weights.

Replica-exchange method using the generalized effective potential

We propose an effective scheme for fast conformational searches by combining the replica exchange method (REM) with the generalized effect potential concept. The present method introduces the "q" value from the effective potential as a coupling parameter. It is found that the new method not only requires a much smaller number of replicas than the conventional REM, but also makes it possible to perform effective conformational sampling of complex systems with correct distributions maintained. The advantage of the present method has been demonstrated with in vacuo alanine dipeptide using a molecular dynamics simulation.

Itineration of the Internet over nonequilibrium stationary states in Tsallis statistics

The cumulative probability distribution of sparseness time interval in the Internet is studied by the method of data analysis. Round-trip time between a local host and a destination host through ten odd routers is measured using the ping command, i.e., doing an echo experiment. The data are found to be well described by q-exponential distributions, which maximize the Tsallis entropy indexed by q less or larger than unity, showing a scale-invariant feature of the system. The network is observed to itinerate over a series of the nonequilibrium stationary states characterized by Tsallis statistics.

Global optimization on an evolving energy landscape

Locating the global minimum of a complex potential energy surface is facilitated by considering a homotopy, namely, a family of surfaces that interpolate continuously from an arbitrary initial potential to the system under consideration. Different strategies can be used to follow the evolving minima. It is possible to enhance the probability of locating the global minimum through a heuristic choice of interpolation schemes and parameters, and the continuously evolving potential landscape reduces the probability of trapping in local minima. In application to a model problem, finding the ground-state configuration and the energy of rare-gas (Lennard-Jones) atomic clusters, we demonstrate the utility and the efficacy of this method.

Gravitational smoothing as a global optimization strategy

An optimization scheme for atomic cluster structures, based on exaggerating the importance of the gravitational force, is introduced. Results are presented for calculations on Lennard-Jones clusters of 13, 38, and 55 atoms, and the 13-atom Morse cluster.

Enhanced sampling in numerical path integration: an approximation for the quantum statistical density matrix based on the nonextensive thermostatistics

Here we examine a proposed approximation for the quantum statistical density matrix motivated by the nonextensive thermostatistics of Tsallis and co-workers. The approximation involves replacing the physical potential energy with an effective one, corresponding to a generalized nonextensive statistical ensemble. We examine the convergence properties of averages calculated using the effective potential, and introduce a related method for enhanced sampling in numerical path integration. As a necessary measure, path integral energy estimators are introduced for potentials that involve explicit temperature dependence. This sampling method is found to be effective for path integral simulations involving broken ergodicity.