Flat-Histogram Dynamics and Optimization in Density of States Simulations of Fluids

Direct free energy evaluation of classical and quantum many-body systems via field-theoretic simulation

The accurate evaluation of free energies within molecular simulations is important to many scientific fields, including fluid and solid phase equilibria, biomolecular condensates, and quantum phase transitions, among others. Unfortunately, free energy estimation is tedious and computationally expensive for molecular models whose degrees of freedom are expressed in particle coordinates. We show that alternative representations of a model as a classical or quantum field theory provide access to a chemical potential operator that can be averaged to yield a direct and low-cost estimate of the Gibbs free energy. The averaging is performed using a “field-theoretic” computer simulation that employs fluctuating fields rather than particles.

Free energy evaluation in molecular simulations of both classical and quantum systems is computationally intensive and requires sophisticated algorithms. This is because free energy depends on the volume of accessible phase space, a quantity that is inextricably linked to the integration measure in a coordinate representation of a many-body problem. In contrast, the same problem expressed as a field theory (auxiliary field or coherent states) isolates the particle number as a simple parameter in the Hamiltonian or action functional and enables the identification of a chemical potential field operator. We show that this feature leads a “direct” method of free energy evaluation, in which a particle model is converted to a field theory and appropriate field operators are averaged using a field-theoretic simulation conducted with complex Langevin sampling. These averages provide an immediate estimate of the Helmholtz free energy in the canonical ensemble and the entropy in the microcanonical ensemble. The method is illustrated for a classical polymer solution, a block copolymer melt exhibiting liquid crystalline and solid mesophases, and a quantum fluid of interacting bosons.

Flat-Histogram Monte Carlo as an Efficient Tool To Evaluate Adsorption Processes Involving Rigid and Deformable Molecules

Monte Carlo simulations are the foundational technique for predicting thermodynamic properties of open systems where the process of interest involves the exchange of particles. Thus, they have been used extensively to computationally evaluate the adsorption properties of nanoporous materials and are critical for the in silico identification of promising materials for a variety of gas storage and chemical separation applications. In this work we demonstrate that a well-known biasing technique, known as "flat-histogram" sampling, can be combined with temperature extrapolation of the free energy landscape to efficiently provide significantly more useful thermodynamic information than standard open ensemble MC simulations. Namely, we can accurately compute the isosteric heat of adsorption and number of particles adsorbed for various adsorbates over an extremely wide range of temperatures and pressures from a set of simulations at just one temperature. We extend this derivation of the temperature extrapolation to adsorbates with intramolecular degrees of freedom when Rosenbluth sampling is employed. Consequently, the working capacity and isosteric heat can be computed for any given combined temperature/pressure swing adsorption process for a large range of operating conditions with both rigid and deformable adsorbates. Continuous thermodynamic properties can be computed with this technique at very moderate computational cost, thereby providing a strong case for its application to the in silico identification of promising nanoporous adsorbents.

Relationship between Pore-size Distribution and Flexibility of Adsorbent Materials: Statistical Mechanics and Future Material Characterization Techniques

Measurement of the pore-size distribution (PSD) via gas adsorption and the so-called "kernel method" is a widely used characterization technique for rigid adsorbents. Yet, standard techniques and analytical equipment are not appropriate to characterize the emerging class of flexible adsorbents that deform in response to the stress imparted by an adsorbate gas, as the PSD is a characteristic of the material that varies with the gas pressure and any other external stresses. Here, we derive the PSD for a flexible adsorbent using statistical mechanics in the osmotic ensemble to draw analogy to the kernel method for rigid materials. The resultant PSD is a function of the ensemble constraints including all imposed stresses and, most importantly, the deformation free energy of the adsorbent material. Consequently, a pressure-dependent PSD is a descriptor of the deformation characteristics of an adsorbent and may be the basis of future material characterization techniques. We discuss how, given a technique for resolving pressure-dependent PSDs, the present statistical mechanical theory could enable a new generation of analytical tools that measure and characterize certain intrinsic material properties of flexible adsorbents via otherwise simple adsorption experiments.

