Widom's Insertion Method Research Articles

Computing the solubility of argon and xenon in molten sodium chloride and potassium chloride salts

Molten salt reactors (MSRs) offer significant advancements in nuclear reactor safety and efficiency by operating at higher temperatures and lower pressures compared to traditional reactors. A critical aspect of MSR operation involves understanding the solubility of fission byproducts, particularly noble gases, in the molten salts used. This study employs molecular dynamics (MD) simulations to compute Henry’s law constants and enthalpies of solvation for argon and xenon in molten sodium chloride (NaCl) and potassium chloride (KCl). We developed a new pairwise potential for the noble gas and salt interactions based on first principles calculations. We then used this potential to calculate Henry’s law constants of the two gases in the molten salts, which were modeled using both a rigid ion model (RIM) and a polarizable ion model (PIM). The solubility calculations, performed using the Widom insertion method, show qualitative agreement with limited experimental data, highlighting the temperature dependence and greater solubility of both gases in KCl compared to NaCl. Additionally, free volume analysis elucidated the role of available space within the molten salts in governing solubility trends. Our findings suggest that PIM trajectories provide more reliable predictions for noble gas solubility than RIM due to their accurate density representation. These results enhance understanding of gas solubility in MSR environments, and the methods can be readily extended to other systems.

Modelling across Multiple Scales to Design Biopolymer Membranes for Sustainable Gas Separations: 1—Atomistic Approach

In this work, we assessed the CO2 and CH4 sorption and transport in copolymers of 3-hydroxybutyrate and 3-hydroxyvalerate (PHBV), which showed good CO2 capture potential in our previous papers, thanks to their good solubility-selectivity, and are potential biodegradable alternatives to standard membrane-separation materials. Experimental tests were carried out on a commercial material containing 8% of 3-hydroxyvalerate (HV), while molecular modelling was used to screen the performance of the copolymers across the entire composition range by simulating structures with 0%, 8%, 60%, and 100% HV, with the aim to provide a guide for the selection of the membrane material. The polymers were simulated using molecular dynamics (MD) models and validated against experimental density, solubility parameters, and X-ray diffraction. The CO2/CH4 solubility-selectivity predicted by the Widom insertion method is in good agreement with experimental data, while the diffusivity-selectivity obtained via mean square displacement is somewhat overestimated. Overall, simulations indicate promising behaviour for the homopolymer containing 100% of HV. In part 2 of this series of papers, we will investigate the same biomaterials using a macroscopic model for polymers and compare the accuracy and performance of the two approaches.

In silico screening of nanoporous materials for urea removal in hemodialysis applications.

The design of miniaturized hemodialysis devices, such as wearable artificial kidneys, requires regeneration of the dialysate stream to remove uremic toxins from water. Adsorption has the potential to capture such molecules, but conventional adsorbents have low urea/water selectivity. In this work, we performed a comprehensive computational study of 560 porous crystalline adsorbents comprising mainly covalent organic frameworks (COFs), as well as some siliceous zeolites, metal organic frameworks (MOFs) and graphitic materials. An initial screening using Widom insertion method assessed the excess chemical potential at infinite dilution for water and urea at 310 K, providing information on the strength and selectivity of urea adsorption. From such analysis it was observed that urea adsorption and urea/water selectivity increased strongly with fluorine content in COFs, while other compositional or structural parameters did not correlate with material performance. Two COFs, namely COF-F6 and Tf-DHzDPr were explored further through Molecular Dynamics simulations. The results agree with those of the Widom method and allow to identify the urea binding sites, the contribution of electrostatic and van der Waals interactions, and the position of preferential urea-urea and urea-framework interactions. This study paves the way for a well-informed experimental campaign and accelerates the development of novel sorbents for urea removal, ultimately advancing on the path to achieve wearable artificial kidneys.

Calculating adsorption isotherms using the two-phase thermodynamic method and molecular dynamics simulations

