On the use of transition matrix methods with extended ensembles

On the calculation of free energies over Hamiltonian and order parameters via perturbation and thermodynamic integration

In this work, complementary formulas are presented to compute free-energy differences via perturbation (FEP) methods and thermodynamic integration (TI). These formulas are derived by selecting only the most statistically significant data from the information extractable from the simulated points involved. On the one hand, commonly used FEP techniques based on overlap sampling leverage the full information contained in the overlapping macrostate probability distributions. On the other hand, conventional TI methods only use information on the first moments of those distributions, as embodied by the first derivatives of the free energy. Since the accuracy of simulation data degrades considerably for high-order moments (for FEP) or free-energy derivatives (for TI), it is proposed to consider, consistently for both methods, data up to second-order moments/derivatives. This provides a compromise between the limiting strategies embodied by common FEP and TI and leads to simple, optimized expressions to evaluate free-energy differences. The proposed formulas are validated with an analytically solvable harmonic Hamiltonian (for assessing systematic errors), an atomistic system (for computing the potential of mean force with coordinate-dependent order parameters), and a binary-component coarse-grained model (for tracing a solid-liquid phase diagram in an ensemble sampled through alchemical transformations). It is shown that the proposed FEP and TI formulas are straightforward to implement, perform similarly well, and allow robust estimation of free-energy differences even when the spacing of successive points does not guarantee them to have proper overlapping in phase space.

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.

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.

Understanding free-energy perturbation calculations through a model of harmonic oscillators: theory and implications to improve the sampling efficiency by molecular simulation

Free-energy perturbation calculation is frequently used to calculate free-energy differences because it is easy to implement and the computation is fast. However, the calculation is subject to large inaccuracies in some circumstances due to the insufficient sampling of the relevant tails of the energy-difference distributions. Here we expand this knowledge of insufficient sampling into a two-dimensional (2D) energy space using a model of harmonic oscillators. We show analytically the relation between the energies of the sampling system and those of the desired target energy spaces, which provide the basis to understand the difficulties in free-energy perturbation calculations. We clarify the reasons of the inaccurate calculation in the different harmonic cases that stem from the spatial separations of the reference and the target energy pairs located in the two-dimensional energy space. The potential-energy space introduced into this 2D energy-space model provides additional clues to improve the sampling efficiency. Based on this understanding, we propose two ways to calculate the free-energy differences using the two schemes of the distribution method. We show that the distribution method implemented in the appropriate energy space--the energy-difference space and the potential-energy space, respectively--can improve the calculation of free energies in different circumstances. This analysis implies that the sampling can be improved if it is directed toward the appropriate region in the potential-energy space, which is easily implemented in various types of free-energy calculations. To test this, we calculate the free-energy surface of alanine dipeptide in gas phase and in aqueous phase, respectively. We demonstrate that the free-energy surface calculation is improved when the biased sampling of the potential energy is integrated into the sampling scheme.

Optimized expanded ensembles for simulations involving molecular insertions and deletions. II. Open systems

In the Grand Canonical, osmotic, and Gibbs ensembles, chemical potential equilibrium is attained via transfers of molecules between the system and either a reservoir or another subsystem. In this work, the expanded ensemble (EXE) methods described in part I [F. A. Escobedo and F. J. Martinez-Veracoechea, J. Chem. Phys. 127, 174103 (2007)] of this series are extended to these ensembles to overcome the difficulties associated with implementing such whole-molecule transfers. In EXE, such moves occur via a target molecule that undergoes transitions through a number of intermediate coupling states. To minimize the tunneling time between the fully coupled and fully decoupled states, the intermediate states could be either: (i) sampled with an optimal frequency distribution (the sampling problem) or (ii) selected with an optimal spacing distribution (staging problem). The sampling issue is addressed by determining the biasing weights that would allow generating an optimal ensemble; discretized versions of this algorithm (well suited for small number of coupling stages) are also presented. The staging problem is addressed by selecting the intermediate stages in such a way that a flat histogram is the optimized ensemble. The validity of the advocated methods is demonstrated by their application to two model problems, the solvation of large hard spheres into a fluid of small and large spheres, and the vapor-liquid equilibrium of a chain system.

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.