Canonical Ensemble Simulations Research Articles

Generalized Green's function molecular dynamics for canonical ensemble simulations.

The need of small integration time steps (∼1fs) in conventional molecular dynamics simulations is an important issue that inhibits the study of physical, chemical, and biological systems in real timescales. Additionally, to simulate those systems in contact with a thermal bath, thermostating techniques are usually applied. In this work, we generalize the Green's function molecular dynamics technique to allow simulations within the canonical ensemble. By applying this technique to one-dimensional systems, we were able to correctly describe important thermodynamic properties such as the temperature fluctuations, the temperature distribution, and the velocity autocorrelation function. We show that the proposed technique also allows the use of time steps one order of magnitude larger than those typically used in conventional molecular dynamics simulations. We expect that this technique can be used in long-timescale molecular dynamics simulations.

Understanding pore-level phenomena of n-alkanes at high pressures

Adsorption of shale related n-alkane fluids in graphite nano-channels is studied. Pure component canonical (NVT) ensemble simulations are performed over a range of isotherms and pore widths. Results for n-hexane are used to validate our approach by comparing them with experimental data in the open literature. The impact of pore throat diameter on adsorption and internal energies of departure is also reported and further demonstrate that there is a carbon chain length dependence on adsorption. All pure component results reported in this work are fit to a Langmuir adsorption model, which provides easy re-use in future studies. Numerical adsorption results for pure component methane adsorption are compared to experimental reservoir conditions. Finally, Ideal Adsorbed Solution Theory (IAST) is applied to the resulting Langmuir isotherms in order to provide some insight into the phase equilibria for applicable shale gas mixtures.

Effect of pore geometry on the compressibility of a confined simple fluid.

Fluids confined in nanopores exhibit properties different from the properties of the same fluids in bulk; among these properties is the isothermal compressibility or elastic modulus. The modulus of a fluid in nanopores can be extracted from ultrasonic experiments or calculated from molecular simulations. Using Monte Carlo simulations in the grand canonical ensemble, we calculated the modulus for liquid argon at its normal boiling point (87.3 K) adsorbed in model silica pores of two different morphologies and various sizes. For spherical pores, for all the pore sizes (diameters) exceeding 2 nm, we obtained a logarithmic dependence of fluid modulus on the vapor pressure. Calculation of the modulus at saturation showed that the modulus of the fluid in spherical pores is a linear function of the reciprocal pore size. The calculation of the modulus of the fluid in cylindrical pores appeared too scattered to make quantitative conclusions. We performed additional simulations at higher temperature (119.6 K), at which Monte Carlo insertions and removals become more efficient. The results of the simulations at higher temperature confirmed both regularities for cylindrical pores and showed quantitative difference between the fluid moduli in pores of different geometries. Both of the observed regularities for the modulus stem from the Tait-Murnaghan equation applied to the confined fluid. Our results, along with the development of the effective medium theories for nanoporous media, set the groundwork for analysis of the experimentally measured elastic properties of fluid-saturated nanoporous materials.

Transferability of cross-interaction pair potentials: Vapor-liquid phase equilibria of n-alkane/nitrogen mixtures using the TAMie force field

The van der Waals contribution of force fields designed for calculating phase equilibria and thermodynamic properties are usually adjusted to experimental data of pure components. Using these force fields for mixtures by applying combining rules for the cross-wise interaction potentials does not always lead to satisfactory results. This study considers the cross-wise van der Waals energy parameters εαβ among individual pairs of (united-)atom groups as adjustable (but subsequently transferable) parameters. The cross-energy parameters are simultaneously adjusted to an experimental training set of several mixtures. An analytic equation of state, the Perturbed-Chain Polar Statistical Associating Fluid Theory (PCP-SAFT), is used for ensuring very swift convergence of the optimization procedure. As an example, we consider the phase equilibria of n-alkane/nitrogen mixtures. The cross-energy parameters (-CH3 to N and -CH2- to N) are adjusted to three mixtures as the training set. We explore the transferability of the so-obtained force field for varying temperatures and for mixtures that were not considered in the training set. Phase equilibria of nine n-alkane/nitrogen mixtures (each at various temperatures) are determined using Monte Carlo simulations in the grand canonical ensemble. Results for vapor-liquid equilibrium calculations show very good agreement with experimental data.

