Generic, hierarchical framework for massively parallel Wang-Landau sampling

Mean Field Initialization of the Annealed Importance Sampling Algorithm for an Efficient Evaluation of the Partition Function Using Restricted Boltzmann Machines

Probabilistic models in physics often require the evaluation of normalized Boltzmann factors, which in turn implies the computation of the partition function Z. Obtaining the exact value of Z, though, becomes a forbiddingly expensive task as the system size increases. A possible way to tackle this problem is to use the Annealed Importance Sampling (AIS) algorithm, which provides a tool to stochastically estimate the partition function of the system. The nature of AIS allows for an efficient and parallel implementation in Restricted Boltzmann Machines (RBMs). In this work, we evaluate the partition function of magnetic spin and spin-like systems mapped into RBMs using AIS. So far, the standard application of the AIS algorithm starts from the uniform probability distribution and uses a large number of Monte Carlo steps to obtain reliable estimations of Z following an annealing process. We show that both the quality of the estimation and the cost of the computation can be significantly improved by using a properly selected mean-field starting probability distribution. We perform a systematic analysis of AIS in both small- and large-sized problems, and compare the results to exact values in problems where these are known. As a result, we propose two successful strategies that work well in all the problems analyzed. We conclude that these are good starting points to estimate the partition function with AIS with a relatively low computational cost. The procedures presented are not linked to any learning process, and therefore do not require a priori knowledge of a training dataset.

A Perspective on the Investigation of Spectroscopy and Kinetics of Complex Molecular Systems with Semiclassical Approaches

In this Perspective we show that semiclassical methods provide a rigorous hierarchical way to study the vibrational spectroscopy and kinetics of complex molecular systems. The time averaged approach to spectroscopy and the semiclassical transition state theory for kinetics, which have been first adopted and then further developed in our group, provide accurate quantum results on rigorous physical grounds and can be applied even when dealing with a large number of degrees of freedom. In spectroscopy, the multiple coherent, divide-and-conquer, and adiabatically switched semiclassical approaches have practically permitted overcoming issues related to the convergence of results. In this Perspective we demonstrate the possibility of studying the semiclassical vibrational spectroscopy of a molecule adsorbed on an anatase (101) surface, a system made of 51 atoms. In kinetics, the semiclassical transition state theory is able to account for anharmonicity and the coupling between the reactive and bound modes. Our group has developed this technique for practical applications involving the study of phenomena like kinetic isotope effect, heavy atom tunneling, and elusive conformer lifetimes. Here, we show that our multidimensional anharmonic quantum approach is able to tackle on-the-fly the thermal kinetic rate constant of a 135 degree-of-freedom system. Overall, semiclassical methods open up the possibility to describe at the quantum mechanical level systems characterized by hundreds of degrees of freedom leading to the accurate spectroscopic and kinetic description of biomolecules and complex molecular systems.

Machine Learning Force Field-Aided Cluster Expansion Approach to Phase Diagram of Alloyed Materials

First-principles approaches based on density functional theory (DFT) have played important roles in the theoretical study of multicomponent alloyed materials. Considering the highly demanding computational cost of direct DFT-based sampling of the configurational space, it is crucial to build efficient and low-cost surrogate Hamiltonian models with DFT accuracy for efficient simulation of alloyed systems with configurational disorder. Recently, the machine learning force field (MLFF) method has been proposed to tackle complicated multicomponent disordered systems. However, the importance of integrating significant physical considerations, including, in particular, convex hull preservation, which is the prerequisite for the accurate prediction of phase diagrams, into the training process of the MLFF remains rarely addressed. In this work, a workflow is proposed to train a convex-hull-preserved (CHP) MLFF for binary alloy systems, based on which the order-disorder phase boundary is predicted by using the Wang-Landau Monte Carlo (WLMC) technique. The predicted values for order-disorder phase transition temperatures agree well with the experiment. The CHP-MLFF is further used to build CE models with the same accuracy as the MLFF and higher efficiency in sampling configurational space. Using the results obtained from the MLFF-based WLMC simulation as a reference, the performances of different schemes for constructing CE models were evaluated in a transparent manner, which revealed the close correlation between the prediction accuracy of ground-state configurations and that of the order-disorder phase transition temperature. This work clearly indicates the great importance of reproducing the convex hull and energetics of ground-state configurations when constructing surrogate Hamiltonians for the statistical modeling of alloyed systems.

Predicting Ion Sequestration in Charged Polymers with the Steepest-Entropy-Ascent Quantum Thermodynamic Framework

The steepest-entropy-ascent quantum thermodynamic framework is used to investigate the effectiveness of multi-chain polyethyleneimine-methylenephosphonic acid in sequestering rare-earth ions (Eu3+) from aqueous solutions. The framework applies a thermodynamic equation of motion to a discrete energy eigenstructure to model the binding kinetics of europium ions to reactive sites of the polymer chains. The energy eigenstructure is generated using a non-Markovian Monte Carlo model that estimates energy level degeneracies. The equation of motion is used to determine the occupation probability of each energy level, describing the unique path through thermodynamic state space by which the polymer system sequesters rare-earth ions from solution. A second Monte Carlo simulation is conducted to relate the kinetic path in state space to physical descriptors associated with the polymer, including the radius of gyration, tortuosity, and Eu-neighbor distribution functions. These descriptors are used to visualize the evolution of the polymer during the sequestration process. The fraction of sequestered Eu3+ ions depends upon the total energy of the system, with lower energy resulting in greater sequestration. The kinetics of the overall sequestration are dependent on the steepest-entropy-ascent principle used by the equation of motion to generate a unique kinetic path from an initial non-equilibrium state.

Phase-space path-integral representation of the quantum density of states: Monte Carlo simulation of strongly correlated soft-sphere fermions