We describe the calculation of adsorption isotherms from molecular dynamics simulations based on the two-phase thermodynamic (2PT) model. The 2PT model developed for bulk fluid phases treats the gas-like components as hard spheres (HSs), which correctly recovers the limiting behaviors of unconfined fluids. We showed that this treatment, however, does not always lead to the correct zero-loading behavior in strongly confining systems. For methane adsorption into zeolite MFI, the HS reference state underestimates entropy by up to 20% at low loadings and leads to an order-of-magnitude increase in the adsorption onset pressure. To fix these issues, we propose the use of ideal adsorbed gas (IAG) as the gas reference model, the properties of which can be computed using the Widom insertion method on an empty adsorbent. We further describe three routes to compute adsorption isotherms from the Helmholtz free energy at different loadings. Comparing against established Monte Carlo (MC) methods, we found that the adsorption isotherms obtained using the IAG reference state agrees to within 40%, which corresponds to deviations of <5% in adsorption free energy. The isotherms calculated using the HS reference state underestimate the adsorption uptake at low to medium loadings in strongly confining systems, but its accuracy improves at higher loadings and as the pore size increases relative to the sorbate diameter. The methods described here provide an alternative approach for computing adsorption isotherms when MC simulations in an open ensemble are undesirable and enable a direct comparison of computed adsorption thermodynamics with experiments.

Widom insertion method in simulations with Ewald summation

We discuss the application of the Widom insertion method for calculation of the chemical potential of individual ions in computer simulations with Ewald summation. Two approaches are considered. In the first approach, an individual ion is inserted into a periodically replicated overall charge neutral system representing an electrolyte solution. In the second approach, an inserted ion is also periodically replicated, leading to the violation of the overall charge neutrality. This requires the introduction of an additional neutralizing background. We find that the second approach leads to a much better agreement with the results of grand canonical Monte Carlo simulation for the total chemical potential of a neutral ionic cluster.

The simulation of a green room-temperature ternary solution of water, methanol and 1-ethyl-3-methyl imidazolium chloride by all-atom Monte Carlo and DFT computational approaches

A highly dense ternary solution of water, methanol and 1-ethyl-3-methyl imidazolium chloride [emim][Cl] was investigated by a Monte Carlo (MC) program written in C language, Density Functional Theory (DFT) and Atoms-In-Molecules (AIM) theory. The all-atom approach in MC calculations of a ternary solution, the use of nonzero LJ parameters for all hydrogen atoms in simulation box and the creation of a new strategy for filling simulation box with many molecules were adopted in this research to yield a comprehensive, accurate and different work. The radial distribution function, the molar heat capacities in constant volume and pressure, the thermal expansivity, and the isothermal compressibility of the ternary system were calculated. The chemical potentials were calculated at different temperatures through Widom insertion method and the behaviors of both ideal and excess parts of chemical potentials were analyzed versus temperature. The value of electric permittivity of solution, the change of number of methanol and water molecules and the replacement of chloride anion with iodide in ionic liquid were also explored through separate simulations. The solvation structure and coordination number were investigated for various systems. The optimized sites for water and methanol around [emim]+ recognized by DFT and the HOMO and LUMO isosurfaces were discussed. Also, the hydrogen bonds between components of ternary solution were determined by AIM theory. It was shown that the results of quantum mechanics calculations confirm those of MC simulations.

Computational Prediction of Water Sorption in Facilitated Transport Membranes

Polymeric facilitated transport membranes (FTMs) have emerged as an innovative class of promising carbon capture technology. Studies have shown that the transport properties of an FTM are significantly influenced by its water uptake. In order to better quantify FTM performance, we herein explore the potential of computational techniques to predict the equilibrium water uptake values of FTMs. Two prediction approaches were examined. First, the water sorption was explicitly simulated by iteratively conducting grand canonical Monte Carlo and molecular dynamics simulations. Second, the water sorption was predicted based on the chemical potential of the adsorbed water, which was calculated using the Widom insertion or continuous fractional component Monte Carlo (CFCMC) method. The chemical potential-based approach with CFCMC demonstrated good prediction of the equilibrium water uptake values of FTMs with poly(N-vinylformamide-co-vinylamine) as the fixed-site carrier and 2-(1-piperazinyl)ethylamine sarcosinate as the mobile carrier. The predicted water uptake values increased with increasing mobile carrier content and were in good agreement with the experimental values. The higher water uptake promoted the diffusion of CO2, N2, and mobile carrier, as well as slightly stifled the sorption of N2. Such an approach significantly contributes to a more comprehensive theoretical evaluation of FTMs.

