Spatial updating grand canonical Monte Carlo algorithms for fluid simulation: generalization to continuous potentials and parallel implementation

Spatial updating grand canonical Monte Carlo algorithms are generalizations of random and sequential updating algorithms for lattice systems to continuum fluid models. The elementary steps, insertions or removals, are constructed by generating points in space either at random (random updating) or in a prescribed order (sequential updating). These algorithms have previously been developed only for systems of impenetrable spheres for which no particle overlap occurs. In this work, spatial updating grand canonical algorithms are generalized to continuous, soft-core potentials to account for overlapping configurations. Results on two- and three-dimensional Lennard-Jones fluids indicate that spatial updating grand canonical algorithms, both random and sequential, converge faster than standard grand canonical algorithms. Spatial algorithms based on sequential updating not only exhibit the fastest convergence but also are ideal for parallel implementation due to the absence of strict detailed balance and the nature of the updating that minimizes interprocessor communication. Parallel simulation results for three-dimensional Lennard-Jones fluids show a substantial reduction of simulation time for systems of moderate and large size. The efficiency improvement by parallel processing through domain decomposition is always in addition to the efficiency improvement by sequential updating.     

Modeling desorption of fluids from disordered mesoporous materials

The desorption mechanism of fluids in disordered mesoporous glasses is studied by Monte Carlo simulations of a coarse-grained lattice model with realistic matrix configurations representative of Vycor. Two methods of simulation are considered: grand canonical ensemble Monte Carlo simulations and dynamic Monte Carlo simulations which mimic the diffusion of the fluid in and out of the material using Kawasaki dynamics. In the grand canonical simulations, cavitation via nucleation of bubbles inside the pores plays the dominant role in determining the fluid configurations along the desorption isotherm. The Kawasaki dynamics simulations indicate that such configurations are achieved dynamically via the gradual advancement of macroscopic front interfaces toward the interior. This is made possible by the bubble nucleation mechanism operating on a length scale that is determined by both the typical pore size and the strength of the solid-fluid interaction.     

Simulating prescribed particle densities in the grand canonical ensemble using iterative algorithms

We present two efficient iterative Monte Carlo algorithms in the grand canonical ensemble with which the chemical potentials corresponding to prescribed (targeted) partial densities can be determined. The first algorithm works by always using the targeted densities in the kT log(rho(i)) (ideal gas) terms and updating the excess chemical potentials from the previous iteration. The second algorithm extrapolates the chemical potentials in the next iteration from the results of the previous iteration using a first order series expansion of the densities. The coefficients of the series, the derivatives of the densities with respect to the chemical potentials, are obtained from the simulations by fluctuation formulas. The convergence of this procedure is shown for the examples of a homogeneous Lennard-Jones mixture and a NaCl-CaCl(2) electrolyte mixture in the primitive model. The methods are quite robust under the conditions investigated. The first algorithm is less sensitive to initial conditions.     

Incorporation of inhomogeneous ion diffusion coefficients into kinetic lattice grand canonical monte carlo simulations and application to ion current calculations in a simple model ion channel

To deal with inhomogeneous diffusion coefficients of ions without altering the lattice spacing in the kinetic lattice grand canonical Monte Carlo (KLGCMC) simulation, an algorithm that incorporates diffusion coefficient variation into move probabilities is proposed and implemented into KLGCMC calculations. Using this algorithm, the KLGCMC simulation method is applied to the calculation of ion currents in a simple model ion channel system. Comparisons of ion currents and ion concentrations from these simulations with Poisson-Nernst-Planck (PNP) results show good agreement between the two methods for parameters where the latter method is expected to be accurate.     

A continuum, O(N) Monte Carlo algorithm for charged particles

We introduce a Monte Carlo algorithm for the simulation of charged particles moving in the continuum. Electrostatic interactions are not instantaneous as in conventional approaches, but are mediated by a constrained, diffusing electric field on an interpolating lattice. We discuss the theoretical justifications of the algorithm and show that it efficiently equilibrates model electrolytes and polar fluids. In order to reduce lattice artifacts that arise from the interpolation of charges to the grid we implement a local, dynamic subtraction algorithm. This dynamic scheme is completely general and can also be used with other Coulomb codes, such as multigrid based methods.     

Accelerating molecular simulations by reversible mapping between local minima

