A study of the multi-canonical Monte Carlo method

Diabat method for polymorph free energies: Extension to molecular crystals

Lattice-switch Monte Carlo and the related diabat methods have emerged as efficient and accurate ways to compute free energy differences between polymorphs. In this work, we introduce a one-to-one mapping from the reference positions and displacements in one molecular crystal to the positions and displacements in another. Two features of the mapping facilitate lattice-switch Monte Carlo and related diabat methods for computing polymorph free energy differences. First, the mapping is unitary so that its Jacobian does not complicate the free energy calculations. Second, the mapping is easily implemented for molecular crystals of arbitrary complexity. We demonstrate the mapping by computing free energy differences between polymorphs of benzene and carbamazepine. Free energy calculations for thermodynamic cycles, each involving three independently computed polymorph free energy differences, all return to the starting free energy with a high degree of precision. The calculations thus provide a force field independent validation of the method and allow us to estimate the precision of the individual free energy differences.

Gibbs free-energy differences between polymorphs via a diabat approach

Polymorph free-energy differences are critical to several applications. A recently proposed diabat interpolation framework estimated free-energy differences between polymorphs by quadratic interpolation of diabats. This work extends the Zwanzig-Bennett relation to the NPT ensemble so that the diabats directly give Gibbs free-energy differences. We also demonstrate how the approach can be used in cases where the diabats are not parabolic. We illustrate the diabat method for Gibbs free-energy difference of zirconium (BCC and HCP phases) and compare it with the conventional lattice switch Monte Carlo approach.

Improved grand canonical sampling of vapour-liquid transitions

Simulation within the grand canonical ensemble is the method of choice for accurate studies of first order vapour-liquid phase transitions in model fluids. Such simulations typically employ sampling that is biased with respect to the overall number density in order to overcome the free energy barrier associated with mixed phase states. However, at low temperature and for large system size, this approach suffers a drastic slowing down in sampling efficiency. The culprits are geometrically induced transitions (stemming from the periodic boundary conditions) which involve changes in droplet shape from sphere to cylinder and cylinder to slab. Since the overall number density does not discriminate sufficiently between these shapes, it fails as an order parameter for biasing through the transitions. Here we report two approaches to ameliorating these difficulties. The first introduces a droplet shape based order parameter that generates a transition path from vapour to slab states for which spherical and cylindrical droplets are suppressed. The second simply biases with respect to the number density in a tetragonal subvolume of the system. Compared to the standard approach, both methods offer improved sampling, allowing estimates of coexistence parameters and vapor-liquid surface tension for larger system sizes and lower temperatures.

Renormalized multicanonical sampling in multiple dimensions

This paper extends to two and three dimensions the recently proposed renormalized multicanonical sampling procedure [D. Yevick, Int. J. Mod. Phys. C 27, 1650033 (2016)IJMPEO0129-183110.1142/S0129183116500339], which generates the density of states of a system from that of a smaller system. The accuracy and efficiency of the method is additionally enhanced through a transition matrix Monte Carlo method that incorporates multicanonical sampling. The performance of the resulting procedure is investigated through calculations of the specific heat of the standard Ising model for square lattices with up to 128^{2} spins and for cubic lattices with up to 24^{3} spins. A comparison of the scaling properties to those of other algorithms based on the transition matrix Monte Carlo method demonstrates that the relative advantage of the renormalized procedure increases rapidly with system size.

Continuous Fractional Component Monte Carlo: An Adaptive Biasing Method for Open System Atomistic Simulations

A new open system Monte Carlo procedure designed to overcome difficulties with insertion and deletion of molecules is introduced. The method utilizes gradual insertions and deletions of molecules through the use of a continuous coupling parameter and an adaptive bias potential. The method draws upon concepts from previous open system molecular dynamics and expanded ensemble Monte Carlo techniques and is applied to both the grand canonical and osmotic ensembles. It is shown to yield correct results for the volumetric properties of the Lennard-Jones fluid and water as well as the phase behavior of the CO2-ethanol binary system.

Quantifying Density Fluctuations in Water at a Hydrophobic Surface: Evidence for Critical Drying