Ab Initio Evaluation of Henry Coefficients Using Importance Sampling.

We present a new algorithm that allows for an efficient evaluation of the Henry coefficient of a guest molecule inside a porous material, which permits to use ab initio energy calculations. The Widom insertion method, which is currently used to compute these Henry coefficients, typically requires millions of energy evaluations. Our new methodology reduces this number by more than 1 order of magnitude, enabling the use of an ab initio potential energy surface. The methodology we propose is reminiscent of the well-known importance sampling technique which is frequently used in Monte Carlo integrations. First, a conventional Widom insertion simulation is performed using a force field. In the second step, the Widom results are used to select a limited number of configurations and only for these configurations the ab initio evaluation of the energy is required. Finally, by appropriately reweighting the latter energies, an accurate estimation of the ab initio Henry coefficient is possible at a moderate computational cost. We apply our methodology to the adsorption of CO2 in Mg-MOF-74, a prototypical case where interactions of a polar guest molecule with unsaturated metal sites dominate the adsorption mechanism. In this case generic force fields such as UFF or Dreiding are inappropriate and the use of ab initio methods is indispensable. In a second case study, we compute Henry coefficients of methane in UiO-66 using different levels of theory. We pay particular attention to the influence of the dispersion corrections and the role of many-body effects. For the final example, we qualitatively investigate adsorption features for a series of functionalized UiO-66 frameworks. Overall the cases we present show that accurate computations of Henry coefficients is extremely challenging, as different levels of theory provide strongly varying results. At the same time ab initio calculations have added value compared to force fields, as they provide a physically more sound description of the adsorption mechanism and in some cases clearly improve correspondence with experiment.

Chemical potential calculations in non-homogeneous liquids.

The numerical computation of chemical potential in dense non-homogeneous fluids is a key problem in the study of confined fluid thermodynamics. To this day, several methods have been proposed; however, there is still need for a robust technique, capable of obtaining accurate estimates at large average densities. A widely established technique is the Widom insertion method, which computes the chemical potential by sampling the energy of insertion of a test particle. Non-homogeneity is accounted for by assigning a density dependent weight to the insertion points. However, in dense systems, the poor sampling of the insertion energy is a source of inefficiency, hampering a reliable convergence. We have recently presented a new technique for the chemical potential calculation in homogeneous fluids. This novel method enhances the sampling of the insertion energy via well-tempered metadynamics, reaching accurate estimates at very large densities. In this paper, we extend the technique to the case of non-homogeneous fluids. The method is successfully tested on a confined Lennard-Jones fluid. In particular, we show that, thanks to the improved sampling, our technique does not suffer from a systematic error that affects the classic Widom method for non-homogeneous fluids, providing a precise and accurate result.

Molecular-level insight of the differences in the diffusion and solubility of penetrants in polypropylene, poly(propylmethylsiloxane) and poly(4-methyl-2-pentyne)