Heat capacities of supercritical fluids via Grand Canonical ensemble simulations

In this work, suitable mathematical relationships to compute isobaric heat capacities from molecular simulations in the Grand Canonical (GC) ensemble are derived and tested via Monte Carlo methods. Using atomistic classical force fields, the residual isobaric heat capacities of pure carbon dioxide (CO2) and pure methanol (MeOH) were obtained at supercritical conditions (with critical properties estimated from a finite-size scaling analysis). The total isobaric heat capacity was determined by combining the residual isobaric heat capacity obtained from molecular simulations with the ideal gas contributions obtained from experimental correlations. Isobaric heat capacities generated from both GC and Isothermal–Isobaric ensemble simulations were compared to predictions from accurate equations of state (EOS)s for CO2 and MeOH at corresponding reduced temperatures and pressures. Isobaric heat capacities calculated from both ensembles were in good agreement with those obtained from the Span and Wagner EOS for CO2 and the IUPAC EOS for MeOH. For comparable computation times, simulations run in the GC ensemble generate results with significantly lower statistical uncertainty than those run in the Isothermal–Isobaric ensemble.

Neural network potentials for dynamics and thermodynamics of gold nanoparticles.

For understanding the dynamical and thermodynamical properties of metal nanoparticles, one has to go beyond static and structural predictions of a nanoparticle. Accurate description of dynamical properties may be computationally intensive depending on the size of nanoparticle. Herein, we demonstrate the use of atomistic neural network potentials, obtained by fitting quantum mechanical data, for extensive molecular dynamics simulations of gold nanoparticles. The fitted potential was tested by performing global optimizations of size selected gold nanoparticles (Aun, 17 ≤ n ≤ 58). We performed molecular dynamics simulations in canonical (NVT) and microcanonical (NVE) ensembles on Au17, Au34, Au58 for a total simulation time of around 3 ns for each nanoparticle. Our study based on both NVT and NVE ensembles indicate that there is a dynamical coexistence of solid-like and liquid-like phases near melting transition. We estimate the probability at finite temperatures for set of isomers lying below 0.5 eV from the global minimum structure. In the case of Au17 and Au58, the properties can be estimated using global minimum structure at room temperature, while for Au34, global minimum structure is not a dominant structure even at low temperatures.

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.

Hydrogen adsorption on graphane: An estimate using ab-initio interaction

The interaction between hydrogen molecule and graphane, material synthesized when a graphene plane is fully functionalized by hydrogen atoms, is assessable by quantum mechanical ab-initio calculations. Therefore for hydrogen, it is possible to estimate the adsorption properties of a porous material similar to activated carbons, the adsorbent surface of which is made of graphane planes instead of graphene or basal graphitic planes. The calculation realized by Monte-Carlo simulations in the grand canonical ensemble shows that the hydrogen adsorption of graphane stays qualitatively similar to that of graphene.

Magic compositions in Pd-Au nanoalloys

We identified new magic compositions of Pd-Au nanolloys of simple symmetries according to the icosahedral and the cubic truncated octahedral structures. Chemically ordered structures as well as segregated ones are characterized using Monte Carlo simulations in semi-grand canonical ensemble and tight-binding semi-empirical potential fitted to ab initio calculations for the Pd-Au system. Ordering and surface (or core) segregation phenomena can be distinguished by their behavior as a function of the temperature: chemical ordering at surface or in the core disappears above a critical temperature (which is lower or equal to the bulk alloy order-disorder temperature), whereas surface/core segregation (core-shell structures) remains at high temperature, although the interfacial profile is less abrupt.