Employing smart Monte Carlo sampling techniques within the grand canonical ensemble, we investigate the properties of water at a model hydrophobic substrate. By reducing the strength of substrate-water attraction, we find that fluctuations in the local number density, quantified by a rigorous definition of the local compressibility χ(z), increase rapidly for distances z within one or two molecular diameters from the substrate as the degree of hydrophobicity, measured by the macroscopic contact angle θ, increases. Our simulations provide evidence for a continuous (critical) drying transition as the substrate-water interaction becomes very weak: cos(θ)→-1. We speculate that the existence of such a transition might account for earlier simulation observations of strongly enhanced density fluctuations.

Relative free enthalpies for point mutations in two proteins with highly similar sequences but different folds

Enveloping distribution sampling was used to calculate free-enthalpy changes associated with single amino acid mutations for a pair of proteins, GA95 and GB95, that show 95% sequence identity yet fold into topologically different structures. Of the L → A, I → F, and L → Y mutations at positions 20, 30, and 45, respectively, of the 56-residue sequence, the first and the last contribute the most to the free-enthalpy difference between the native and non-native sequence-structure combinations, in agreement with the experimental findings for this protein pair. The individual free-enthalpy changes are almost sequence-independent in the four-strand/one-helix structure, the stable form of GB95, while in the three-helix bundle structure, the stable form of GA95, an interplay between residues 20 and 45 is observed.

Transitions between imperfectly ordered crystalline structures: a phase switch Monte Carlo study

A model for two-dimensional colloids confined laterally by "structured boundaries" (i.e., ones that impose a periodicity along the slit) is studied by Monte Carlo simulations. When the distance D between the confining walls is reduced at constant particle number from an initial value D(0), for which a crystalline structure commensurate with the imposed periodicity fits, to smaller values, a succession of phase transitions to imperfectly ordered structures occur. These structures have a reduced number of rows parallel to the boundaries (from n to n-1 to n-2, etc.) and are accompanied by an almost periodic strain pattern, due to "soliton staircases" along the boundaries. Since standard simulation studies of such transitions are hampered by huge hysteresis effects, we apply the phase switch Monte Carlo method to estimate the free energy difference between the structures as a function of the misfit between D and D(0), thereby locating where the transitions occur in equilibrium. For comparison, we also obtain this free energy difference from a thermodynamic integration method: The results agree, but the effort required to obtain the same accuracy as provided by phase switch Monte Carlo would be at least three orders of magnitude larger. We also show for a situation where several "candidate structures" exist for a phase, that phase switch Monte Carlo can clearly distinguish the metastable structures from the stable one. Finally, applying the method in the conjugate statistical ensemble (where the normal pressure conjugate to D is taken as an independent control variable), we show that the standard equivalence between the conjugate ensembles of statistical mechanics is violated.

Thermodynamic stability of the cubatic phase of hard cut spheres evaluated by expanded ensemble simulations

The system of hard cut spheres (disk-shaped particles formed by symmetrically truncating the end caps of a sphere) exhibits an intriguing "cubatic" phase with cubic orientational symmetry. However, it is unclear whether this phase is metastable with respect to the columnar phase. We attempt to provide an answer to this question by carrying out free energy calculations by the expanded ensemble Monte Carlo method. We conclude that there may be a very small region of cubatic stability in the vicinity of the isotropic-cubatic phase transition, but that that transition would need to be determined more accurately to obtain a definitive answer. We also comment on the efficacy of the expanded ensemble method for these kinds of calculations.

Excess entropy scaling of dynamic quantities for fluids of dumbbell-shaped particles

