Widom Particle Insertion Method Research Articles

Boundary-Monte Carlo Method for Neutral and Charged Confined Fluids

In this work, we describe a new Monte Carlo (MC) simulation method to investigate highly coupled fluids in confined geometries at a constant chemical potential. This method is based on so-called multi-scale Hamiltonian methods, wherein the chemical potential is determined using a more amenable Hamiltonian for a fluid in an "outer" region, which facilitates standard methods, such as grand canonical MC simulations or Widom's particle insertion method. The (inner region) fluid of interest is placed in diffusive contact with the simpler outer fluid via a boundary zone wherein the Hamiltonian is transformed. The current method utilizes an ideal fluid for the outer regions, which allows for implicit rather than explicit simulations. Only the boundary and inner region need explicit consideration; hence, the nomenclature used is boundary-Monte Carlo. We illustrate the utility of the method for simple neutral and charged fluids in cylindrical and planar pores. In the latter case, we use a dense room-temperature ionic liquid model and illustrate how the boundary zone establishes a proper Donnan equilibrium between inner and outer fluids in the presence of charged planar electrodes. Thus, the method allows direct calculation of properties such as the differential capacitance, without the need for additional difficult calculations of the requisite Donnan potential.

The Gibbs free energy of cavity formation in a diverse set of solvents.

The concept of the formation of a solute-sized cavity in a solvent is widely used in the theories of solvation processes; however, most of the studies of cavity formation using atomistic simulations were limited to water and hydrocarbon models. We calculated the Gibbs free energy of cavity formation ΔcavG for a structurally diverse set of 23 common organic solvents. For the calculation, molecular dynamics simulations of solvent boxes were conducted, and the Widom particle insertion method was applied. The results obtained with two different force fields for the same solvent were in good agreement with each other in most cases. The obtained cavity size dependences of ΔcavG allowed ranking the solvents by the free energy cost of creation of a cavity with a certain size. Surprisingly, this cost was somewhat higher in glycerol, formamide, and ethylene glycol than in water. In general, higher values of ΔcavG are observed for the solvents with a branched network of intermolecular hydrogen bonds and strongly polar aprotic solvents. The numerical results can be used to improve the accuracy of the calculation of the cavity term in non-aqueous continuum solvation models.

Local Grand Canonical Monte Carlo Simulation Method for Confined Fluids.

We describe a new local grand canonical Monte Carlo method to treat fluids in pores in chemical equilibrium with a reference bulk. The method is applied to Lennard-Jones particles in pores of different geometry and is shown to be much more accurate and efficient than other techniques such as traditional grand canonical simulations or Widom's particle insertion method. It utilizes a penalty potential to create a gas phase, which is in equilibrium with a more dense liquid component in the pore. Grand canonical Monte Carlo moves are employed in the gas phase, and the system then maintains chemical equilibrium by "diffusion" of particles. This creates an interface, which means that the confined fluid needs to occupy a large enough volume so that this is not an issue. We also applied the method to confined charged fluids and show how it can be used to determine local electrostatic potentials in the confined fluid, which are properly referenced to the bulk. This precludes the need to determine the Donnan potential (which controls electrochemical equilibrium) explicitly. Prior approaches have used explicit bulk simulations to measure this potential difference, which are significantly costly from a computational point of view. One outcome of our analysis is that pores of finite cross-section create a potential difference with the bulk via a small but nonzero linear charge density, which diminishes as ∼1/ln(L), where L is the pore length.

Chemical potential of a test hard sphere of variable size in hard-sphere fluid mixtures.

A detailed comparison between the Boublík-Mansoori-Carnahan-Starling-Leland (BMCSL) equation of state of hard-sphere mixtures is made with Molecular Dynamics (MD) simulations of the same compositions. The Labík and Smith simulation technique [S. Labík and W. R. Smith, Mol. Simul. 12, 23-31 (1994)] was used to implement the Widom particle insertion method to calculate the excess chemical potential, βμ0ex, of a test particle of variable diameter, σ0, immersed in a hard-sphere fluid mixture with different compositions and values of the packing fraction, η. Use is made of the fact that the only polynomial representation of βμ0ex which is consistent with the limits σ0 → 0 and σ0 → ∞ has to be of the cubic form, i.e., c0(η)+c¯1(η)σ0/M1+c¯2(η)(σ0/M1)2+c¯3(η)(σ0/M1)3, where M1 is the first moment of the distribution. The first two coefficients, c0(η) and c¯1(η), are known analytically, while c¯2(η) and c¯3(η) were obtained by fitting the MD data to this expression. This in turn provides a method to determine the excess free energy per particle, βaex, in terms of c¯2, c¯3, and the compressibility factor, Z. Very good agreement between the BMCSL formulas and the MD data is found for βμ0ex, Z, and βaex for binary mixtures and continuous particle size distributions with the top-hat analytic form. However, the BMCSL theory typically slightly underestimates the simulation values, especially for Z, differences which the Boublík-Carnahan-Starling-Kolafa formulas and an interpolation between two Percus-Yevick routes capture well in different ranges of the system parameter space.