The Wigner formulation of quantum mechanics is used to derive a path-integral representation of the quantum density of states (DOS) of strongly correlated fermions in the canonical ensemble. A path-integral Monte Carlo approach for the simulation of DOS and other thermodynamic functions is suggested. The derived Wigner function in the phase space resembles the Maxwell-Boltzmann distribution but allows for quantum effects. We consider a three-dimensional quantum system of strongly correlated soft-sphere fermions at different densities and temperatures. The calculated properties include the DOS, momentum distribution functions, spin-resolved radial distribution functions, potentials of mean force, and related energy levels obtained from the Bohr-Sommerfeld condition. We observe sharp peaks on DOS and momentum distribution curves, which are explained by the appearance of fermionic bound states.

Predicting Polymer Brush Behavior in Solvents Using the Steepest-Entropy-Ascent Quantum Thermodynamic Framework

The steepest-entropy-ascent quantum thermodynamic (SEAQT) framework is utilized to study the effects of temperature on polymer brushes. The brushes are represented by a discrete energy spectrum, and energy degeneracies obtained through the replica-exchange Wang-Landau algorithm. The SEAQT equation of motion is applied to the density of states to establish a unique kinetic path from an initial thermodynamic state to a stable equilibrium state. The kinetic path describes the brush's evolution in state space, as it interacts with a thermal reservoir. The predicted occupation probabilities along the kinetic path are used to determine the expected thermodynamic and structural properties. The polymer density profile of a polystyrene brush in cyclohexane solvent is predicted using the equation of motion, and it agrees qualitatively with the experimental density profiles. The Flory-Huggins parameter chosen to describe brush-solvent interactions affects the solvent distribution in the brush but has a minimal impact on the polymer density profile. Three types of nonequilibrium kinetic paths with differing amounts of entropy production are considered: a heating path, a cooling path, and a heating-cooling path. Properties such as tortuosity, radius of gyration, brush density, solvent density, and brush chain conformations are calculated for each path.

Predicting non-equilibrium folding behavior of polymer chains using the steepest-entropy-ascent quantum thermodynamic framework

The steepest-entropy-ascent quantum thermodynamic (SEAQT) framework is used to explore the influence of heating and cooling on polymer chain folding kinetics. The framework predicts how a chain moves from an initial non-equilibrium state to stable equilibrium along a unique thermodynamic path. The thermodynamic state is expressed by occupation probabilities corresponding to the levels of a discrete energy landscape. The landscape is generated using the Replica Exchange Wang-Landau method applied to a polymer chain represented by a sequence of hydrophobic and polar monomers with a simple hydrophobic-polar amino acid model. The chain conformation evolves as energy shifts among the levels of the energy landscape according to the principle of steepest entropy ascent. This principle is implemented via the SEAQT equation of motion. The SEAQT framework has the benefit of providing insight into structural properties under non-equilibrium conditions. Chain conformations during heating and cooling change continuously without sharp transitions in morphology. The changes are more drastic along non-equilibrium paths than along quasi-equilibrium paths. The SEAQT-predicted kinetics are fitted to rates associated with the experimental intensity profiles of cytochrome c protein folding with Rouse dynamics.

Predicting defect stability and annealing kinetics in two-dimensional PtSe2using steepest entropy ascent quantum thermodynamics

The steepest-entropy-ascent quantum thermodynamic (SEAQT) framework was used to calculate the stability of a collection of point defects in 2D PtSe₂and predict the kinetics with which defects rearrange during thermal annealing. The framework provides a non-equilibrium, ensemble-based framework with a self-consistent link between mechanics (both quantum and classical) and thermodynamics. It employs an equation of motion derived from the principle of steepest entropy ascent (maximum entropy production) to predict the time evolution of a set of occupation probabilities that define the states of a system undergoing a non-equilibrium process. The system is described by a degenerate energy landscape of eigenvalues, and the entropy is found from the occupation probabilities and the eigenlevel degeneracies. Scanning tunneling microscopy was used to identify the structure and distribution of point defects observed experimentally in a 2D PtSe₂film. A catalog of observed defects includes six unique point defects (vacancies and anti-site defects on Pt and Se sublattices) and twenty combinations of multiple point defects in close proximity. The defect energies were estimated with density functional theory, while the degeneracies, or density of states, for the 2D film with all possible combinations or arrangements of cataloged defects was constructed using a non-Markovian Monte-Carlo approach (i.e. the Replica-Exchange-Wang-Landau algorithm (Vogelet al2013Phys. Rev. Lett.110210603)) with a q-state Potts model. The energy landscape and associated degeneracies were determined for a 2D PtSe₂film two molecules thick and30×30unit cells in area (total of 5400 atoms). The SEAQT equation of motion was applied to the energy landscape to determine how an arbitrary density and arrangement of the six defect types evolve during annealing. Two annealing processes were modeled: heating from 77 K (-196 ∘C) to 523 K (250 ^∘C) and isothermal annealing at 523 K. The SEAQT framework predicted defect configurations, which were consistent with experimental STM images.

Configurational entropy, transition rates, and optimal interactions for rapid folding in coarse-grained model proteins

Under certain conditions, the dynamics of coarse-grained models of solvated proteins can be described using a Markov state model, which tracks the evolution of populations of configurations. The transition rates among states that appear in the Markov model can be determined by computing the relative entropy of states and their mean first passage times. In this paper, we present an adaptive method to evaluate the configurational entropy and the mean first passage times for linear chain models with discontinuous potentials. The approach is based on event-driven dynamical sampling in a massively parallel architecture. Using the fact that the transition rate matrix can be calculated for any choice of interaction energies at any temperature, it is demonstrated how each state's energy can be chosen such that the average time to transition between any two states is minimized. The methods are used to analyze the optimization of the folding process of two protein systems: the crambin protein and a model with frustration and misfolding. It is shown that the folding pathways for both systems are comprised of two regimes: first, the rapid establishment of local bonds, followed by the subsequent formation of more distant contacts. The state energies that lead to the most rapid folding encourage multiple pathways, and they either penalize folding pathways through kinetic traps by raising the energies of trapping states or establish an escape route from the trapping states by lowering free energy barriers to other states that rapidly reach the native state.

Effects of lattice constraints in coarse-grained protein models