Use of molecular dynamics simulations to estimate the solubility of menadione in supercritical CO2 using Chrastil's model

The binary system of menadione in explicit supercritical carbon dioxide (SC-CO2) was studied using molecular dynamics (MD) simulations, with the objective to understand the nature of interactions between menadione and SC-CO2 at different temperatures and pressures in order to complement experimental solubility measurements. A force field was developed for menadione and tested by comparing computed and experimental monoclinic crystal structures at 283 K and 0.01 MPa. Lattice parameters obtained from anisotropic isothermal-isobaric MD simulations agreed reasonably well with experimental values, with an average absolute relative deviation (AARD%) less than 7%. A previously validated force field for SC-CO2 was used, and simple mixing rules were used to describe cross interactions. Canonical ensemble MD simulations were used to estimate the association number for CO2 about menadione and the enthalpy required to form a SC-CO2 solvate complex with menadione as a function of temperature and CO2 density. Spatial distribution functions were computed to better understand the nature of the molecular-level interactions between menadione and SC-CO2 as well as between associating menadione molecules. This work is the first part of a study that uses MD simulations as the main tool to represent a binary system. The MD methodologies of this work will be applied to represent our future studies of menadione derivatives in SC-CO2.

Phase behaviour of a continuous shouldered well model fluid. A grand canonical Monte Carlo study.

The phase behavior of the continuous shouldered well model fluid proposed by Franzese [J. Mol. Liq. 136 (2007) 267] was examined using the Monte Carlo computer simulations in the grand canonical ensemble. The essential parts of the vapour-liquid and liquid-liquid coexistence envelopes were obtained. The Widom lines departing from coexistence envelopes were calculated using maxima of the fluctuations of the number of particles as a function of chemical potential along various isotherms. The region embracing anomalies in the properties of the model was located using the approximate criterion that involves the excess pair entropy.. The temperature of maximum density line was built by performing canonical Monte Carlo simulations. Our results are consistent with previous results from molecular dynamics constant pressure-constant temperature simulations and provide wider insight into the phase behavior of the model by using the chemical potential as the external parameter.

Monte Carlo Simulation of Near- and Supercritical Hexane Fluid and Physisorption Phase Behavior

Monte Carlo (MC) simulations were used to determine the near- and supercritical fluid phase behavior of hexane and the adsorption thermodynamics of hexane on cobalt. Canonical ensemble MC simulations were used to compute the equation of state and chemical potentials, and Gibbs ensemble MC simulations (GEMC) were used to compute fluid phase coexistence. Overall, MC results with the TraPPE-UA force field were successful in representing near- and supercritical hexane thermodynamic behavior, and represent an improvement over treatments with a classical cubic equation of state. In addition, GEMC simulations predict a critical point for TraPPE-UA hexane of 509 K, 35 bar, and 0.22 g/mL, in excellent agreement with critical constants for real hexane. With regard to physisorption, hexane physisorption onto a (0001) cobalt surface at 523.15 K was modeled using grand canonical Monte Carlo (GCMC). GCMC excess adsorption results show crossover from adsorption to depletion around a bulk density of 0.1 g/mL and a global...

Prediction of Radon-222 Phase Behavior by Monte Carlo Simulation

Histogram-reweighting Monte Carlo simulations in the grand canonical ensemble are used to determine saturated liquid and vapor densities, vapor pressures, heats of vaporization, and compressibility factors for radon-222 from the normal boiling point to the critical point. An optimized intermolecular potential is developed by fitting parameters to reproduce experimental vapor pressures and the critical temperature. Vapor pressures are reproduced by simulation to within 2.2% of experiment, while the critical temperature and normal boiling point are reproduced exactly within the combined uncertainty of simulation and experimental data. The predictions of simulation are used to evaluate the reliability of reported experimental liquid densities, and those predicted by equations of state. While the predictions of simulation are in agreement with equations of state for the liquid density, significant differences are observed between the predictions of simulation and experiment. These results suggest that previously reported experimental liquid densities for radon-222 may be erroneous.