A new framework is presented for performing Monte Carlo simulations of condensed matter based on a recently developed bijective mapping between local energy minima. The framework is used to implement a range of new multiparticle Monte Carlo moves, which are investigated by simulating atomic Lennard-Jones fluids in the canonical and grand canonical ensembles. Important aspects of the simulation protocol and their effect on performance are analyzed in detail. Using the mapping accelerates the simulations by many orders of magnitude when compared to the equivalent moves without the mapping, and leads to particularly efficient configurational sampling at low temperatures and high densities. The method appears to be suitable for adapting to quantitative simulations of more complex molecular systems over long effective time scales.     

Modified Monte Carlo simulation of adsorption in a spatially inhomogeneous porous system

A modified Monte Carlo method in conjunction with the canonical and grand canonical ensembles is proposed for simulating adsorption in spatially inhomogeneous porous systems. Unlike the traditional Monte Carlo simulation in terms of the grand canonical ensemble, the simulation for the regions of pore space having no direct communication with the bulk phase is performed in local conditions of the canonical ensemble.     

Monte carlo simulations of micellization in model ionic surfactants: application to sodium dodecyl sulfate

A lattice model for ionic surfactants with explicit counterions is proposed for which the micellization behavior can be accurately determined from grand canonical Monte Carlo simulations. The model is characterized by a few parameters that can be adjusted to represent various linear surfactants with ionic headgroups. The model parameters have a clear physical interpretation and can be obtained from experimental data unrelated to micellization, namely, geometric information and solubilities of tail segments. As a specific example, parameter values for sodium dodecyl sulfate were obtained by optimizing for the solubility of hydrocarbons in water and the structural properties of dodecane. The critical micelle concentration (cmc), average aggregation number, degree of counterion binding, and their dependence on temperature were determined from histogram reweighting grand canonical Monte Carlo simulations and were compared to experimental results. The model gives the correct trend and order of magnitude for all quantities but underpredicts the cmc and aggregation number. We suggest ways to modify the model that may improve agreement with experimental values.     

Parallel Markov chain Monte Carlo simulations

With strict detailed balance, parallel Monte Carlo simulation through domain decomposition cannot be validated with conventional Markov chain theory, which describes an intrinsically serial stochastic process. In this work, the parallel version of Markov chain theory and its role in accelerating Monte Carlo simulations via cluster computing is explored. It is shown that sequential updating is the key to improving efficiency in parallel simulations through domain decomposition. A parallel scheme is proposed to reduce interprocessor communication or synchronization, which slows down parallel simulation with increasing number of processors. Parallel simulation results for the two-dimensional lattice gas model show substantial reduction of simulation time for systems of moderate and large size.     

Thermodynamic characterization of fluids confined in heterogeneous pores by monte carlo simulations in the grand canonical and the isobaric-isothermal ensembles

Materials presenting nanoscale porosity are able to condense gases in their structure. This "capillary condensation" phenomenon has been studied for more than one century. Theoretical models help to understand experimental results but fail in explaining all experimental features. Most of the time, the difficulties in making quantitative or even qualitative predictions are due to the geometric complexity of the porous materials, such as large pore size distribution, chemical heterogeneities, or pore interconnections. Numerical calculations (lattice gas models or molecular simulations) are of considerable interest to calculate the adsorption properties of a fluid confined in a porous model with characteristic sizes up to several tens of nanometers. For instance, the grand canonical Monte Carlo method allows one to compute the average amount of fluid adsorbed in the porous model as a function of the temperature and the chemical potential of the fluid. However, the grand potential, necessary for a complete characterization of the system, is not a direct output of the algorithm. It is shown in this paper that the use of the isobaric-isothermal (NPT) ensemble allows one to circumvent this problem; that is, it is possible to get in one single Monte Carlo run the absolute grand potential for any given thermodynamic state of the fluid. A simplified thermodynamic integration scheme is then used to evaluate the grand potential over the whole isotherm branch passing through this initially given point. Since the usual NPT technique is a priori limited to homogeneous pores, it is proposed, for the first time, to generalize this procedure to a pore presenting a chemical heterogeneity along its axis. The new method gives the same results as the previous for homogeneous pores and allows new predictions for chemically heterogeneous pores. Comparison with the full integration scheme shows that the proposed direct calculation is faster since it avoids multiple Monte Carlo runs and more precise because it avoids the possible cumulative errors of the integration procedure.     