Evaluation of the grand-canonical partition function using expanded Wang-Landau simulations. IV. Performance of many-body force fields and tight-binding schemes for the fluid phases of silicon

We extend Expanded Wang-Landau (EWL) simulations beyond classical systems and develop the EWL method for systems modeled with a tight-binding Hamiltonian. We then apply the method to determine the partition function and thus all thermodynamic properties, including the Gibbs free energy and entropy, of the fluid phases of Si. We compare the results from quantum many-body (QMB) tight binding models, which explicitly calculate the overlap between the atomic orbitals of neighboring atoms, to those obtained with classical many-body (CMB) force fields, which allow to recover the tetrahedral organization in condensed phases of Si through, e.g., a repulsive 3-body term that favors the ideal tetrahedral angle. Along the vapor-liquid coexistence, between 3000 K and 6000 K, the densities for the two coexisting phases are found to vary significantly (by 5 orders of magnitude for the vapor and by up to 25% for the liquid) and to provide a stringent test of the models. Transitions from vapor to liquid are predicted to occur for chemical potentials that are 10%-15% higher for CMB models than for QMB models, and a ranking of the force fields is provided by comparing the predictions for the vapor pressure to the experimental data. QMB models also reveal the formation of a gap in the electronic density of states of the coexisting liquid at high temperatures. Subjecting Si to a nanoscopic confinement has a dramatic effect on the phase diagram with, e.g. at 6000 K, a decrease in liquid densities by about 50% for both CMB and QMB models and an increase in vapor densities between 90% (CMB) and 170% (QMB). The results presented here provide a full picture of the impact of the strategy (CMB or QMB) chosen to model many-body effects on the thermodynamic properties of the fluid phases of Si.

Role of hydrophobic hydration in protein stability: a 3D water-explicit protein model exhibiting cold and heat denaturation

We investigate the microscopic mechanism of cold and heat denaturation using a 3D lattice model of a hydrated protein in which water is represented explicitly. The water model, which incorporates directional bonding and tetrahedral geometry, captures many aspects of water thermodynamics and properly describes hydrophobic hydration around apolar solutes because the hydrogen bonding rules in the model were gleaned from off-lattice atomistic simulations of water around representative protein structures. By incorporating local chain stiffness in the protein model, a homopolymer can fold into a β-hairpin. It is shown that the homopolymer can be folded by either attractive interactions between the monomers or as a direct consequence of the entropic cost of forming interfacial hydrogen bonds in the solvent. However, cold denaturation is not observed if the collapse transition is induced by intramolecular attractions. We further find that it is the changes in hydrophobic hydration with decreasing temperature that drive cold unfolding and that the overall process is enthalpically driven, whereas heat denaturation is entropically driven.

Density of states-based molecular simulations

One of the central problems in statistical mechanics is that of finding the density of states of a system. Knowledge of the density of states of a system is equivalent to knowledge of its fundamental equation, from which all thermodynamic quantities can be obtained. Over the past several years molecular simulations have made considerable strides in their ability to determine the density of states of complex fluids and materials. In this review we discuss some of the more promising approaches proposed in the recent literature along with their advantages and limitations.

Accuracy and convergence of the Wang-Landau sampling algorithm

We present estimations of the accuracy and convergence of the Wang-Landau algorithm. Both accuracy and the related length of the Monte Carlo run depend on the modification parameter f and the density of states. The analytical solution obtained for the two-level system was checked numerically on the two-dimensional Ising model. Although the two-level system is a very simple one, it appears that the proposed solution describes the generic features of the Wang-Landau algorithm. The estimations should prove useful in Monte Carlo calculations of protein folding, first-order transitions, and other systems with a rough energy landscape.

Origin of the decoupling of oxygen and silicon dynamics in liquid silica as expressed by its potential energy landscape