We use molecular simulation to study the ability of entropy scaling relationships to describe the kinetic properties of two Lennard-Jones dumbbell models. We begin by examining the excess entropy, the key quantity used to correlate dynamic properties within entropy scaling strategies. We compute the thermodynamic excess entropy as well as contributions to the two-body excess entropy stemming from translational and orientational intermolecular correlations. Our results indicate that the total two-body contribution accounts for more than 70% of the thermodynamic excess entropy at all state conditions explored. For the two dumbbell models studied here, the orientational component of the two-body excess entropy dominates at moderate and high fluid densities. We next investigate the relationships between kinetic properties and various contributions to the excess entropy. Four dynamic properties are considered: translational and rotational diffusivities, a characteristic relaxation time for rotational motion, and a collective relaxation time stemming from analysis of the coherent intermediate-scattering function. We find that the thermodynamic excess entropy provides the best metric for describing kinetic properties. For each of the dynamic properties considered, reduced data collapse onto a common curve when expressed as a function of the thermodynamic excess entropy. The likelihood of a two-body contribution to the excess entropy serving as a reliable scaling variable is linked to the extent to which it correlates with the thermodynamic excess entropy. The total two-body term contributes significantly to the excess entropy, and therefore this quantity generally serves as a suitable scaling variable.

Accelerating simulation of metastable decay

We present a new adaptation of umbrella sampling technique in studying nucleation, allowing for a large free energy barrier to be overcome within a time scale accessible to simulation. The method requires a set of order parameters that can distinguish a metastable state and a more stable state toward which the system relaxes, but does not demand that the order parameters be adequate for describing the process around the top of the free energy barrier. To illustrate its effectiveness, the method was implemented using Monte Carlo simulation and was applied to study bubble nucleation in single component truncated and shifted Lennard-Jones fluids. Under the conditions studied, the rate of nucleation, measured in Monte Carlo cycle, was enhanced by a factor ranging from 10(2) to 10(7). The simulation also yields various thermodynamic and structural information on the system during nucleation.

Alignment methods for biased multicanonical sampling

The efficiency of the multicanonical procedure can be significantly improved by applying an additional bias to the numerically generated sample space. However, results obtained by biasing in different sampling regions cannot in general be accurately combined, since their relative normalization coefficient is not known precisely. We demonstrate that for overlapping biasing regions a simple iterative procedure can be employed to determine the required coefficients.

Examining the accuracy of ideal adsorbed solution theory without curve-fitting using transition matrix Monte Carlo simulations

Ideal adsorbed solution theory (IAST) is a well-known approach to predicting multicomponent adsorption isotherms in microporous materials from experimental or simulation data for single-component adsorption. A limitation in practical applications of IAST is that useful calculations often require extrapolation of fitted single-component isotherms beyond the range for which data are available. We introduce a molecular simulation approach in which the intrinsic accuracy of IAST can be examined in a context that avoids any need to perform curve fitting with single-component data. Our approach is based on using transition matrix Monte Carlo to define single-component adsorption isotherms for arbitrary bulk-phase pressures from a single simulation. We apply our approach to several light gas mixtures in silica zeolites and a carbon nanotube to examine the intrinsic accuracy of IAST for these model systems.

Bayesian update method for adaptive weighted sampling

Exploring conformational spaces is still a challenging task for simulations of complex systems. One way to enhance such a task is weighted sampling, e.g., by assigning high weights to regions that are rarely sampled. It is, however, difficult to estimate adequate weights beforehand, and therefore adaptive methods are desired. Here we present a method for adaptive weighted sampling based on Bayesian inference. Within the framework of Bayesian inference, we develop an update scheme in which the information from previous data is stored in a prior distribution which is then updated to a posterior distribution according to new data. The method proposed here is particularly well suited for distributed computing, in which one must deal with rapid influxes of large amounts of data.

Accurate estimation of the density of states from Monte Carlo transition probability data

This study develops an efficient approach for calculating the density of states from energy transition probability matrices generated from extended sampling Monte Carlo simulations. Direct and iterative variants of the method are shown to achieve high accuracy when applied to the two-dimensional Ising model for which the density of states function can be determined exactly. They are also used to calculate the density of states of lattice protein and Lennard-Jones models which generate more complex nonzero matrix structures. Whereas the protein simulations test the method on a system exhibiting a rugged free energy landscape, the Lennard-Jones calculations highlight implementation details that arise in applications to continuous energy systems. Density of states results for these two systems agree with estimates from multiple histogram reweighting, demonstrating that the new method provides an alternative approach for computing the thermodynamic properties of complex systems.