We compare and contrast folding behavior in several coarse-grained protein models, both on- and off-lattice, in an attempt to uncover the effect of lattice constraints in these kinds of models. Using modern, extended ensemble Monte Carlo methods-Wang-Landau sampling, multicanonical sampling, replica-exchange Wang-Landau sampling, and replica-exchange multicanonical sampling, we investigate the thermodynamic and structural behavior of the protein Crambin within the context of the hydrophobic-polar, hydrophobic-"neutral"-polar (H0P), and semi-flexible H0P model frameworks. We uncover the folding process in all cases; all models undergo, at least, the two major structural transitions observed in nature-the coil-globule collapse and the folding transition. As the complexity of the model increases, these two major transitions begin to split into multi-step processes, wherein the lattice coarse-graining has a significant impact on the details of these processes. The results show that the level of structural coarse-graining is coupled to the level of interaction coarse-graining.

Calculation of the residual entropy of Ice Ih by Monte Carlo simulation with the combination of the replica-exchange Wang-Landau algorithm and multicanonical replica-exchange method

We estimated the residual entropy of Ice Ih by the recently developed simulation protocol, namely, the combination of the replica-exchange Wang-Landau algorithm and multicanonical replica-exchange method. We employed a model with the nearest neighbor interactions on the three-dimensional hexagonal lattice, which satisfied the ice rules in the ground state. The results showed that our estimate of the residual entropy is in accordance with various previous results. In this article, we not only give our latest estimate of the residual entropy of Ice Ih but also discuss the importance of the uniformity of a random number generator in Monte Carlo simulations.

First-principles calculation of the configurational energy density of states for a solid-state ion conductor with a variant of the Wang and Landau algorithm

In this work, a variant of the Wang and Landau algorithm for calculation of the configurational energy density of states is proposed. The algorithm was developed for the purpose of using first-principles simulations, such as density functional theory, to calculate the partition function of disordered sublattices in crystal materials. The expensive calculations of first-principles methods make a parallel algorithm necessary for a practical computation of the configurational energy density of states within a supercell approximation of a solid-state material. The algorithm developed in this work is tested with the two-dimensional (2d) Ising model to bench mark the algorithm and to help provide insight for implementation to a materials science application. Tests with the 2d Ising model revealed that the algorithm has good performance compared to the original Wang and Landau algorithm and the 1/t algorithm, in particular the short iteration performance. A proof of convergence is presented within an adiabatic assumption, and the analysis is able to correctly predict the time dependence of the modification factor to the density of states. The algorithm was then applied to the lithium and lanthanum sublattice of the solid-state lithium ion conductor Li_{0.5}La_{0.5}TiO_{3}. This was done to help understand the disordered nature of the lithium and lanthanum. The results find, overall, that the algorithm performs very well for the 2d Ising model and that the results for Li_{0.5}La_{0.5}TiO_{3} are consistent with experiment while providing additional insight into the lithium and lanthanum ordering in the material. The primary result is that the lithium and lanthanum become more mixed between layers along the c axis for increasing temperature. In part, the simulation of the disordered Li_{0.5}La_{0.5}TiO_{3} system serves as a benchmark for what size systems are currently and in the near future practical to calculate with density functional theory methods.

Ising universality in the two-dimensional Blume-Capel model with quenched random crystal field

Using high-precision Monte Carlo simulations based on a parallel version of the Wang-Landau algorithm and finite-size scaling techniques, we study the effect of quenched disorder in the crystal-field coupling of the Blume-Capel model on a square lattice. We mainly focus on the part of the phase diagram where the pure model undergoes a continuous transition, known to fall into the universality class of a pure Ising ferromagnet. A dedicated scaling analysis reveals concrete evidence in favor of the strong universality hypothesis with the presence of additional logarithmic corrections in the scaling of the specific heat. Our results are in agreement with an early real-space renormalization-group study of the model as well as a very recent numerical work where quenched randomness was introduced in the energy exchange coupling. Finally, by properly fine tuning the control parameters of the randomness distribution we also qualitatively investigate the part of the phase diagram where the pure model undergoes a first-order phase transition. For this region, preliminary evidence indicate a smoothing of the transition to second-order with the presence of strong scaling corrections.

Flat-histogram method comparison on the two-dimensional Ising model

We compare the convergence of several flat-histogram methods applied to the two-dimensional Ising model, including the recently introduced stochastic approximation with a dynamic update factor (SAD) method. We compare this method to the Wang-Landau (WL) method, the 1/t variant of the WL method, and standard stochastic approximation Monte Carlo (SAMC). In addition, we consider a procedure WL followed by a "production run" with fixed weights that refines the estimation of the entropy. We find that WL followed by a production run does converge to the true density of states, in contrast to pure WL. Three of the methods converge robustly: SAD, 1/t-WL, and WL followed by a production run. Of these, SAD does not require a priori knowledge of the energy range. This work also shows that WL followed by a production run performs superior to other forms of WL while ensuring both ergodicity and detailed balance.

Wang-Landau algorithm as stochastic optimization and its acceleration

We show that the Wang-Landau algorithm can be formulated as a stochastic gradient descent algorithm minimizing a smooth and convex objective function, of which the gradient is estimated using Markov chain Monte Carlo iterations. The optimization formulation provides us another way to establish the convergence rate of the Wang-Landau algorithm, by exploiting the fact that almost surely the density estimates (on the logarithmic scale) remain in a compact set, upon which the objective function is strongly convex. The optimization viewpoint motivates us to improve the efficiency of the Wang-Landau algorithm using popular tools including the momentum method and the adaptive learning rate method. We demonstrate the accelerated Wang-Landau algorithm on a two-dimensional Ising model and a two-dimensional ten-state Potts model.

Logarithmic finite-size scaling correction to the leading Fisher zeros in the p-state clock model: A higher-order tensor renormalization group study