The oxygen and silicon dynamics in silica are compared via computer simulations. In agreement with experimental data and previous simulations a decoupling of oxygen and silicon dynamics is observed upon cooling. The origin of this decoupling is studied in the framework of the potential energy landscape. From analysis of the transition features between neighboring superstructures of minima, denoted metabasins, the differences between the oxygen and silicon dynamics can be quantified. The decoupling can be explicitly related to the presence of generalized rotational processes, giving rise to oxygen but not to silicon displacement. Closer analysis of these processes yields important insight into the nature of the potential energy landscape of silica. The physical picture of relaxation processes in silica, obtained in previous work for oxygen dynamics, is consistent with decoupling effects, elucidated here.

Performances of Wang-Landau algorithms for continuous systems

The relative performances of different implementations of the Wang-Landau method are assessed on two classes of systems with continuous degrees of freedom, namely, two polypeptides and two atomic Lennard-Jones clusters. Parallel tempering Monte Carlo simulations serve as a reference, and we pay particular attention to the variations of the multiplicative factor f during the course of the simulation. For the systems studied, the Wang-Landau method is found to be of comparable accuracy as parallel tempering, but has significant difficulties in reproducing low-temperature transitions exhibited by the Lennard-Jones clusters at low temperature. Using a complementary order parameter and calculating a two-dimensional joint density of states significantly improves the situation, especially for the notoriously difficult LJ(38) system. However, while parallel tempering easily converges for LJ(31), we have not been able to get data of comparable accuracy with Wang-Landau multicanonical sampling.

Efficient combination of Wang–Landau and transition matrix Monte Carlo methods for protein simulations

An efficient combination of the Wang-Landau and transition matrix Monte Carlo methods for protein and peptide simulations is described. At the initial stage of simulation the algorithm behaves like the Wang-Landau algorithm, allowing to sample the entire interval of energies, and at the later stages, it behaves like transition matrix Monte Carlo method and has significantly lower statistical errors. This combination allows to achieve fast convergence to the correct values of density of states. We propose that the violation of TTT identities may serve as a qualitative criterion to check the convergence of density of states. The simulation process can be parallelized by cutting the entire interval of simulation into subintervals. The violation of ergodicity in this case is discussed. We test the algorithm on a set of peptides of different lengths and observe good statistical convergent properties for the density of states. We believe that the method is of general nature and can be used for simulations of other systems with either discrete or continuous energy spectrum.

Computational characterization of the sequence landscape in simple protein alphabets

We characterize the "sequence landscapes" in several simple, heteropolymer models of proteins by examining their mutation properties. Using an efficient flat-histogram Monte Carlo search method, our approach involves determining the distribution in energy of all sequences of a given length when threaded through a common backbone. These calculations are performed for a number of Protein Data Bank structures using two variants of the 20-letter contact potential developed by Miyazawa and Jernigan [Miyazawa S, Jernigan WL. Macromolecules 1985;18:534], and the 2-monomer HP model of Lau and Dill [Lau KF, Dill KA. Macromolecules 1989;22:3986]. Our results indicate significant differences among the energy functions in terms of the "smoothness" of their landscapes. In particular, one of the Miyazawa-Jernigan contact potentials reveals unusual cooperative behavior among its species' interactions, resulting in what is essentially a set of phase transitions in sequence space. Our calculations suggest that model-specific features can have a profound effect on protein design algorithms, and our methods offer a number of ways by which sequence landscapes can be quantified.

A unified methodological framework for the simulation of nonisothermal ensembles

A general framework is developed for the simulation of nonisothermal statistical-mechanical ensembles. This framework is intended to synthesize the formulation of advanced Monte Carlo simulation methods such as multihistogram reweighting, replica-exchange methods, and expanded ensemble techniques so that they can be applied to different nonisothermal ensembles. Using Lennard-Jones systems as test cases, novel implementations of these methods are demonstrated with different ensembles including the microcanonical, isobaric-isoenthalpic, and isobaric-semigrand ensembles. In particular, it is shown that the use of multiensemble methods allows the efficient simulation of microcanonical density of states, entropies, vapor-liquid and solid-liquid equilibrium for pure component systems, and fluid-phase coexistence for binary mixtures. In these applications, comparisons are also presented that highlight the advantages of the proposed multiensemble implementations over alternative methods used before.