Determination of Free Energy from Chemical Potentials: Application of the Expanded Ensemble Method

Alchemical Free-Energy Calculations by Multiple-Replica λ-Dynamics: The Conveyor Belt Thermodynamic Integration Scheme

A new method is proposed to calculate alchemical free-energy differences based on molecular dynamics (MD) simulations, called the conveyor belt thermodynamic integration (CBTI) scheme. As in thermodynamic integration (TI), K replicas of the system are simulated at different values of the alchemical coupling parameter λ. The number K is taken to be even, and the replicas are equally spaced on a forward-turn-backward-turn path, akin to a conveyor belt (CB) between the two physical end-states; and as in λ-dynamics (λD), the λ-values associated with the individual systems evolve in time along the simulation. However, they do so in a concerted fashion, determined by the evolution of a single dynamical variable Λ of period 2π controlling the advance of the entire CB. Thus, a change of Λ is always associated with K/2 equispaced replicas moving forward and K/2 equispaced replicas moving backward along λ. As a result, the effective free-energy profile of the replica system along Λ is periodic of period 2 πK^-1, and the magnitude of its variations decreases rapidly upon increasing K, at least as K^-1 in the limit of large K. When a sufficient number of replicas is used, these variations become small, which enables a complete and quasi-homogeneous coverage of the λ-range by the replica system, without application of any biasing potential. If desired, a memory-based biasing potential can still be added to further homogenize the sampling, the preoptimization of which is computationally inexpensive. The final free-energy profile along λ is calculated similarly to TI, by binning of the Hamiltonian λ-derivative as a function of λ considering all replicas simultaneously, followed by quadrature integration. The associated quadrature error can be kept very low owing to the continuous and quasi-homogeneous λ-sampling. The CBTI scheme can be viewed as a continuous/deterministic/dynamical analog of the Hamiltonian replica-exchange/permutation (HRE/HRP) schemes or as a correlated multiple-replica analog of the λD or λ-local elevation umbrella sampling (λ-LEUS) schemes. Compared to TI, it shares the advantage of the latter schemes in terms of enhanced orthogonal sampling, i.e. the availability of variable-λ paths to circumvent conformational barriers present at specific λ-values. Compared to HRE/HRP, it permits a deterministic and continuous sampling of the λ-range, is expected to be less sensitive to possible artifacts of the thermo- and barostating schemes, and bypasses the need to carefully preselect a λ-ladder and a swapping-attempt frequency. Compared to λ-LEUS, it eliminates (or drastically reduces) the dead time associated with the preoptimization of a biasing potential. The goal of this article is to provide the mathematical/physical formulation of the proposed CBTI scheme, along with an initial application of the method to the calculation of the hydration free energy of methanol.

A method for computing the solubility limit of solids: application to sodium chloride in water and alcohols

We present an adaptable method to compute the solubility limit of solids by molecular simulation, which avoids the difficulty of reference state calculations. In this way, the method is highly adaptable to molecules of complex topology. Results are shown for solubility calculations of sodium chloride in water and light alcohols at atmospheric conditions. The pseudosupercritical path integration method is used to calculate the free energy of the solid and gives results that are in good agreement with previous studies that reference the Einstein crystal. For the solution phase calculations, the self-adaptive Wang-Landau transition-matrix Monte Carlo method is used within the context of an expanded isothermal-isobaric ensemble. The method shows rapid convergence properties and the uncertainty in the calculated chemical potential was 1% or less for all cases. The present study underpredicts the solubility limit of sodium chloride in water, suggesting a shortcoming of the molecular models. Importantly, the proper trend for the chemical potential in various solvents was captured, suggesting that relative solubilities can be computed by the method.

A generalized direct-particle-deletion scheme for the calculation of chemical potential and solubilities of small- and medium-sized molecules in amorphous polymers

The development, validation, and first applications of a generalized version of an inverse Widom method are described. It permits the calculation of solubility coefficients for molecules as large as, e.g., benzene in all polymers for which reasonable forcefield parameters exist. Predicting the solubility is a key to the knowledge-based design of materials utilized to solve permeability related problems. For long time, particle insertion methods, such as the Widom method, were the only way to predict solubilities from molecular models, but they, in most cases, only worked well for rather small penetrants (e.g., H2, O2, N2). Therefore, a few years ago, a new particle deletion algorithm "DPD" was introduced by Boulougouris, Economou, and Theodorou to overcome this problem in principle. The related computer code was, however, only applicable to special, relatively simple model systems. As application examples for the generalized version described here, solubility calculations for nitrogen, oxygen, and benzene in poly(dimethyl siloxane) are presented.

Toward efficient chemical potential calculations by expanded ensemble simulations; to make the free energy pathway fairly level

A scheme is suggested of how to construct good bias potentials ("balancing factors") to be used in expanded ensemble (EE) calculations of chemical potentials of solutions. A combination of two strategies are used: (i) to use a pathway for particle insertions that avoids large variations in free energy and (ii) to use calculated free energy derivatives to construct a bias potential that makes the pathway fairly level. Only a few very short simulations are needed to accomplish the latter, and then, a full EE simulation is done to obtain the chemical potential. By practical calculations of the chemical potential of benzene, cyclohexane, and benzylamine in water, it is shown that this method is at least equally efficient to the recent adaptive EE (AEE) method by Aberg et al. (J. Chem. Phys. 2004, 120, 3370). Furthermore, the new method provides an alternative strategy that complements existing EE methods.

Monte Carlo Simulations of Gas Solubility in the Ionic Liquid 1-n-Butyl-3-methylimidazolium Hexafluorophosphate

The Henry's constants of water, carbon dioxide, ethane, ethene, methane, oxygen, and nitrogen are computed in the ionic liquid 1-n-butyl-3-methylimidazolium hexafluorophosphate ([bmim][PF(6)]) using test particle insertion and expanded ensemble Monte Carlo methods. The partial molar enthalpy and partial molar entropy of solvation are also computed for water, carbon dioxide, and oxygen. The results from the simulations are compared against experimental data from the literature. In addition, the accuracy and precision of the two methods in determining the Henry's constant are examined. Local organization of the ionic liquid around a solute molecule is analyzed, and the interactions responsible for the experimentally observed solubility trends are identified.

Isomolar semigrand ensemble molecular dynamics: Development and application to liquid-liquid equilibria

An extended system molecular dynamics method for the isomolar semigrand ensemble (fixed number of particles, pressure, temperature, and fugacity fraction) is developed and applied to the calculation of liquid-liquid equilibria (LLE) for two Lennard-Jones mixtures. The method utilizes an extended system variable to dynamically control the fugacity fraction xi of the mixture by gradually transforming the identity of particles in the system. Two approaches are used to compute coexistence points. The first approach uses multiple-histogram reweighting techniques to determine the coexistence xi and compositions of each phase at temperatures near the upper critical solution temperature. The second approach, useful for cases in which there is no critical solution temperature, is based on principles of small system thermodynamics. In this case a coexistence point is found by running N-P-T-xi simulations at a common temperature and pressure and varying the fugacity fraction to map out the difference in chemical potential between the two species A and B (mu(A)-mu(B)) as a function of composition. Once this curve is known the equal-distance/equal-area criterion is used to determine the coexistence point. Both approaches give results that are comparable to those of previous Monte Carlo (MC) simulations. By formulating this approach in a molecular dynamics framework, it should be easier to compute the LLE of complex molecules whose intramolecular degrees of freedom are often difficult to properly sample with MC techniques.