We investigate the finite-size-scaling (FSS) behavior of the leading Fisher zero of the partition function in the complex temperature plane in the p-state clock models of p=5 and 6. We derive the logarithmic finite-size corrections to the scaling of the leading zeros which we numerically verify by performing the higher-order tensor renormalization group (HOTRG) calculations in the square lattices of a size up to 128×128 sites. The necessity of the deterministic HOTRG method in the clock models is noted by the extreme vulnerability of the numerical leading zero identification against stochastic noises that are hard to be avoided in the Monte Carlo approaches. We characterize the system-size dependence of the numerical vulnerability of the zero identification by the type of phase transition, suggesting that the two transitions in the clock models are not of an ordinary first- or second-order type. In the direct FSS analysis of the leading zeros in the clock models, we find that their FSS behaviors show excellent agreement with our predictions of the logarithmic corrections to the Berezinskii-Kosterlitz-Thouless ansatz at both of the high- and low-temperature transitions.

Efficient simulation protocol for determining the density of states: Combination of replica-exchange Wang-Landau method and multicanonical replica-exchange method

By combining two generalized-ensemble algorithms, the replica-exchange Wang-Landau (REWL) method and the multicanonical replica-exchange method (MUCAREM), we propose an effective simulation protocol to determine the density of states with high accuracy. The new protocol is referred to as REWL-MUCAREM, and REWL is first performed and then MUCAREM is performed. In order to verify the effectiveness of our protocol, we performed simulations of a square-lattice Ising model using the three methods, namely REWL, MUCAREM, and REWL-MUCAREM. The results showed that the density of states obtained by REWL-MUCAREM is more accurate than that is estimated by the two methods separately.

Morphology and thermodynamics of polymers with monofunctional hydrogen bonding ends in dilute and semidilute concentration

Rheological properties of supramolecular polymers (SMPs) depend on their equilibrium structure including the size, the number, and the topology of aggregates. A polymer with a hydrogen bonding (H-bonding) motif at both ends is one widely used precursor to build SMPs. Due to the complex interplay between chain stiffness, H-bonding interaction, polarity along a chain, and polymer conformational entropy, it is difficult to theoretically predict the structure of SMPs. In this work we investigate thermodynamics of SMPs with H-bonding ends in a wide range of densities. A replica exchange stochastic approximation Monte Carlo method with coarse-grained models for polyethylene and polybuthylene glycols is used. Our simulation shows that SMPs have two morphological transition lines with increasing temperature, a ring-linear transition, and a linear-free chain transition. The latter is a thermodynamic transition and turns out to be continuous. Comparing the two different spacers, we find that ring-linear transition temperatures differ from each other at the constant volume fraction due to different looping probabilities, which can be calculated from the average polymer size by mean field. However, the linear-free chain transition temperatures are similar because the entropic penalty to form a hydrogen bond mainly depends on the probability of finding H-bonding groups in a system, which is the same for both systems at a given volume fraction.

First-principles study of order-disorder transitions in multicomponent solid-solution alloys

In this review, we will focus on the recent development of the order-disorder transition in metallic materials. The past decades have witnessed fast development in the first-principles methodologies and their applications to ordering transitions in multi-component alloys, particularly the high-entropy alloys. The driving force for the proceedings comes from (i) the advance of algorithms and increasingly cheaper hardware, and also (ii) the great passion to model alloys with increasing number of components. The review starts with a brief introduction of the history for the ordering transitions. More detailed scientific proceedings prior to the 1970s had been well summarized in Krivoglaz and Smirnov (1965 The Theory of Order-Disorder in Alloys (New York: Elsevier)) and Stoloff and Davies (1968 Prog. Mater. Sci. 13 1-84). In the second part, the methods to study the ordering transitions, primarily on the theoretic methods are introduced. These will include (i) KKR-CPA method and supercell methods for energetic calculations; and (ii) thermodynamic and statistical methods to compute the transition temperatures. The third part will focus on representative applications in alloys, including our own work and many others. This part supplies the primary information of this review to the readers. The fourth part will summarize the connections between ordering transitions and broader physical properties (e.g. the mechanical properties). In the last part, some concluding remarks and perspectives will be given.

Parallel Implementation of Semiclassical Transition State Theory

This paper presents the parsctst code, an efficient parallel implementation of the semiclassical transition state theory (SCTST) for reaction rate constant calculations. Parsctst is developed starting from a previously presented approach for the computation of the vibrational density of states of fully coupled anharmonic molecules ( Nguyen et al. Chem. Phys. Lett. 2010 , 499 , 915 ). The parallel implementation makes it practical to tackle reactions involving more than 100 fully coupled anharmonic vibrational degrees of freedom and also includes multidimensional tunneling effects. After describing the pseudocode and demonstrating its computational efficiency, we apply the new code for estimating the rate constant of the proton transfer isomerization reaction of the 2,4,6-tri- tert-butylphenyl to 3,5-di- tert-butylneophyl. Comparison with both theoretical and experimental results is presented. Parsctst code is user-friendly and provides a significant computational time saving compared to serial calculations. We believe that parsctst can boost the application of SCTST as an alternative to the basic transition state theory for accurate kinetics modeling not only in combustion or atmospheric chemistry, but also in organic synthesis, where bigger reactive systems are encountered.

Direct coupling of first-principles calculations with replica exchange Monte Carlo sampling of ion disorder in solids

We demonstrate the feasibility of performing sufficient configurational sampling of disordered oxides directly from first-principles without resorting to the use of fitted models such as cluster expansion. This is achieved by harnessing the power of modern-day cluster supercomputers using the replica exchange Monte Carlo method coupled directly with structural relaxation and energy calculation performed by density functional codes. The idea is applied successfully to the calculation of the temperature-dependence of the degree of inversion in the cation sublattice of MgAl₂O₄ spinel oxide. The possibility of bypassing fitting models will lead to investigation of disordered systems where cluster expansion is known to perform badly, for example, systems with large lattice deformation due to defects, or systems where long-range interactions dominate such as electrochemical interfaces.

Elucidating thermal behavior, native contacts, and folding funnels of simple lattice proteins using replica exchange Wang-Landau sampling