On the use of Bennett's acceptance ratio method in multi-canonical-type simulations

A common strategy for mapping coexistence curves is to employ multi-canonical (MUCA) sampling to simulate along a macrostate path connecting two phases. Central to this approach is the task of accurately calculating the importance weights used in the MUCA procedure, which are needed for both effective sampling and accurate determination of phase boundaries. The purpose of this study is to develop a strategy for determining the importance weights that is built upon Bennett's optimized acceptance ratio method. This approach is shown to be closely related to transition matrix schemes, and is used to compute the vapor-liquid equilibrium of a Lennard-Jones fluid and the liquid-liquid equilibrium of a n-hexane/n-perfluorohexane mixture. For the Lennard-Jones system, the importance weights as a function of the number of particles "N" (at fixed temperature and volume) are obtained by using Bennett's method to estimate free energy differences between N and N+1 particle systems over the desired range of N values. In this application, the method is found to perform slightly better than a related transition matrix scheme. For the n-hexane/n-perfluorohexane liquid mixture, the method is designed to obtain weights as a function of composition (for fixed temperature, pressure, and total number of particles); in this case, the method is found to outperform the Gibbs ensemble approach.

Universal features and tail analysis of the order-parameter distribution of the two-dimensional Ising model: an entropic sampling Monte Carlo study

We present a numerical study of the order-parameter probability density function (PDF) of the square Ising model for lattices with linear sizes L = 80 - 140. A recent efficient entropic sampling scheme, combining the Wang-Landau and broad histogram methods and based on the high levels of the Wang-Landau process in dominant energy subspaces is employed. We find that for large lattices there exists a stable window of the scaled order-parameter in which the full ansatz including the pre-exponential factor for the tail regime of the universal PDF is well obeyed. This window is used to estimate the equation of state exponent and to observe the behavior of the universal constants implicit in the functional form of the universal PDF. The probability densities are used to estimate the universal Privman-Fisher coefficient and to investigate whether one could obtain reliable estimates of the universal constants controlling the asymptotic behavior of the tail regime.

Simulation of the density of states in isothermal and adiabatic ensembles

This paper provides a unified treatment of the fundamental methods used to obtain the density of states via molecular simulations with isothermal ensembles (IEs) and adiabatic ensembles (AEs). Our analysis and results show that provides a natural bridge to go back and forth between IE and AE simulation data. They also underline the difference between the density of states of potential energy macrostates and that of total energy macrostates Omega, even though both provide access to the thermodynamic properties of the system. Visited-states approaches and transition matrix methods are described and applied to the Lennard-Jones fluid to target omega and Omega as functions of energy and volume macrostates. It is shown that one can obtain omega via a generalized acceptance-ratio formula that is applicable regardless of the conditions at which the ensemble is simulated. In this way, one can obtain while performing conventional IE or AE simulations, and do it at no extra cost and with a higher accuracy than is achievable with histogram methods.

On the use of transition matrix methods with extended ensembles

Different extended ensemble schemes for non-Boltzmann sampling (NBS) of a selected reaction coordinate lambda were formulated so that they employ (i) "variable" sampling window schemes (that include the "successive umbrella sampling" method) to comprehensibly explore the lambda domain and (ii) transition matrix methods to iteratively obtain the underlying free-energy eta landscape (or "importance" weights) associated with lambda. The connection between "acceptance ratio" and transition matrix methods was first established to form the basis of the approach for estimating eta(lambda). The validity and performance of the different NBS schemes were then assessed using as lambda coordinate the configurational energy of the Lennard-Jones fluid. For the cases studied, it was found that the convergence rate in the estimation of eta is little affected by the use of data from high-order transitions, while it is noticeably improved by the use of a broader window of sampling in the variable window methods. Finally, it is shown how an "elastic" window of sampling can be used to effectively enact (nonuniform) preferential sampling over the lambda domain, and how to stitch the weights from separate one-dimensional NBS runs to produce a eta surface over a two-dimensional domain.