CHEMICAL POTENTIALS OF HARD-CORE MOLECULES BY A STEPWISE INSERTION METHOD

ABSTRACT A molecular simulation algorithm was implemented to calculate chemical potentials of hard-core molecular systems at high densities. The method is based on the Widom particle insertion method and the step-function character of free energy variations. The algorithm was evaluated for hard-sphere mixtures at infinite dilution approximation by varying the solute/solvent diameter ratio, for systems with reduced densities from 0.1 to 0.8. The proposed methodology was verified by comparing simulations of trimers diluted in spheres and of single-component dimer systems with results from the literature. Then, the method was applied to mixtures of hard-spheres and dimers at several conditions regarding composition, reduced density, and bond-length/diameter ratio. The results were used to validate equations of state from the literature. The proposed approach was able to obtain accurate chemical potentials for different hard-core molecular mixtures. Lower uncertainties were obtained when comparing with traditional methods, especially at high densities.

Equation of state and Helmholtz free energy for the atomic system of the repulsive Lennard-Jones particles.

Simple and accurate expressions are presented for the equation of state (EOS) and absolute Helmholtz free energy of a system composed of simple atomic particles interacting through the repulsive Lennard-Jones potential model in the fluid and solid phases. The introduced EOS has 17 and 22 coefficients for fluid and solid phases, respectively, which are regressed to the Monte Carlo (MC) simulation data over the reduced temperature range of 0.6≤T*≤6.0 and the packing fraction range of 0.1 ≤ η ≤ 0.72. The average absolute relative percent deviation in fitting the EOS parameters to the MC data is 0.06 and 0.14 for the fluid and solid phases, respectively. The thermodynamic integration method is used to calculate the free energy using the MC simulation results. The Helmholtz free energy of the ideal gas is employed as the reference state for the fluid phase. For the solid phase, the values of the free energy at the reduced density equivalent to the close-packed of a hard sphere are used as the reference state. To check the validity of the predicted values of the Helmholtz free energy, the Widom particle insertion method and the Einstein crystal technique of Frenkel and Ladd are employed. The results obtained from the MC simulation approaches are well agreed to the EOS results, which show that the proposed model can reliably be utilized in the framework of thermodynamic theories.

Chemical potential of a test hard sphere of variable size in a hard-sphere fluid.

The Labík and Smith Monte Carlo simulation technique to implement the Widom particle insertion method is applied using Molecular Dynamics (MD) instead to calculate numerically the insertion probability, P0(η,σ0), of tracer hard-sphere (HS) particles of different diameters, σ0, in a host HS fluid of diameter σ and packing fraction, η, up to 0.5. It is shown analytically that the only polynomial representation of -ln⁡P0(η,σ0) consistent with the limits σ0→0 and σ0→∞ has necessarily a cubic form, c0(η)+c1(η)σ0/σ+c2(η)(σ0/σ)2+c3(η)(σ0/σ)3. Our MD data for -ln⁡P0(η,σ0) are fitted to such a cubic polynomial and the functions c0(η) and c1(η) are found to be statistically indistinguishable from their exact solution forms. Similarly, c2(η) and c3(η) agree very well with the Boublík-Mansoori-Carnahan-Starling-Leland and Boublík-Carnahan-Starling-Kolafa formulas. The cubic polynomial is extrapolated (high density) or interpolated (low density) to obtain the chemical potential of the host fluid, or σ0→σ, as βμex=c0+c1+c2+c3. Excellent agreement between the Carnahan-Starling and Carnahan-Starling-Kolafa theories with our MD data is evident.

Chemical potential and entropy in monodisperse and polydisperse hard-sphere fluids using Widom's particle insertion method and a pore size distribution-based insertion probability.