Prediction of fluid phase equilibria and interfacial tension of triangle-well fluids using transition matrix Monte Carlo

The triangle-well (TW) potential is a simple model which is able to capture the essence of the intermolecular attraction in real molecules. Transition matrix Monte Carlo simulations in the grand canonical ensemble (GC-TMMC) are performed to investigate the role of the range of attraction on the features of fluid phase equilibria. As the TW potential range increases, the vapour–liquid coexistence curves shift towards a higher temperature range with the critical temperature and pressure increasing, and the critical density values decreasing. These GC-TMMC results are in excellent agreement with the predictions of Gibbs ensemble Monte Carlo and replica exchange Monte Carlo (REMC) simulations reported in literature. Using the GC-TMMC method, the vapour pressures are also computed directly from the particle number probability distributions (PNPDs). It has been noted in literature that the surface tension values are computationally more expensive and difficult to determine than other coexistence properties using molecular simulations. The PNPDs from GC-TMMC simulations along with Binder’s formalism allow for the calculation of the interfacial tension with relative ease. Also, our simulation generated results for the interfacial tension are in good agreement with the literature data obtained using REMC (via the virial route) and the plots of our interfacial tension values as a function of temperature are smooth unlike the literature data.

Predicting the adsorption of n-perfluorohexane (n-C6F14) on BAM-P109 activated carbon using an ab initio force field

The adsorption isotherm of n-perfluorohexane on BAM-P109 activated carbon was predicted using an ab initio force field and a slit-like pore model representing the activated carbon. The force field parameters were derived from data calculated using quantum mechanics density functional theory with dispersion energy correction, and the total adsorption amount was calculated as the sum of the adsorption amounts of microscopic and mesoscopic pores. It was found that the Grand Canonical Monte Carlo method can efficiently calculate adsorption on microscopic pores but not on mesoscopic pores. Thus, we divided the mesoscopic pores into surface and free-volume regions and calculated the density and thickness of the adsorbed layers and the density of the vapor phase using canonical ensemble simulations. Using this protocol, we can calculate the adsorption isotherms in good agreement with the experimental data.

Study of critical dynamics in fluids via molecular dynamics in canonical ensemble.

With the objective of understanding the usefulness of thermostats in the study of dynamic critical phenomena in fluids, we present results for transport properties in a binary Lennard-Jones fluid that exhibits liquid-liquid phase transition. Various collective transport properties, calculated from the molecular dynamics (MD) simulations in canonical ensemble, with different thermostats, are compared with those obtained from MD simulations in microcanonical ensemble. It is observed that the Nosé-Hoover and dissipative particle dynamics thermostats are useful for the calculations of mutual diffusivity and shear viscosity. The Nosé-Hoover thermostat, however, as opposed to the latter, appears inadequate for the study of bulk viscosity.

Micro-structural Change During Nucleation: From Nucleus To Bicontinuous Morphology.

Although the microstructure of coexistence phase provides direct insights of the nucleation mechanism and their change is substantial in the phase transition, their study is limited due to the lack of suitable tools capturing the thermodynamically unstable transient states. We resolve this problem in computational study by introducing a generalized canonical ensemble simulation and investigate the morphological change of the nucleus during the water evaporation and condensation. We find that at very low pressure, where the transition is first order, classical nucleation theory holds approximately. A main nucleus is formed in the supersaturation near spinodal, and the overall shape of the nucleus is finite and compact. On increasing the pressure of the system, more nuclei are formed even before spinodal. They merge into a larger nuclei with a smaller free energy penalty to form ramified shapes. We suggest order parameters to describe the extent of fluctuation, and their relation to the free energy profile.

On the microscopic origin of the hysteresis loop in closed end pores – Adsorbate restructuring