Thermodynamic properties of lattice hard-sphere models

Thermodynamic properties of several lattice hard-sphere models were obtained from grand canonical histogram- reweighting Monte Carlo simulations. Sphere centers occupy positions on a simple cubic lattice of unit spacing and exclude neighboring sites up to a distance sigma. The nearestneighbor exclusion model, sigma = radical2, was previously found to have a second-order transition. Models with integer values of sigma = 1 or 2 do not have any transitions. Models with sigma = radical3 and sigma = 3 have weak first-order fluid-solid transitions while those with sigma = 2 radical2, 2 radical3, and 3 radical2 have strong fluid-solid transitions. Pressure, chemical potential, and density are reported for all models and compared to the results for the continuum, theoretical predictions, and prior simulations when available.     

Boltzmann bias grand canonical Monte Carlo

We derive an efficient method for the insertion of structured particles in grand canonical Monte Carlo simulations of adsorption in very confining geometries. We extend this method to path integral simulations and use it to calculate the isotherm of adsorption of hydrogen isotopes in narrow carbon nanotubes (two-dimensional confinement) and slit pores (one-dimensional confinement) at the temperatures of 20 and 77 K, discussing its efficiency by comparison to the standard path integral grand canonical Monte Carlo algorithm. We use this algorithm to perform multicomponent simulations in order to calculate the hydrogen isotope selectivity for adsorption in narrow carbon nanotubes and slit pores at finite pressures. The algorithm described here can be applied to the study of adsorption of real oligomers and polymers in narrow pores and channels.     

Lattice gas 2D/3D equilibria: chemical potentials and adsorption isotherms with correct critical points

A priori information is used to derive the chemical potential as a function of density and temperature for 2D and 3D lattice systems. The functional form of this equation of state is general in terms of lattice type and dimensionality, though it contains critical temperature and critical density as parameters which depend on lattice type and dimensionality. The adsorption isotherm is derived from equilibrium between two-dimensional and three-dimensional phases. Theoretical predictions are in excellent agreement with grand canonical Monte Carlo simulations.     

A computer simulation study of stick-slip transitions in water films confined between model hydrophilic surfaces. 1. Monolayer films

The shear behavior of monolayer water films confined in a slit-like pore between hydrophilic walls is simulated in the quasistatic regime using the grand canonical Monte Carlo technique. Each wall is represented as a hexagonal lattice of force sites that interact with water through an orientation-dependent hydrogen-bonding potential. When the walls are in registry, the water oxygen atoms form either a crystal- or fluid-like structure, depending on the period of the wall's lattice. In both cases, however, the monolayer structure is orientationally disordered. Both the crystal- and fluid-like monolayers prove to be capable of experiencing well-defined stick-slip transitions, with the largest yield stress occurring in the crystal-like case. Beyond the yield point, the crystal-like monolayers "melt", but their structure and molecular motion differ in many respects from those characteristic of normal fluids.     

Molecular Modeling of Adsorption in Activated Carbon: Comparison of Monte Carlo Simulations with Experiment

We present Mont Carlo computer simulation results for a molecular model of fluids adsorbed in porous carbon materials. The model carbon used is based on the platelet model for carbon of Segarra and Glandt (1994). The model we use has a single basal plane per platelet and the structure is isotropic, disordered, with weak short-range correlations between the platelets. We have performed grand canonical Monte Carlo simulations of the adsorption isotherms for methane, ethane, and their mixtures in this model carbon. We find generally good agreement with experimental and the mixture results are quite accurately described by the ideal adsorbed solution theory. An exception to this is the behavior for the mixture at the highest pressures. In this case the experimental data show significant deviations from ideal adsorbed solution theory and the simulation results.     

Density functional theory of adsorption in pillared slit-like pores

We propose a density functional theory to describe adsorption of Lennard-Jones fluid in pillared slit like pores. Specifically, the pillars are built of chains that are bonded by their ends to the opposite pore walls. The approach we propose combines theory of quenched-annealed systems and theory of nonuniform fluids involving chain molecules. We compare the results of theoretical predictions with grand canonical ensemble Monte Carlo simulations and compute theoretical capillary condensation phase diagrams for several model systems.