We estimate the excess chemical potential Δμ and excess entropy per particle Δs of computer-generated, monodisperse and polydisperse, frictionless hard-sphere fluids. For this purpose, we utilize the Widom particle insertion method, which for hard-sphere systems relates Δμ to the probability to successfully (without intersections) insert a particle into a system. This insertion probability is evaluated directly for each configuration of hard spheres by extrapolating to infinity the pore radii (nearest-surface) distribution and integrating its tail. The estimates of Δμ and Δs are compared to (and comply well with) predictions from the Boublík-Mansoori-Carnahan-Starling-Leland equation of state. For polydisperse spheres, we employ log-normal particle radii distributions with polydispersities δ = 0.1, 0.2, and 0.3.

Upper bound on the Edwards entropy in frictional monodisperse hard-sphere packings.

We extend the Widom particle insertion method [B. Widom, J. Chem. Phys., 1963, 39, 2808-2812] to determine an upper bound sub on the Edwards entropy in frictional hard-sphere packings. sub corresponds to the logarithm of the number of mechanically stable configurations for a given volume fraction and boundary conditions. To accomplish this, we extend the method for estimating the particle insertion probability through the pore-size distribution in frictionless packings [V. Baranau, et al., Soft Matter, 2013, 9, 3361-3372] to the case of frictional particles. We use computer-generated and experimentally obtained three-dimensional sphere packings with volume fractions φ in the range 0.551-0.65. We find that sub has a maximum in the vicinity of the Random Loose Packing Limit φRLP = 0.55 and decreases then monotonically with increasing φ to reach a minimum at φ = 0.65. Further on, sub does not distinguish between real mechanical stability and packings in close proximity to mechanical stable configurations. The probability to find a given number of contacts for a particle inserted in a large enough pore does not depend on φ, but it decreases strongly with the contact number.

Crystal nuclei in melts: a Monte Carlo simulation of a model for attractive colloids

As a model for a suspension of hard-sphere-like colloidal particles where small non-adsorbing dissolved polymers create a depletion attraction, we introduce an effective colloid–colloid potential closely related to the Asakura–Oosawa model, but that does not have any discontinuities. In simulations, this model straightforwardly allows the calculation of the pressure from the virial formula, and the phase transition in the bulk from the liquid to crystalline solid can be accurately located from a study where a stable coexistence of a crystalline slab with a surrounding liquid phase occurs. For this model, crystalline nuclei surrounded by fluid are studied both by identifying the crystal–fluid interface on the particle level (using suitable bond orientational order parameters to distinguish the phases) and by ‘thermodynamic’ means, i.e. the latter method amounts to compute the enhancement of chemical potential and pressure relative to their coexistence values. We show that the chemical potential can be obtained from simulating thick films, where one wall with a rather long-range repulsion is present, since near this wall, the Widom particle insertion method works, exploiting the fact that the chemical potential in the system is homogeneous. Finally, the surface excess free energy of the nucleus is obtained, for a wide range of nuclei volumes. From this method, it is established that classical nucleation theory works, showing that for the present model, the anisotropy of the interface excess free energy of crystals and their resulting non-spherical shape has only a very small effect on the barrier.

Solubility of Sodium in Sodium Chloride: A Density Functional Theory Molecular Dynamics Study

We present the first density functional theory (DFT) molecular dynamics (MD) study of metal-molten-salt solutions using the isobaric-isothermal (NpT) ensemble. We study solutions of excess sodium in molten sodium chloride of concentrations ranging from 0 to 11% Na. We use the Widom particle insertion method to calculate the chemical potential of sodium across this concentration range, and find the solubility limit to be between 3 to 6%, in excellent agreement with the experimental solubility limit. We use the same particle insertion method to calculate the standard reduction potential of sodium in molten sodium chloride as 3.9 ± 0.6 V, in good agreement with the experimental value of 3.18 V. We demonstrate the robustness of our DFT-MD particle insertion method, and discuss applications of our method to modeling electrowinning recycling processes of rare earth metals.

Monte Carlo tests of nucleation concepts in the lattice gas model