We present here compelling evidence for adsorbate restructuring as the microscopic origin of hysteresis in closed end pores. Our argument is based on molecular simulations in the grand canonical (GCE) and meso-canonical (MCE) ensembles, for argon adsorbed at 87K in slit mesopores of different topologies (open end, closed end and closed pores). The MCE isotherms have sigmoidal van der Waals type loops, while the GCE isotherms exhibit hysteresis loops, which are confined between the gas-like and liquid-like spinodal points. One interesting feature, not previously recognized, is that the condensation in the MCE isotherms for pores of different topologies is identical, and it occurs at the same chemical potential, which is the coexistence chemical potential of two phases in equilibrium: the liquid-like condensate and the gas-like phase. At the onset of the condensation, the liquid-like condensate is a thin liquid bridge formed by taking molecules from the metastable adsorbed film, resulting in a thinner and stable adsorbed film. As more molecules are added into the pore, the liquid bridge thickens along the axial direction and both menisci advance toward the pore mouth. As the menisci approach the pore mouth, the liquid condensate is progressively restructured to become more cohesive. When the pressure is reduced in the grand canonical ensemble (open system), a lower pressure is required to disrupt this cohesive structure, resulting in hysteresis in the closed end pores. The classical theories, which conclude that isotherms for closed end pores must be reversible, are in error because they assume that the interfacial energy parameter (the product of the surface tension and the liquid molar volume) has the same value as in a bulk liquid.

An efficientmethod to determine chemical potential of mixtures in the isothermal and isobaric bulk phase with kineticMonte Carlo simulation

ABSTRACTA kinetic Monte Carlo (kMC) scheme in the NPT ensemble (constant number of molecules, pressure and temperature) has been developed to determine accurate chemical potentials for all components in a homogeneous mixture. The simulation requires two moves: (1) a displacement move and (2) a volume change move. In the former, the mobility rate of a selected molecule is determined by its interaction with all the other molecules in the system and is moved to a random position within the simulation box, according to the Rosenbluth algorithm, without any rejections (entropic sampling). The volume change move is decided by a comparison between either the instant pressure or the partial average pressure (with long-range correction) and the specified pressure and is carried out much less frequently than the displacement move. We applied this NPT scheme to a number of mixtures in both the gaseous and liquid phases, and show that the derived chemical potentials are accurate and reproducible. The method is recommended for obtaining chemical potentials in mixtures that are required as input in a grand canonical ensemble simulation.

Grand canonical Monte Carlo using solvent repacking: Application to phase behavior of hard disk mixtures.

A new "solvent repacking Monte Carlo" strategy for performing grand canonical ensemble simulations in condensed phases is introduced and applied to the study of hard-disk systems. The strategy is based on the configuration-bias approach, but uses an auxiliary biasing potential to improve the efficiency of packing multiple solvent particles in the cavity formed by removing one large solute. The method has been applied to study the coexistence of ordered and isotropic phases in three binary mixtures of hard disks with a small mole fraction (xL < 0.02) of the larger "solute" component. A chemical potential of 12.81 ± 0.01 kBT was found to correspond to the freezing transition of the pure hard disk "solvent." Simulations permitted the study of partitioning of large disks between ordered and isotropic phases, which showed a distinct non-monotonic dependence on size; the isotropic phase was enriched approximately 10-fold, 20-fold, and 5-fold over the coexisting ordered phases at diameter ratios d = 1.4, 2.5, and 3, respectively. Mixing of large and small disks within both phases near coexistence was strongly non-ideal in spite of the dilution. Structures of systems near coexistence were analyzed to determine correlations between large disks' positions within each phase, the orientational correlation length of small disks within the fluid phases, and the nature of translational order in the ordered phase. The analyses indicate that the ordered phase coexists with an isotropic phase resembling a nanoemulsion of ordered domains of small disks, with large disks enriched at the disordered domain interfaces.