In order to design optimal polymer membranes it is important to have a molecular level understanding of penetrant transport in such membranes. As such, molecular dynamics simulation was employed to study the diffusion and solubility of gas molecules in polypropylene (PP), poly(propylmethylsiloxane) (PPMS) and poly (4-methyl-2-pentyne) (PMP). After the structures of PP, PPMS and PMP were established and relaxed, the average amplitudes of the oscillation of the main and branch chains of PP, PPMS and PMP were evaluated with a proposed analysis method; in addition the cavity size distributions of PP, PPMS and PMP were determined. PPMS has the largest average cavity size and accessible cavity fraction for penetrants, followed by PMP and PP, which has the smallest. The logarithm plot of mean-squared displacements (MSDs) versus time for the transport of methane in the three polymers revealed three regimes, namely the ballistic, the subdiffusive and the Fickian diffusive regimes. The ballistic regime of PMP is longer than that of PP and shorter than that of PPMS, because the average cavity size and accessible cavity fraction of PMP are larger than that of PP and smaller than that of PPMS and the penetrants have less space in PP and more space in PPMS to move freely than they do in PMP before they hit any matrix units. Statistical tests showed that the random walk on a fractal (RWF) model was the most appropriate model to explain the subdiffusion phenomena. The subdiffusive regime is induced by the trap of penetrants in cavities before they get the chances to jump from cavities to cavities nearby. Next the diffusivities of methane and n-butane were estimated from the plot of MSDs versus time for the penetrant transport in the Fickian diffusive regime. The diffusivities of penetrants are much smaller in PP than they are in PPMS and PMP and these values are larger in PPMS than they are in PMP. With both unbiased and biased Widom insertion methods, the solubility coefficients K of gas molecules in the polymers were calculated, though biased Widom insertion method has significantly faster calculation speed. PPMS has the largest solubilities of methane and n-butane, followed by PMP and PP, which has the smallest. Though the permeabilities of both methane and n-butane in PPMS are larger than the corresponding values in PMP, the selectivity of n-butane over methane in PPMS is lower than that in PMP. The research results may be employed in membrane design.

Determination of the Residual Entropy of Simple Mixtures by Monte Carlo Simulation.

Extensive Monte Carlo simulations were performed for (neon + krypton) mixtures for temperatures between 200 K and 600 K and pressures up to 1 GPa, using Lennard-Jones potentials to describe the intermolecular interactions. The residual entropies were obtained via Widom's insertion method, as well as via an integration technique. At high pressures, the residual entropy is, to a very good approximation, a linear function of λa-1, which is the reciprocal value of the average Monte Carlo displacement parameter that gives the acceptance ratio a for translational moves. The slope of this linear function varies linearly with the mole fraction and is related to the effective collision diameters of the molecules. As the displacement parameter is available during Monte Carlo simulations of fluids, its linearity with the residual entropy can be used to compute this properties with negligible computational effort at high densities, when particle insertion methods become unreliable.

On the consistency of NVT, NPT, μVT and Gibbs ensembles in the framework of kinetic Monte Carlo – Fluid phase equilibria and adsorption of pure component systems

This paper aims to show the consistency between simulations of fluid phase properties, obtained with various ensembles, developed within the framework of kinetic Monte Carlo (kMC) simulation: NVT (canonical), NPT (isothermal-isobaric systems), μVT (grand canonical) and Gibbs ensembles, to ensure the reliability of the kMC methodology. The advantages of the kMC scheme, as compared to the conventional Metropolis Monte Carlo, are: (1) accurate determination of the chemical potentials compared to the Widom insertion method, and (2) a rejection-free algorithm, making the implementation of the kMC scheme simpler. For internal consistency in a grand canonical ensemble simulation, we have developed a means to calculate the intrinsic chemical potential of the system accurately, which must be the same (within statistical error of the simulation) as the specified chemical potential to ensure convergence to equilibrium. We test the consistency of canonical (NVT-kMC) and grand canonical (GC-kMC) ensembles for argon adsorption at 87K and 120K in a uniform open-ended slit pore, and hence derive governing factors affecting hysteresis in the isotherm and the microscopic mechanisms of condensation and evaporation.

Probing gas adsorption in MOFs using an efficient ab initio widom insertion Monte Carlo method.