We studied the folding behavior of two coarse-grained, lattice models, the HP (hydrophobic-polar) model and the semi-flexible H0P model, whose 124 monomer long sequences were derived from the protein Ribonuclease A. Taking advantage of advanced parallel computing techniques, we applied replica exchange Wang-Landau sampling and calculated the density of states over the models entire energy ranges to high accuracy. We then determined both energetic and structural quantities in order to elucidate the folding behavior of each model completely. As a result of sufficiently long sequences and model complexity, yet computational accessibility, we were able to depict distinct free energy folding funnels for both models. In particular, we found that the HP model folds in a single-step process with a very highly degenerate native state and relatively flat low temperature folding funnel minimum. By contrast, the semi-flexible H0P model folds via a multi-step process and the native state is almost four orders of magnitude less degenerate than that for the HP model. In addition, for the H0P model, the bottom of the free energy folding funnel remains rough, even at low temperatures.

The role of chain-stiffness in lattice protein models: A replica-exchange Wang-Landau study

Using Monte Carlo simulations, we investigate simple, physically motivated extensions to the hydrophobic-polar lattice protein model for the small (46 amino acid) protein Crambin. We use two-dimensional replica-exchange Wang-Landau sampling to study the effects of a bond angle stiffness parameter on the folding and uncover a new step in the collapse process for particular values of this stiffness parameter. A physical interpretation of the folding is developed by analysis of changes in structural quantities, and the free energy landscape is explored. For these special values of stiffness, we find non-degenerate ground states, a property that is consistent with behavior of real proteins, and we use these unique ground states to elucidate the formation of native contacts during the folding process. Through this analysis, we conclude that chain-stiffness is particularly influential in the low energy, low temperature regime of the folding process once the lattice protein has partially collapsed.

Complex free-energy landscapes in biaxial nematic liquid crystals and the role of repulsive interactions: A Wang-Landau study

General quadratic Hamiltonian models, describing the interaction between liquid-crystal molecules (typically with D_{2h} symmetry), take into account couplings between their uniaxial and biaxial tensors. While the attractive contributions arising from interactions between similar tensors of the participating molecules provide for eventual condensation of the respective orders at suitably low temperatures, the role of cross coupling between unlike tensors is not fully appreciated. Our recent study with an advanced Monte Carlo technique (entropic sampling) showed clearly the increasing relevance of this cross term in determining the phase diagram (contravening in some regions of model parameter space), the predictions of mean-field theory, and standard Monte Carlo simulation results. In this context, we investigated the phase diagrams and the nature of the phases therein on two trajectories in the parameter space: one is a line in the interior region of biaxial stability believed to be representative of the real systems, and the second is the extensively investigated parabolic path resulting from the London dispersion approximation. In both cases, we find the destabilizing effect of increased cross-coupling interactions, which invariably result in the formation of local biaxial organizations inhomogeneously distributed. This manifests as a small, but unmistakable, contribution of biaxial order in the uniaxial phase. The free-energy profiles computed in the present study as a function of the two dominant order parameters indicate complex landscapes. On the one hand, these profiles account for the unusual thermal behavior of the biaxial order parameter under significant destabilizing influence from the cross terms. On the other, they also allude to the possibility that in real systems, these complexities might indeed be inhibiting the formation of a low-temperature biaxial order itself-perhaps reflecting the difficulties in their ready realization in the laboratory.

Partition function zeros of the p-state clock model in the complex temperature plane

We investigate the partition function zeros of the two-dimensional p-state clock model in the complex temperature plane by using the Wang-Landau method. For p=5, 6, 8, and 10, we propose a modified energy representation to enumerate exact irregular energy levels for the density of states without any binning artifacts. Comparing the leading zeros between different p's, we provide strong evidence that the upper transition at p=6 is indeed of the Berezinskii-Kosterlitz-Thouless (BKT) type in contrast to the claim of the previous Fisher zero study [Phys. Rev. E 80, 042103 (2009)10.1103/PhysRevE.80.042103]. We find that the leading zeros of p=6 at the upper transition collapse onto the zero trajectories of the larger p's including the XY limit while the finite-size behavior of p=5 differs from the converged behavior of p≥6 within the system sizes examined. In addition, we argue that the nondivergent specific heat in the BKT transition is responsible for the small partition function magnitude that decreases exponentially with increasing system size near the leading zero, fundamentally limiting access to large systems in search for zeros with an estimator under finite statistical fluctuations.

Solubility prediction from first principles: a density of states approach

Solubility is a fundamental property of widespread significance. Its accurate prediction remains a major challenge. We present a novel, efficient approach to solubility prediction for molecules over a range of conditions based on density of states.

Convex hulls of random walks in higher dimensions: A large-deviation study

The distribution of the hypervolume V and surface ∂V of convex hulls of (multiple) random walks in higher dimensions are determined numerically, especially containing probabilities far smaller than P=10^{-1000} to estimate large deviation properties. For arbitrary dimensions and large walk lengths T, we suggest a scaling behavior of the distribution with the length of the walk T similar to the two-dimensional case and behavior of the distributions in the tails. We underpin both with numerical data in d=3 and d=4 dimensions. Further, we confirm the analytically known means of those distributions and calculate their variances for large T.

Macroscopically constrained Wang-Landau method for systems with multiple order parameters and its application to drawing complex phase diagrams

A generalized approach to Wang-Landau simulations, macroscopically constrained Wang-Landau, is proposed to simulate the density of states of a system with multiple macroscopic order parameters. The method breaks a multidimensional random-walk process in phase space into many separate, one-dimensional random-walk processes in well-defined subspaces. Each of these random walks is constrained to a different set of values of the macroscopic order parameters. When the multivariable density of states is obtained for one set of values of fieldlike model parameters, the density of states for any other values of these parameters can be obtained by a simple transformation of the total system energy. All thermodynamic quantities of the system can then be rapidly calculated at any point in the phase diagram. We demonstrate how to use the multivariable density of states to draw the phase diagram, as well as order-parameter probability distributions at specific phase points, for a model spin-crossover material: an antiferromagnetic Ising model with ferromagnetic long-range interactions. The fieldlike parameters in this model are an effective magnetic field and the strength of the long-range interaction.

Implications of Occupational Disorder on Ion Mobility in Li4Ti5O12 Battery Materials

