Widom Particle Insertion Research Articles

PACMAN: A Robust Partial Atomic Charge Predicter for Nanoporous Materials Based on Crystal Graph Convolution Networks.

We report a fast and easy method (PACMAN) to assign partial atomic charges on metal-organic framework (MOF) and covalent-organic framework (COF) crystal structures based on graph convolution networks (GCNs) trained on >1.8 million high-fidelity partial atomic charge data obtained from the Quantum Metal-Organic Framework (QMOF) database. The developed model shows outstanding performance, achieving a mean absolute error (MAE) of 0.0055 e (test set performance) while maintaining consistency with DDEC6, Bader, and CM5 charges across diverse chemistry and topologies of MOFs and COFs. We find that the new method accurately assigns partial atomic charges for ion-containing nanoporous materials, which has not been possible in previous machine learning (ML) models. Grand canonical Monte Carlo (GCMC) simulation results for CO2 and N2 uptakes and the Widom particle insertion calculation for Henry's law constant of water results based on PACMAN and the original DDEC6 charges show excellent agreements compared to other ML models reported in the literature. The runtime analysis of the new method demonstrates that the partial atomic charges of MOF and COF structures with up to 500 atoms can be obtained in less than 10 s. An easy-to-use web interface has been developed to facilitate the adoption of the developed model.

Canonical titration simulations.

We present a Monte Carlo approach for performing titration simulations in the canonical ensemble. The standard constant pH (cpH) simulation methods are intrinsically grand canonical, allowing us to study the protonation state of molecules only as a function of pH in the reservoir. Due to the Donnan potential between a system and an (implicit) reservoir of a semi-grand canonical simulation, the pH of the reservoir can be significantly different from that of an isolated system, for an identical protonation state. The new titration method avoids this difficulty by using the canonical reactive Monte Carlo algorithm to calculate the protonation state of macromolecules as a function of the total number of protons present inside the simulation cell. The pH of an equilibrated system is then calculated using a new surface insertion Widom algorithm, which bypasses the difficulties associated with the bulk Widom particle insertion for intermediate and high pH values. To properly treat the long range Coulomb force, we use the Ewald summation method, showing the importance of the Bethe potential for calculating the pH of canonical systems.

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.

Quasi-Universal Solubility Behavior of Light Gases in Imidazolium-Based Ionic Liquids with Varying Anions: A Molecular Dynamics Simulation Study.

In this work, the temperature-dependent solvation behavior of a number of important light gases, such as carbon dioxide, xenon, krypton, argon, oxygen, methane, nitrogen, neon, and hydrogen, in two important imidazolium-based ionic liquids (ILs) of the type 1-n-alkyl-3-methylimidazolium hexafluorophosphate ([Cnmim][PF6]) and 1-n-alkyl-3-methylimidazolium tetrafluoroborate ([CnmimBF4]) with varying chain lengths (n = 2, 4, 6, and 8) are investigated using molecular dynamics simulations for a temperature range between 300 and 500 K at a pressure of 1 bar. The aim of this work is first to propose a reliable estimate for the temperature-dependent solubility behavior of (very) light gases, e.g., hydrogen and nitrogen, where reported experimental data are inconsistent. Moreover, we would like to rationalize the common features of the temperature-dependent solvation of light gases for various imidazolium-based ionic liquids. For the selected solute gases in our simulated imidazolium-based ILs, we applied the potential distribution theorem using both Bennet's overlapping distribution method (ODM) and Widom's particle insertion technique to determine the temperature-dependent solvation free energies with good statistical accuracy. We observed from the simulations that the quantity of the solvation free energy of a gas molecule and its temperature derivatives are connected in regard to each other at a chosen reference temperature. This trend was observed for all the studied light gases. Moreover, the computed solvation enthalpies of all gases obey an enthalpy-entropy compensation behavior, which is almost identical for all the studied ILs. Based on this observation, we report a correlation between the temperature-dependent solubility behavior of light gases in various ILs at their reference state so that we are now able to semiquantitively predict the temperature-dependent solubility behavior of a certain gas in various imidazolium-based ionic liquids based on a single solubility value of that gas in one of the ILs at a certain temperature.

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.