We propose a novel biased Widom insertion method that can efficiently compute the Henry coefficient, KH , of gas molecules inside porous materials exhibiting strong adsorption sites by employing purely DFT calculations. This is achieved by partitioning the simulation volume into strongly and weakly adsorbing regions and selectively biasing the Widom insertion moves into the former region. We show that only few thousands of single point energy calculations are necessary to achieve accurate statistics compared to many hundreds of thousands or millions of such calculations in conventional random insertions. The methodology is used to compute the Henry coefficient for CO2 , N2 , CH4 , and C2 H2 in M-MOF-74(M = Zn and Mg), yielding good agreement with published experimental data. Our results demonstrate that the DFT binding energy and the heat of adsorption are not accurate enough indicators to rank the guest adsorption properties at the Henry regime. © 2016 Wiley Periodicals, Inc.

Chemical potential calculations in dense liquids using metadynamics

The calculation of chemical potential has traditionally been a challenge in atomistic simulations. One of the most used approaches is Widom's insertion method in which the chemical potential is calculated by periodically attempting to insert an extra particle in the system. In dense systems this method fails since the insertion probability is very low. In this paper we show that in a homogeneous fluid the insertion probability can be increased using metadynamics. We test our method on a supercooled high density binary Lennard-Jones fluid. We find that we can obtain efficiently converged results even when Widom's method fails.

MCCCS Towhee: a tool for Monte Carlo molecular simulation

The history of the Monte Carlo for complex chemical systems Towhee open source Monte Carlo molecular simulation tool is discussed. A proof is given that the Widom insertion method for computing the chemical potential is formally correct even when combined with the most general version of arbitrary trial distribution configurational-bias Monte Carlo. A simulation strategy for computing single component vapour–liquid phase coexistence curves is presented as a guide for inexperienced practitioners of Monte Carlo simulations. A review of papers that cite the Towhee code is presented. The paper concludes with a discussion about releasing and sustaining a simulation package that uses an open source software license.

Solubility of solids in supercritical fluid using the hard-body expanded virial equation of state

The solubility of solids in supercritical fluids was examined using the hard-body expanded virial equation of state (HBE-VEOS). Due to the semi-soft core potential function, the developed EOS is analytic and has a closed form, which makes it appropriate for the engineering level. Molecular shape was considered in that the molecules were standardized as hard convex bodies (HCB) such as ellipsoids or spherocylinders. Distance-angle decoupling approximation was used in the integration of the cluster integrals to derive a simple equation for the second virial coefficient. Using the semi-soft core potential function, isothermal-isobaric Monte Carlo simulations were performed to generate a pVT diagram of pure carbon dioxide and canonical Monte Carlo simulations using the Widom insertion method were performed to calculate the residual chemical potentials for infinite dilute systems. Experimental solubility data of naphthalene, benzoic acid and phenanthrene in carbon dioxide were compared with the calculation of the developed EOS mapping of the molecules into ellipsoids. We confirmed from the results that second order expansion in the HBE-VEOS is the optimal combination in terms of both simplicity and accuracy for industrial supercritical extraction of solids.

An NPT Monte Carlo Molecular Simulation-Based Approach to Investigate Solid-Vapor Equilibrium: Application to Elemental Sulfur-H2S System

In this work, a method to estimate solid elemental sulfur solubility in pure and gas mixtures using Monte Carlo (MC) molecular simulation is proposed. This method is based on Isobaric-Isothermal (NPT) ensemble and the Widom insertion technique for the gas phase and a continuum model for the solid phase. This method avoids the difficulty of having to deal with high rejection rates that are usually encountered when simulating using Gibbs ensemble. The application of this method is tested with a system made of pure hydrogen sulfide gas (H2S) and solid elemental sulfur. However, this technique may be used for other solid-vapor systems provided the fugacity of the solid phase is known (e.g., through experimental work). Given solid fugacity at the desired pressure and temperature, the mole fraction of the solid dissolved in gas that would be in chemical equilibrium with the solid phase might be obtained. In other words a set of MC molecular simulation experiments is conducted on a single box given the pressure and temperature and for different mole fractions of the solute. The fugacity of the gas mixture is determined using the Widom insertion method and is compared with that predetermined for the solid phase until one finds the mole fraction which achieves the required fugacity. In this work, several examples of MC have been conducted and compared with experimental data. The Lennard-Jones parameters related to the sulfur molecule model (ɛ, σ) have been optimized to achieve better match with the experimental work.