Lithium-titanium-oxide (Li₄Ti₅O₁₂, LTO) is unique among battery materials due to its exceptional cyclability and high rate capability. This performance is believed to derive at least partly from the occupational disorder introduced via mixed Li/Ti occupancy in the LTO spinel-like structure. We explore the vast configuration space accessible during high-temperature LTO synthesis by Monte Carlo sampling and indeed find lowest-energy structures to be characterized by a high degree of microscopic inhomogeneity. Dynamical simulations in corresponding configurations reveal the dominant fraction of Li ions to be immobile on nanosecond time scales. However, Ti antisite-like defects stabilized by the configurational disorder give rise to a novel correlated ion diffusion mechanism. The resulting fast but localized diffusion could be a key element in the sudden rise in conductivity found in LTO in the early stages of charging and questions the validity of ion mobility measurements for this and other configurationally disordered materials.

A single-walker approach for studying quasi-nonergodic systems

The jump-walking Monte-Carlo algorithm is revisited and updated to study the equilibrium properties of systems exhibiting quasi-nonergodicity. It is designed for a single processing thread as opposed to currently predominant algorithms for large parallel processing systems. The updated algorithm is tested on the Ising model and applied to the lattice-gas model for sorption in aerogel at low temperatures, when dynamics of the system is critically slowed down. It is demonstrated that the updated jump-walking simulations are able to produce equilibrium isotherms which are typically hidden by the hysteresis effect characteristic of the standard single-flip simulations.

Nanothermodynamics of iron clusters: Small clusters, icosahedral and fcc-cuboctahedral structures

The study of the thermodynamics and structures of iron clusters has been carried on, focusing on small clusters and initial icosahedral and fcc-cuboctahedral structures. Two combined tools are used. First, energy intervals are explored by the Monte Carlo algorithm, called σ-mapping, detailed in the work of Soudan et al. [J. Chem. Phys. 135, 144109 (2011), Paper I]. In its flat histogram version, it provides the classical density of states, g_p(E_p), in terms of the potential energy of the system. Second, the iron system is described by a potential which is called "corrected EAM" (cEAM), explained in the work of Basire et al. [J. Chem. Phys. 141, 104304 (2014), Paper II]. Small clusters from 3 to 12 atoms in their ground state have been compared first with published Density Functional Theory (DFT) calculations, giving a complete agreement of geometries. The series of 13, 55, 147, and 309 atom icosahedrons is shown to be the most stable form for the cEAM potential. However, the 147 atom cluster has a special behaviour, since decreasing the energy from the liquid zone leads to the irreversible trapping of the cluster in a reproducible amorphous state, 7.38 eV higher in energy than the icosahedron. This behaviour is not observed at the higher size of 309 atoms. The heat capacity of the 55, 147, and 309 atom clusters revealed a pronounced peak in the solid zone, related to a solid-solid transition, prior to the melting peak. The corresponding series of 13, 55, and 147 atom cuboctahedrons has been compared, underscoring the unstability towards the icosahedral structure. This unstability occurs clearly in several steps for the 147 atom cluster, with a sudden transformation at a transition state. This illustrates the concerted icosahedron-cuboctahedron transformation of Buckminster Fuller-Mackay, which is calculated for the cEAM potential. Two other clusters of initial fcc structures with 24 and 38 atoms have been studied, as well as a 302 atom cluster. Each one relaxes towards a more stable structure without regularity. The 38 atom cluster exhibits a nearly glassy relaxation, through a cascade of six metastable states of long life. This behaviour, as that of the 147 atom cluster towards the amorphous state, shows that difficulties to reach ergodicity in the lower half of the solid zone are related to particular features of the potential energy landscape, and not necessarily to a too large size of the system. Comparisons of the cEAM iron system with published results about Lennard-Jones systems and DFT calculations are made. The results of the previous clusters have been combined with that of Paper II to plot the cohesive energy E_c and the melting temperature T_m in terms of the cluster atom number N_at. The N_at^-1/3 linear dependence of the melting temperature (Pawlow law) is observed again for N_at > 150. In contrast, for N_at < 150, the curve diverges strongly from the Pawlow law, giving it an overall V-shape, with a linear increase of T_m when N_at goes from 55 to 13 atoms. Surprisingly, the 38 atom cluster is anomalously below the overall curve.

Entropy of diluted antiferromagnetic Ising models on frustrated lattices using the Wang-Landau method

We use a Monte Carlo simulation to study the diluted antiferromagnetic Ising model on frustrated lattices including the pyrochlore lattice to show the dilution effects. Using the Wang-Landau algorithm, which directly calculates the energy density of states, we accurately calculate the entropy of the system. We discuss the nonmonotonic dilution concentration dependence of residual entropy for the antiferromagnetic Ising model on the pyrochlore lattice, and compare it to the generalized Pauling approximation proposed by Ke et al. [Phys. Rev. Lett. 99, 137203 (2007)PRLTAO0031-900710.1103/PhysRevLett.99.137203]. We also investigate other frustrated systems, the antiferromagnetic Ising model on the triangular lattice and the kagome lattice, demonstrating the difference in the dilution effects between the system on the pyrochlore lattice and that on other frustrated lattices.

Subphase transitions in first-order aggregation processes

In this paper, we investigate the properties of aggregation transitions in the context of generic coarse-grained homopolymer systems. By means of parallel replica-exchange Monte Carlo methods, we perform extensive simulations of systems consisting of up to 20 individual oligomer chains with five monomers each. Using the tools of the versatile microcanonical inflection-point analysis, we show that the aggregation transition is a first-order process consisting of a sequence of subtransitions between intermediate structural phases. We unravel the properties of these intermediate phases by collecting and analyzing their individual contributions towards the density of states of the system. The central theme of this systematic study revolves around translational entropy and its role in the striking phenomena of missing intermediate phases. We conclude with a brief discussion of the scaling properties of the transition temperature and the latent heat.

Characterizing folding funnels with replica exchange Wang-Landau simulation of lattice proteins