The conventional theory of homogeneous and heterogeneous nucleation in a supersaturated vapor is tested by Monte Carlo simulations of the lattice gas (Ising) model with nearest-neighbor attractive interactions on the simple cubic lattice. The theory considers the nucleation process as a slow (quasistatic) cluster (droplet) growth over a free energy barrier ΔF(*), constructed in terms of a balance of surface and bulk term of a critical droplet of radius R(*), implying that the rates of droplet growth and shrinking essentially balance each other for droplet radius R=R(*). For heterogeneous nucleation at surfaces, the barrier is reduced by a factor depending on the contact angle. Using the definition of physical clusters based on the Fortuin-Kasteleyn mapping, the time dependence of the cluster size distribution is studied for quenching experiments in the kinetic Ising model and the cluster size ℓ(*) where the cluster growth rate changes sign is estimated. These studies of nucleation kinetics are compared to studies where the relation between cluster size and supersaturation is estimated from equilibrium simulations of phase coexistence between droplet and vapor in the canonical ensemble. The chemical potential is estimated from a lattice version of the Widom particle insertion method. For large droplets it is shown that the physical clusters have a volume consistent with the estimates from the lever rule. Geometrical clusters (defined such that each site belonging to the cluster is occupied and has at least one occupied neighbor site) yield valid results only for temperatures less than 60% of the critical temperature, where the cluster shape is nonspherical. We show how the chemical potential can be used to numerically estimate ΔF(*) also for nonspherical cluster shapes.

Direct first-principles chemical potential calculations of liquids

We propose a scheme that drastically improves the efficiency of Widom's particle insertion method by efficiently sampling cavities while calculating the integrals providing the chemical potentials of a physical system. This idea enables us to calculate chemical potentials of liquids directly from first-principles without the help of any reference system, which is necessary in the commonly used thermodynamic integration method. As an example, we apply our scheme, combined with the density functional formalism, to the calculation of the chemical potential of liquid copper. The calculated chemical potential is further used to locate the melting temperature. The calculated results closely agree with experiments.

Field induced gradient simulations: A high throughput method for computing chemical potentials in multicomponent systems

We present a simulation method for direct computation of chemical potentials in multicomponent systems. The method involves application of a field to generate spatial gradients in the species number densities at equilibrium, from which the chemical potential of each species is theoretically estimated. A single simulation yields results over a range of thermodynamic states, as in high throughput experiments, and the method remains computationally efficient even at high number densities since it does not involve particle insertion at high densities. We illustrate the method by Monte Carlo simulations of binary hard sphere mixtures of particles with different sizes in a gravitational field. The results of the gradient Monte Carlo method are found to be in good agreement with chemical potentials computed using the classical Widom particle insertion method for spatially uniform systems.

Calculation of chemical potentials of chain molecules by the incremental gauge cell method

The gauge cell Monte Carlo method is extended to calculations of the incremental chemical potentials and free energies of linear chain molecules. The method was applied to chains of Lennard-Jones beads with stiff harmonic bonds up to 500 monomers in length. We show that the suggested method quantitatively reproduces the modified Widom particle insertion method of Kumar et al. [S. K. Kumar, I. Szleifer, and A. Z. Panagiotopoulos, Phys. Rev. Lett. 66(22), 2935 (1991)], and is by an order of magnitude more efficient for long chains in terms of the computational time required for the same accuracy of chemical potential calculations. The chain increment ansatz, which suggests that the incremental chemical potential is independent of the chain length, was tested at different temperatures. We confirmed that the ansatz holds only for coils above the θ temperature. Special attention is paid to the effects of the magnitude of adsorption potential and temperature on the behavior of single chains in confinements that are comparable in size with the free chain radius of gyration. At sufficiently low temperatures, the dependence of the incremental chemical potential on the chain length in wetting pores is superficially similar to a capillary condensation isotherm, reflecting monolayer formation following by pore volume filling, as the chain length increases. We find that the incremental gauge cell method is an accurate and efficient technique for calculations of the free energies of chain molecules in bulk systems and nanoconfinements alike. The suggested method may find practical applications, such as modeling polymer partitioning on porous substrates and dynamics of chain translocation into nanopores.

Solute effects on supercritical fluid

Supercritical carbon dioxide has been widely used as solvent, although its solubility is much lower than that of typical organic solvents in liquid state in general. Entrainer effects are known and utilised, in which adding low-weight molecules, such as alcohols (e.g. methanol and ethanol) and hydrocarbons (e.g. pentane and hexane), increases the solubility and the chemical reactivity. In this study, we investigated the change of local fluid structure and thermodynamic properties with molecular dynamics simulations, focusing on the entrainer effect. A monatomic fluid mixture model based on the Lennard-Jones model potential was adopted for simplicity; the parameters of which roughly correspond to those of carbon dioxide (solvent) and ethanol (solute). A series of simulations were performed in the vicinity of the gas–liquid critical point (CP). The radial distribution function and the coordination number reveal the enhancement of the solvent local density. The solvent chemical potential was evaluated with the Widom's particle insertion method. The chemical potential change by adding solutes is smallest in magnitude at the CP, qualitatively agreeing with recent experimental analysis and interpretation with a small-angle neutron scattering method.