A thermodynamic study of the mid-density scheme to determine the equilibrium phase transition in cylindrical pores

We present a thermodynamic analysis of the mid-density scheme that was introduced recently as a simple alternative to more complex methods (such as the gauge cell and the thermodynamic integration procedures) for determining the equilibrium phase transition in pores. A construction of configurations having densities falling between those of the low- and high-density states inside the hysteresis loop is carried out with a canonical Monte Carlo simulation (NVT, where the number of particles N, the system volume V and the temperature T in the system are constant), and it is found that this ‘nucleation’ state has two coexisting phases of alternating liquid bridges and gas cavities. This state is maintained at a threshold chemical potential, μeq (obtained with the Widom insertion method), which is the coexistence chemical potential. When this ‘nucleation’ state is in contact with an infinite reservoir (grand canonical ensemble) of chemical potential less than μeq, the system will evolve to a low-density state. On the other hand, when it is exposed to a chemical potential greater than μeq, the system will evolve to a high-density state. We test this method with a number of examples of cylindrical pores, and investigate the effects of the pore size, the simulation box length for the case of infinite pores, the pore length for the case of finite pores and the temperature. When the pore is finite in length, the coexistence is no longer a first-order transition, like the one observed in an infinitely long pore, but instead shows a second-order transition with a continuous, but sharp change from a low-density state to a high-density state.

Atomistic Simulations of Macromolecular Crowding

Protein folding, binding, and aggregation are known to be affected by macromolecular crowding [1], which is an integral part of intracellular environments. Realistic modeling of intracellular crowding requires computer simulations. In direct simulations of test proteins mixed with crowders, it has only been practical to represent the proteins at a coarse-grained level. Our recently developed “postprocessing” approach has made it possible to represent test proteins at the atomic level [2, 3]. In this approach, the motions of a test protein and those of the crowders are followed in two separate simulations. The effects of crowding are then modeled by calculating Δμ, the crowding-induced change in the chemical potential of the test protein. For a repulsive type of protein-crowder interactions, Δμ is related to the fraction, f, of allowed placements of the test protein into a box of crowders. For spherical crowders, two methods have been devised to calculate f. The first is an efficient implementation of Widom's insertion method. The second is a theoretical prediction, which uses the volume, surface area, and linear size defined on a “crowder-exclusion” surface. These methods have been applied to study crowding effects on protein folding, binding, and internal dynamics. We have also tested the postprocessing approach against direct simulations of folding-unfolding and open-to-closed transitions in the presence of crowders. To calculate f for atomistic crowders, we have just devised an algorithm based on fast Fourier transform. These computational tools establish a solid foundation for realistic modeling of intracellular environments.[1] H.-X. Zhou, G. Rivas, and A. P. Minton, Annu Rev Biophys 37, 375 (2008).[2] S. Qin and H.-X. Zhou, Biophys J 97, 12 (2009).[3] J. Batra, K. Xu, S. Qin, and H.-X. Zhou, Biophys J 97, 906 (2009).

Barrier properties of small gas molecules in amorphous cis-1,4-polybutadiene estimated by simulation

The solubility and diffusivity of small gas molecules in amorphous cis-1,4-polybutadiene were estimated in the temperature range 250–400 K, using both molecular dynamics simulations and the transition state theory implementation of the Widom insertion method. A comparison of the methods is given and the results obtained are compared with available experimental data. The accuracy in predicting diffusion and solubility coefficients of a range of small gas molecules in an amorphous polymer of both methods is good when compared to experimental values. The effects of the temperature and the models size are also examined. Selectivity to oxygen and nitrogen are estimated for the various models studied as well. This work shows the potential of computational methods for the prediction of physical properties of industrial importance like selectivity.