We have studied the folding of ribonuclease A by mapping it onto coarse-grained lattice protein models. With replica exchange Wang-Landau sampling, we calculated the free energy vs end-to-end distance as a function of temperature. A mapping to the famous hydrophobic-polar (HP) model shows a relatively shallow folding funnel and flat free energy minimum, reflecting the high degeneracy of the ground state. In contrast, extending the HP model with an additional "neutral" monomer type (i.e., a mapping to the three-letter H0P model) has a well developed, rough free energy funnel with a low degeneracy ground state. In both cases, folding funnels are asymmetric with temperature dependent shape.

Magnetic phase transition in coupled spin-lattice systems: A replica-exchange Wang-Landau study

Coupled, dynamical spin-lattice models provide a unique test ground for simulations investigating the finite-temperature magnetic properties of materials under the direct influence of the lattice vibrations. These models are constructed by combining a coordinate-dependent interatomic potential with a Heisenberg-like spin Hamiltonian, facilitating the treatment of both the atomic coordinates and the spins as explicit phase variables. Using a model parameterized for bcc iron, we study the magnetic phase transition in these complex systems via the recently introduced, massively parallel replica-exchange Wang-Landau Monte Carlo method. Comparison with the results obtained from rigid lattice (spin-only) simulations shows that the transition temperature as well as the amplitude of the peak in the specific heat curve is marginally affected by the lattice vibrations. Moreover, the results were found to be sensitive to the particular choice of interatomic potential.

From network reliability to the Ising model: A parallel scheme for estimating the joint density of states

Network reliability is the probability that a dynamical system composed of discrete elements interacting on a network will be found in a configuration that satisfies a particular property. We introduce a reliability property, Ising feasibility, for which the network reliability is the Ising model's partition function. As shown by Moore and Shannon, the network reliability can be separated into two factors: structural, solely determined by the network topology, and dynamical, determined by the underlying dynamics. In this case, the structural factor is known as the joint density of states. Using methods developed to approximate the structural factor for other reliability properties, we simulate the joint density of states, yielding an approximation for the partition function. Based on a detailed examination of why naïve Monte Carlo sampling gives a poor approximation, we introduce a parallel scheme for estimating the joint density of states using a Markov-chain Monte Carlo method with a spin-exchange random walk. This parallel scheme makes simulating the Ising model in the presence of an external field practical on small computer clusters for networks with arbitrary topology with ∼10^{6} energy levels and more than 10^{308} microstates.

Dilute Semiflexible Polymers with Attraction: Collapse, Folding and Aggregation

We review the current state on the thermodynamic behavior and structural phases of self- and mutually-attractive dilute semiflexible polymers that undergo temperature-driven transitions. In extreme dilution, polymers may be considered isolated, and this single polymer undergoes a collapse or folding transition depending on the internal structure. This may go as far as to stable knot phases. Adding polymers results in aggregation, where structural motifs again depend on the internal structure. We discuss in detail the effect of semiflexibility on the collapse and aggregation transition and provide perspectives for interesting future investigations.

Multidimensional stochastic approximation Monte Carlo

Stochastic Approximation Monte Carlo (SAMC) has been established as a mathematically founded powerful flat-histogram Monte Carlo method, used to determine the density of states, g(E), of a model system. We show here how it can be generalized for the determination of multidimensional probability distributions (or equivalently densities of states) of macroscopic or mesoscopic variables defined on the space of microstates of a statistical mechanical system. This establishes this method as a systematic way for coarse graining a model system, or, in other words, for performing a renormalization group step on a model. We discuss the formulation of the Kadanoff block spin transformation and the coarse-graining procedure for polymer models in this language. We also apply it to a standard case in the literature of two-dimensional densities of states, where two competing energetic effects are present g(E_{1},E_{2}). We show when and why care has to be exercised when obtaining the microcanonical density of states g(E_{1}+E_{2}) from g(E_{1},E_{2}).

Nonconvergence of the Wang-Landau algorithms with multiple random walkers

This paper discusses some convergence properties in the entropic sampling Monte Carlo methods with multiple random walkers, particularly in the Wang-Landau (WL) and 1/t algorithms. The classical algorithms are modified by the use of m-independent random walkers in the energy landscape to calculate the density of states (DOS). The Ising model is used to show the convergence properties in the calculation of the DOS, as well as the critical temperature, while the calculation of the number π by multiple dimensional integration is used in the continuum approximation. In each case, the error is obtained separately for each walker at a fixed time, t; then, the average over m walkers is performed. It is observed that the error goes as 1/sqrt[m]. However, if the number of walkers increases above a certain critical value m>m_{x}, the error reaches a constant value (i.e., it saturates). This occurs for both algorithms; however, it is shown that for a given system, the 1/t algorithm is more efficient and accurate than the similar version of the WL algorithm. It follows that it makes no sense to increase the number of walkers above a critical value m_{x}, since it does not reduce the error in the calculation. Therefore, the number of walkers does not guarantee convergence.

Simulating Replica Exchange: Markov State Models, Proposal Schemes, and the Infinite Swapping Limit

Replica exchange molecular dynamics is a multicanonical simulation technique commonly used to enhance the sampling of solvated biomolecules on rugged free energy landscapes. While replica exchange is relatively easy to implement, there are many unanswered questions about how to use this technique most efficiently, especially because it is frequently the case in practice that replica exchange simulations are not fully converged. A replica exchange cycle consists of a series of molecular dynamics steps of a set of replicas moving under different Hamiltonians or at different thermodynamic states followed by one or more replica exchange attempts to swap replicas among the different states. How the replica exchange cycle is constructed affects how rapidly the system equilibrates. We have constructed a Markov state model of replica exchange (MSMRE) using long molecular dynamics simulations of a host-guest binding system as an example, in order to study how different implementations of the replica exchange cycle can affect the sampling efficiency. We analyze how the number of replica exchange attempts per cycle, the number of MD steps per cycle, and the interaction between the two parameters affects the largest implied time scale of the MSMRE simulation. The infinite swapping limit is an important concept in replica exchange. We show how to estimate the infinite swapping limit from the diagonal elements of the exchange transition matrix constructed from MSMRE "simulations of simulations" as well as from relatively short runs of the actual replica exchange simulations.