Investigation on the Solubility of SO2 and CO2 in Imidazolium-Based Ionic Liquids Using NPT Monte Carlo Simulation

The solubility of sulfur dioxide (SO(2)) and carbon dioxide (CO(2)) at P = 1 bar in a series of imidazolium-based room-temperature ionic liquids (RTILs) is calculated by Monte Carlo simulation in NPT ensemble using the OPLS-UA force field and Widom particle insertion method. The studied ILs were 1-butyl-3-methylimidazolium ([bmim](+)) tetrafluoroborate ([BF(4)](-)), [bmim](+) hexafluorophosphate ([PF(6)](-)), [bmim](+) bromide ([Br](-)), [bmim](+) nitrate ([NO(3)](-)), [bmim](+) bis-(trifluoromethyl) sulfonylimide ([Tf2N](-)), and 1-ethyl-3-methylimidazolium tetrafluoroborate ([emim][BF(4)]). To validate the simulations, the liquid density of studied ILs and the solubility of CO(2) in [bmim][PF(6)] was compared with corresponding experimental and theoretical studies reported in the literature, and a good agreement was obtained. The results of SO(2) solubility demonstrate that the SO(2) gas has the highest solubility in [bmim][NO(3)] and [bmim][Br] ILs and the lowest solubility in [bmim][PF(6)]. To describe the solubility order of polar gases such as SO(2) and nonpolar gases like CO(2), we have simulated the SO(2)/IL and CO(2)/IL mixtures which made possible to investigate the interaction of solute molecules with anions and cations in the liquid phase. We introduced the ratio of solute-IL interaction over cation-anion interaction energy density as an index for solubility of gases in ILs. The results show that the proposed index can describe the solubility order of SO(2) as well as CO(2) and it might be used as an alternative to standard methods of infinite dilution Henry's constant calculations when the solubility order is desired.

Analyzing the components of the free-energy landscape in a calcium selective ion channel by Widom's particle insertion method

The selectivity filter of the L-type calcium channel works as a Ca(2+) binding site with a very large affinity for Ca(2+) versus Na(+). Ca(2+) replaces half of the Na(+) ions in the filter even when these ions are present in 1 μM and 30 mM concentrations in the bath, respectively. The energetics of this strong selectivity is analyzed in this paper. We use Widom's particle insertion method to compute the space-dependent profiles of excess chemical potential in our grand canonical Monte Carlo simulations. These profiles define the free-energy landscape for the various ions. Following Gillespie [Biophys. J. 94, 1169 (2008)], the difference of the excess chemical potentials for the two competing ions defines the advantage that one of the ions has over the other in the competition for space in the crowded selectivity filter. These advantages depend on ionic bath concentrations: the ion that is present in the bath in larger quantity (Na(+)) has the "number" advantage which is balanced by the free-energy advantage of the other ion (Ca(2+)). The excess chemical potentials are decomposed into hard sphere exclusion and electrostatic components. The electrostatic terms correspond to interactions with the mean electric field produced by ions and induced charges as well to ionic correlations beyond the mean field description. Dielectrics are needed to produce micromolar Ca(2+) versus Na(+) selectivity in the L-type channel. We study the behavior of these terms with changes in bath concentrations of ions, charges, and diameters of ions, as well as geometrical parameters such as radius of the pore and the dielectric constant of the protein. Ion selectivity in calcium binding proteins probably has a similar mechanism.

Pair interactions between complex mesoscopic particles from Widom's particle-insertion method

We demonstrate that Widom's particle insertion technique provides a convenient and efficient method to determine the effective pair interaction between complex, composite soft-matter particles in the zero-density limit. By means of three different test systems, i.e. amphiphilic dendrimers, electrostatic polymers and colloids coated with electrostatic polymers, we demonstrate the validity and the power of the presented method.

Entropy of the hard-core model

The entropy of the hard-core model is analyzed subject to dimension, system size and particle density. Inspired by the exactly solvable one-dimensional case and the results from a virial expansion, we suggest a scaling relation for the two-dimensional and three-dimensional case. These relations are compared with numerical results, obtained by a modified Widom's particle insertion method. The notion of finite size effects, relevant for atomic glass-forming systems, is discussed in the context of the present results.