The chemical potential of a dipole in dipolar solvent at infinite dilution: Mean spherical approximation and Monte Carlo simulation

A new analytical expression was derived for the chemical potential of a hard sphere dipole in hard sphere dipole fluid at infinite dilution of the solute using the mean spherical approximation (MSA). A set of Monte Carlo (MC) simulations has been carried out to investigate the scope of applicability of the derived equation. The mean reaction field (MRF) approach was used in our MC computations. Two different MC methods (Widom particle insertion and thermodynamic integration) were applied for obtaining the chemical potential change associated with the dipole creation at the solute particle to provide adequate accuracy of the MC simulations. Also, corresponding changes in the mean potential energy were calculated by direct method and by thermodynamic integration. The solvation energies have been obtained for the systems of dipolar hard spheres with reduced dipole moment 1.0 at the reduced densities 0.2, 0.5, and 0.8. Computations have been made for solute particles with the reduced dipole moment varied from 0.0 to 1.5 and the hard sphere diameter varied from 0.5 to 2.0. The variation of those quantities with the molecular parameters was analyzed and compared with the MSA equation and Kirkwood classical expressions. It was found that the MSA calculations agree relatively well with MC simulations at densities less than 0.5 and solute dipole moment less than 1.0.

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.

Computational Study of ZIF-8 and ZIF-67 Performance for Separation of Gas Mixtures

ZIF-67, a modification of ZIF-8 framework through Zn substitution with Co, is tested for the first time for the separation of ethylene/ethane mixture using molecular simulations. The framework consists of cages connected with narrow apertures, which exhibit flexibility through a swelling motion, allowing for relatively large penetrants to diffuse. ZIF-67 demonstrates an enhanced separation for the specific mixture. Various computational techniques are employed (conventional molecular dynamics and Monte Carlo simulations, umbrella sampling, and Widom particle insertion), and the separation mechanism is investigated in terms of sorption and diffusion, for both ZIF-8 and ZIF-67. The stiffer bonding of Co with the adjacent N atoms results in a tighter structure and an aperture with smaller size and lower swelling amplitude than ZIF-8. The diffusion results show a clear dependency of the kinetic-driven separation on the aperture flexibility of the different frameworks. The diffusivities of different sized mole...

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.

The test-particle induced inhomogeneous direct correlation functions and extensions of Widom's theorem: Impacts on the incremental chemical potentials and high-order correlation functions

We develop the potential distributions of several test particles to obtain a hierarchy of the nonuniform singlet direct correlation functions (s-DCFs). These correlation functions are interpreted as the segmental chemical potentials or works of insertion of successive test particles in a classical fluid. The development has several interesting consequences: (i) it extends the Widom particle insertion formula to higher-order theorems, the first member gives the chemical potential as in the original theorem, the second member gives the incremental energy for dimer formation, with higher members giving the energies for forming trimers, tetramers, etc. (ii) The second and third order s-DCFs can be related to the cavity distribution functions y((2)) and y((3)) in the liquid-state theory. Thus we can express the triplet cavity function y((3)) in terms of these s-DCFs in an exact form. This enables us to calculate, as an illustration of the above theoretical developments, the numerical values of the s-DCFs via Monte Carlo (MC) simulation data on hard spheres. We use these data to critically analyze the commonly used approximations, the Kirkwood superposition (KSA) and the linear approximation (LA) for triplet correlation functions. An improved rule over KSA and LA is proposed for triplet hard spheres in the rolling-contact configurations. (iii) The s-DCFs are naturally suited for analyzing the chain-incremental Ansatz or hypothesis in the calculation of the chemical potentials of polymeric chain molecules. The first few segments of a polymer chain have been shown from extensive Monte Carlo simulations to not obey this Ansatz. By examining the insertion energies of successive segments through the s-DCFs, we are able to quantitatively decipher the decay of the segmental chemical potentials for at least the first three segments. Comparison with MC data on 4-mer and 8-mer hard-sphere fluids shows commensurate behavior with the s-DCFs. In addition, an analytical density functional theory is derived, through the potential distribution theorem, for obtaining these nonuniform direct correlation functions.

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.