Thermodynamics and structure of macromolecules from flat-histogram Monte Carlo simulations

Over the last decade flat-histogram Monte Carlo simulations, especially multi-canonical and Wang-Landau simulations, have emerged as a strong tool to study the statistical mechanics of polymer chains. These investigations have focused on coarse-grained models of polymers on the lattice and in the continuum. Phase diagrams of chains in bulk as well as chains attached to surfaces were studied, for homopolymers as well as for protein-like models. Also, aggregation behavior in solution of these models has been investigated. We will present here the theoretical background for these simulations, explain the algorithms used and discuss their performance and give an overview over the systems studied with these methods in the literature, where we will limit ourselves to studies of coarse-grained model systems. Implementations of these algorithms on parallel computers will be also briefly described. In parallel to the development of these simulation methods, the power of a micro-canonical analysis of such simulations has been recognized, and we present the current state of the art in applying the micro-canonical analysis to phase transitions in nanoscopic polymer systems.

Towards an Optimal Flow: Density-of-States-Informed Replica-Exchange Simulations

Replica exchange (RE) is one of the most popular enhanced-sampling simulations technique in use today. Despite widespread successes, RE simulations can sometimes fail to converge in practical amounts of time, e.g., when sampling around phase transitions, or when a few hard-to-find configurations dominate the statistical averages. We introduce a generalized RE scheme, density-of-states-informed RE, that addresses some of these challenges. The key feature of our approach is to inform the simulation with readily available, but commonly unused, information on the density of states of the system as the RE simulation proceeds. This enables two improvements, namely, the introduction of resampling moves that actively move the system towards equilibrium and the continual adaptation of the optimal temperature set. As a consequence of these two innovations, we show that the configuration flow in temperature space is optimized and that the overall convergence of RE simulations can be dramatically accelerated.

Confinement effects upon the separation of structural transitions in linear systems with restricted bond fluctuation ranges

By means of advanced parallel replica-exchange Monte Carlo methods we examine the influence of elasticity and confinement on the structural transitions of linear systems with restricted bonded interaction. For this purpose, we adopt a model for coarse-grained flexible polymers of finite length in the dilute regime. Hyperphase diagrams are constructed using energy-dependent canonical quantities to demonstrate the effects of the changes in the range of the confined interaction on the liquid and solid structural phases. With increasing bonded interaction range we observe the disappearance of the liquid phase and the fusion of the gas-liquid (or Θ) and the liquid-solid transitions. One of the most remarkable features, the liquid-gas transition, changes from second to first order if the confined interaction range exceeds a threshold that separates polymeric from nonpolymeric systems. The notoriously difficult sampling of the entropically suppressed conformations in the region of very strong first-order transitions is improved by using multiple Gaussian modified ensembles.

Reexamination of the mean-field phase diagram of biaxial nematic liquid crystals: Insights from Monte Carlo studies

Investigations of the phase diagram of biaxial liquid-crystal systems through analyses of general Hamiltonian models within the simplifications of mean-field theory (MFT), as well as by computer simulations based on microscopic models, are directed toward an appreciation of the role of the underlying molecular-level interactions to facilitate its spontaneous condensation into a nematic phase with biaxial symmetry. Continuing experimental challenges in realizing such a system unambiguously, despite encouraging predictions from MFT, for example, are requiring more versatile simulational methodologies capable of providing insights into possible hindering barriers within the system, typically gleaned through its free-energy dependences on relevant observables as the system is driven through the transitions. The recent paper from this group [Kamala Latha et al., Phys. Rev. E 89, 050501(R) (2014)], summarizing the outcome of detailed Monte Carlo simulations carried out employing an entropic sampling technique, suggested a qualitative modification of the MFT phase diagram as the Hamiltonian is asymptotically driven toward the so-called partly repulsive regions. It was argued that the degree of (cross) coupling between the uniaxial and biaxial tensor components of neighboring molecules plays a crucial role in facilitating a ready condensation of the biaxial phase, suggesting that this could be a plausible factor in explaining the experimental difficulties. In this paper, we elaborate this point further, providing additional evidence from curious variations of free-energy profiles with respect to the relevant orientational order parameters, at different temperatures bracketing the phase transitions.

First-order phase transition and tricritical scaling behavior of the Blume-Capel model: A Wang-Landau sampling approach

We investigate the tricritical scaling behavior of the two-dimensional spin-1 Blume-Capel model by using the Wang-Landau method of measuring the joint density of states for lattice sizes up to 48×48 sites. We find that the specific heat deep in the first-order area of the phase diagram exhibits a double-peak structure of the Schottky-like anomaly appearing with the transition peak. The first-order transition curve is systematically determined by employing the method of field mixing in conjunction with finite-size scaling, showing a significant deviation from the previous data points. At the tricritical point, we characterize the tricritical exponents through finite-size-scaling analysis including the phenomenological finite-size scaling with thermodynamic variables. Our estimation of the tricritical eigenvalue exponents, yt=1.804(5), yg=0.80(1), and yh=1.925(3), provides the first Wang-Landau verification of the conjectured exact values, demonstrating the effectiveness of the density-of-states-based approach in finite-size scaling study of multicritical phenomena.

Dynamical traps in Wang-Landau sampling of continuous systems: Mechanism and solution

We study the mechanism behind dynamical trappings experienced during Wang-Landau sampling of continuous systems reported by several authors. Trapping is caused by the random walker coming close to a local energy extremum, although the mechanism is different from that of the critical slowing-down encountered in conventional molecular dynamics or Monte Carlo simulations. When trapped, the random walker misses the entire or even several stages of Wang-Landau modification factor reduction, leading to inadequate sampling of the configuration space and a rough density of states, even though the modification factor has been reduced to very small values. Trapping is dependent on specific systems, the choice of energy bins, and the Monte Carlo step size, making it highly unpredictable. A general, simple, and effective solution is proposed where the configurations of multiple parallel Wang-Landau trajectories are interswapped to prevent trapping. We also explain why swapping frees the random walker from such traps. The efficacy of the proposed algorithm is demonstrated.