Transition Matrix Monte Carlo Simulations Research Articles

Efficiency Comparison of Single- and Multiple-Macrostate Grand Canonical Ensemble Transition-Matrix Monte Carlo Simulations

Recent interest in parallelizing flat-histogram transition-matrix Monte Carlo simulations in the grand canonical ensemble, due to its demonstrated effectiveness in studying phase behavior, self-assembly and adsorption, has led to the most extreme case of single-macrostate simulations, where each macrostate is simulated independently with ghost particle insertions and deletions. Despite their use in several studies, no efficiency comparisons of these single-macrostate simulations have been made with multiple-macrostate simulations. We show that multiple-macrostate simulations are up to 3 orders of magnitude more efficient than single-macrostate simulations, which demonstrates the remarkable efficiency of flat-histogram biased insertions and deletions, even with low acceptance probabilities. Efficiency comparisons were made for supercritical fluids and vapor-liquid equilibrium of bulk Lennard-Jones and a three-site water model, self-assembling patchy trimer particles and adsorption of a Lennard-Jones fluid confined in a purely repulsive porous network, using the open source simulation toolkit FEASST. By directly comparing with a variety of Monte Carlo trial move sets, this efficiency loss in single-macrostate simulations is attributed to three related reasons. First, ghost particle insertions and deletions in single-macrostate simulations incur the same computational expense as grand canonical ensemble trials in multiple-macrostate simulations, yet ghost trials do not reap the sampling benefit from propagating the Markov chain to a new microstate. Second, single-macrostate simulations lack macrostate change trials that are biased by the self-consistently converging relative macrostate probability, which is a major component of flat histogram simulations. Third, limiting a Markov chain to a single macrostate reduces sampling possibilities. Existing parallelization methods for multiple-macrostate flat-histogram simulations are shown to be more efficient than parallel single-macrostate simulations by approximately an order of magnitude or more in all systems investigated.

Effect of shape anisotropy on the precipitation of dimeric nanoparticles.

We use grand canonical transition matrix Monte Carlo simulations to study the precipitation of dimeric nanoparticles. The dimers are composed of two particles having different chemical features and separated by a fixed distance. The non-attractive and attractive parts of the dimer are modeled using hard-sphere and square-well potentials, respectively. The shape anisotropy is altered by changing the relative sizes of the two particles. We observe that the stability of the nanosuspension increases with the increase in the size of the non-attractive part of the dimer. The precipitates of dimers having larger non-attractive parts have lower packing densities, contain large cavities, and show evidence of self-assembly in the bulk and on the surface. We also use the results from our simulations and the classical nucleation theory to study the kinetics of precipitation. At a given temperature and relative supersaturation, the rate of homogeneous nucleation increases with the increase in the size of the non-attractive parts. Finally, we use an example to show how our results can guide the design of nanosuspensions containing chemically anisotropic dimers that are stable under particular conditions.

Prediction of Thermodynamic Properties of Levulinic Acid via Molecular Simulation Techniques.

Second-generation biofuels are a complex mixture of organic compounds that can be further processed to hydrocarbon fuels and other valuable chemicals. One such chemical is levulinic acid (IUPAC name: 4-oxo pentanoic acid), which is a highly versatile ketoacid obtained from cellulose present in agricultural byproducts. For oxygen-containing compounds that decompose at elevated temperatures and pressures, determining the vapor–liquid equilibria data at high temperatures via the experimental route may be challenging. The molecular simulation approach is a cost-effective tool to obtain the necessary data while also allowing us to understand the microscopic origins of macroscopic observable properties. We have employed the transferable potential for phase equilibria-united atom force field to describe the interactions in this system with the parameters for a torsional interaction that are not reported in the literature (levulinic acid is a ketoacid) being determined from density functional theory calculations. We have verified our parameterization via density computations in the isothermal–isobaric ensemble and by comparing our simulation results with the corresponding data from experiments reported in the literature. We have performed grand-canonical transition-matrix Monte Carlo simulations in the temperature range from 580 to 690 K to estimate the vapor–liquid coexistence curves in the temperature–density plane and the Clapeyron plots. From this data, the critical point (TC = 755 K, ρC = 285.4 kg/m3, and PC = 30.57 bar) has been estimated, and this may be used as input to the equations of state employed in process simulation software for design of industrial separation processes including those for “biorefining”. As levulinic acid is a “ketoacid”, hydrogen bonding occurs, and the liquid phase structure has also been studied using radial distribution functions.

Study of fluid phase equilibria of triangle-well mixtures in bulk and under confinement using grand canonical transition matrix Monte Carlo

Grand canonical transition matrix Monte Carlo simulations are used to determine the fluid phase diagrams, at two different temperatures, in bulk and under confinement of weakly attractive slit pores, for four different types of binary mixtures: mixture 1 is symmetric, with a weakened interaction between the unlike species (unlike the other three mixtures), mixture 2 with asymmetry in size of the molecules, mixture 3 with asymmetry in the depth of the potential well and, finally, mixture 4 where both the size and well-depth are different. The triangle-well (TW) potential is employed to model the fluid-fluid and wall-fluid interactions as this potential provides simple, qualitatively accurate and practical representation of the interactions present in real systems. The simple nature of the TW models enables us to obtain a qualitative understanding of the effects of the fluid-fluid potential parameters and the effect of confinement on the vapor-liquid coexistence behavior of the four different types of TW fluid mixtures described above. Significant differences in the vapor-liquid coexistence behavior of the confined fluid from the bulk phase behavior are observed. Comparison of the bulk and confined systems indicate that the intersection of the liquid-liquid phase envelope with the vapor-liquid equilibrium (VLE) curve observed in the bulk phase of mixture 1 does not occur under confinement with the heteroazeotrope being replaced by VLE showing an azeotrope at lower pressures and liquid-liquid coexistence at higher pressure. A significant decrease in the area of the vapor + liquid coexistence region is noted for mixture 2. The supercritical behavior of component 2 is noted for the confined systems in mixtures 3 and 4 at both temperatures, however, the pressure-composition diagram in the bulk system shows supercritical behavior only at the higher temperature. Canonical Monte Carlo simulations were subsequently performed on the confined systems to obtain the energy values and the density and composition profiles along the pore width at the equilibrium points. The latter data have been utilized to explain the phase diagrams of the confined mixtures and correlate the differences in the phase diagrams to the variation in the potential parameters.

Critical temperature estimation of bulk and confined atomic fluid using vapour−liquid interfacial free energy

The critical temperatures of bulk and confined atomic fluids are investigated using values of the vapour−liquid interfacial free energy of coexistence obtained from grand-canonical transition-matrix Monte Carlo simulations using a histogram reweighting technique. The temperature corresponding to zero interfacial free energy of coexistence is the estimated critical temperature for the system under investigation. Slit width of confined atomic fluid for this investigation is varied from 40 to 1 fluid particle diameter. The obtained critical temperatures have shown nonlinear monotonic trend with inverse of slit width. Moreover, five different linear regimes of critical temperature are observed in the studied range of slit width. Interestingly, in the slit width range of less than two fluid particle diameters, the critical temperature approaches two-dimensional value and remains approximately indifferent with a decrease in slit width up to one fluid particle diameter. This investigation also reveals that the critical temperature of bulk and confined atomic fluid estimated using the vapour−liquid interfacial free energy of coexistence is within reasonable accuracy with that obtained using the simplified form of scaling law.

Predicting Protein Interactions of Concentrated Globular Protein Solutions Using Colloidal Models.

Protein interactions of α-chymotrypsinogen A (aCgn) were quantified using light scattering from low to high protein concentrations. Static light scattering (SLS) was used to determine the excess Rayleigh ratio (Rex) and osmotic second virial coefficients (B22) as a function of pH and total ionic strength (TIS). Repulsive (attractive) protein-protein interactions (PPI) were observed at pH 5 (pH 7), with decreasing repulsions (attractions) upon increasing TIS. Simple colloidal potential of mean force models (PMF) that account for short-range nonelectrostatic attractions and screened electrostatic interactions were used to fit model parameters from data for B22 vs TIS at both pH values. The parameters and PMF models from low-concentration conditions were used as the sole input to transition matrix Monte Carlo simulations to predict high concentration Rex behavior. At conditions where PPI are repulsive to slightly attractive, experimental Rex data at high concentrations could be predicted quantitatively by the simulations. However, accurate predictions were challenging when PPI were strongly attractive due to strong sensitivity to changes in PMF parameter values. Additional simulations with higher-resolution coarse-grained molecular models suggest an approach to qualitatively predict cases when anisotropic surface charge distributions will lead to overall attractive PPI at low ionic strength, without assumptions regarding electrostatic "patches" or multipole expansions.

Corresponding states vapour–liquid phase equilibria of confined square–well fluid

Corresponding states vapour–liquid phase equilibria of confined square-well fluid are studied by means of grand-canonical transition-matrix Monte Carlo simulation and histogram reweighting method. In this study, square-well fluid is considered under hard and attractive slit pore confinements ranging from 1.5 to 40 molecular diameters. Corresponding states vapour–liquid phase coexistence envelopes display insignificant effect of wall−fluid interaction for slit pore confinements ranging from 1.5 to 3 molecular diameters. On the other hand, significant effect of wall−fluid interaction on the corresponding state coexistence envelope is observed for slit pore confinements ranging from 4 to 40 molecular diameters. Moreover, at a given slit width, shrinking in corresponding state coexistence envelope is observed with increase in the wall−fluid interaction. However, in the larger slit pore width of 30 to 40 molecular diameters, shrinking in the corresponding state vapour–liquid coexistence envelopes become indifferent with the stronger wall−fluid interactions studied in this work. Structural behaviour of coexisting phases in slit pores are also investigated through local density profiles, to understand the overall behaviour of corresponding states coexistence envelopes. Fluctuating positive and negative deviations in the corresponding state spreading pressure with respect to corresponding bulk value is observed for studied wall−fluid interactions and slit pore confinements.

Modulus-pressure equation for confined fluids.

Ultrasonic experiments allow one to measure the elastic modulus of bulk solid or fluid samples. Recently such experiments have been carried out on fluid-saturated nanoporous glass to probe the modulus of a confined fluid. In our previous work [G. Y. Gor et al., J. Chem. Phys., 143, 194506 (2015)], using Monte Carlo simulations we showed that the elastic modulus K of a fluid confined in a mesopore is a function of the pore size. Here we focus on the modulus-pressure dependence K(P), which is linear for bulk materials, a relation known as the Tait-Murnaghan equation. Using transition-matrix Monte Carlo simulations we calculated the elastic modulus of bulk argon as a function of pressure and argon confined in silica mesopores as a function of Laplace pressure. Our calculations show that while the elastic modulus is strongly affected by confinement and temperature, the slope of the modulus versus pressure is not. Moreover, the calculated slope is in a good agreement with the reference data for bulk argon and experimental data for confined argon derived from ultrasonic experiments. We propose to use the value of the slope of K(P) to estimate the elastic moduli of an unknown porous medium.

Fluid phase equilibria of triangle-well fluids confined inside slit pores: A transition matrix Monte Carlo simulation study

The phase behaviour of triangle-well (TW) fluids confined inside three types of slit pores: slit pores with or without uniform attraction and slit pores with square-well (SW) type pore-fluid interactions is studied using molecular simulations. The TW potential is chosen to model the fluid-fluid interactions as this model provides a good compromise between realism and simplicity. The simple nature of the models enables us to obtain a qualitative understanding of the effects of confinement on the vapour liquid coexistence curves (VLCCs) of the TW fluids with variable well-width. We note the significant differences in the VLCCs of the confined fluid from the bulk phase behaviour; and also the dependence of the VLCCs on the type of slit pore. These differences have been explained on the basis of the plots of local density along the axis perpendicular to the walls of the slit pore.

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.

Prediction of fluid-phase behavior of symmetrical binary Yukawa fluids using transition matrix Monte Carlo

Binary fluid mixtures display a wide array of fluid phase behavior ranging from simple vapor-liquid equilibrium diagrams to more complicated topology like azeotropy and heteroazeotropy. As a first step towards understanding the phenomena of equilibrium of binary molecular systems, the properties of binary monoatomic fluids have been studied in this work. An additional simplification made in this work is that of a symmetrical binary system where all similar molecules interact via the hard-core Yukawa (HCY) potential, u(r), while the dissimilar molecules interact via a potential, δu(r) where δ is a scalar parameter. The value of δ here is 0.75, which leads to the dissimilar molecules showing a tendency to dislike one another and the resulting phase diagram is such that the mixing-demixing line intersects the vapor-liquid equilibrium curve away from the vapor-liquid critical point. The effect of the variable potential range parameter of the HCY potential on the topology of the phase diagrams is investigated using grand canonical transition matrix Monte Carlo simulations and the conditions of temperature and pressure at which the system exhibits azeotropy and heteroazeotropy are ascertained. We find that the densities of the mixtures, as predicted by our simulation at equimolar concentrations, are in close agreement with the self-consistent Ornstein-Zernike approximation results of Schöll-Paschinger and co-workers.

Estimating the Asymptotics of Solid Partitions

We study the asymptotic behavior of solid partitions using transition matrix Monte Carlo simulations. If $p_3(n)$ denotes the number of solid partitions of an integer $n$, we show that $\lim_{n\rightarrow\infty} n^{-3/4} \log p_3(n)\sim 1.822\pm 0.001$. This shows clear deviation from the value $1.7898$, attained by MacMahon numbers $m_3(n)$, that was conjectured to hold for solid partitions as well. In addition, we find estimates for other sub-leading terms in $\log p_3(n)$. In a pattern deviating from the asymptotics of line and plane partitions, we need to add an oscillatory term in addition to the obvious sub-leading terms. The period of the oscillatory term is proportional to $n^{1/4}$, the natural scale in the problem. This new oscillatory term might shed some insight into why partitions in dimensions greater than two do not admit a simple generating function.

A comparative study of critical temperature estimation of atomic fluid and chain molecules using fourth-order Binder cumulant and simplified scaling laws

Grand-canonical transition-matrix Monte Carlo simulation with histogram reweighting and finite-size scaling technique are used to calculate fourth-order Binder cumulant of order parameter along the vapour–liquid coexistence line to calculate the critical temperature of bulk and confined square-well fluid in slit pore of two pore sizes. Further, this approach is utilised to estimate the critical temperatures of relatively more complex fluids such as n-alkanes confined in graphite and mica slit pores of different slit widths. The estimated critical temperatures are compared with the critical temperature obtained for the same systems using simplified form of the scaling law. This investigation reveals that critical temperatures of simple and complex fluids in bulk state and under confinement, estimated using the scaling law, are within reasonable accuracy with that obtained using more accurate and rigorous approach of fourth-order Binder's cumulant.

Surface phase transitions of multiple-site associating fluids

Surface phase transitions of Lennard–Jones (LJ) based two- and four-site associating fluids have been studied for various associating strengths using grand-canonical transition matrix Monte Carlo simulations. Our results suggest that, in the case of a smooth surface, represented by a LJ 9-3-type potential, multiple-site associating fluids display a prewetting transition within a certain temperature range. However, the range of the prewetting transition decreases with increasing associating strength and increasing number of sites on the fluid molecules. With the addition of associating sites on the surface, a quasi-2D vapor–liquid transition may appear, which is observed at a higher surface site density for weaker associating fluids. The prewetting transition at lower associating strength is found to shift towards the quasi-2D vapor–liquid transition with increasing surface site density. However, for highly associating fluids, the prewetting transition is still intact, but shifts slightly towards the lower temperature range. Adsorption isotherms, chemical potentials and density profiles are used to characterize surface phase transitions.

Vapour–liquid phase equilibria of simple fluids confined in patterned slit pores

We present the influence of surface heterogeneity on the vapour–liquid phase behaviour of square-well fluids in slit pores using grand-canonical transition-matrix Monte Carlo simulations along with the histogram-reweighting method. Properties such as phase coexistence envelopes, critical properties and local density profiles of the confined SW fluid are reported for chemically and physically patterned slit surfaces. It is observed that in the chemically patterned pores, fluid–fluid and surface attraction parameters along with the width of attractive and inert stripes play fundamentally different roles in the phase coexistence and critical properties. On the other hand, pillar gap and height significantly affect the vapour–liquid equilibria in the physically patterned slit pores. We also present the effect of chemically and physically patterned slit surfaces on the spreading pressure.

On the Use of Excess Entropy Scaling to Describe the Dynamic Properties of Water

We use molecular dynamics and transition-matrix Monte Carlo simulation to study the ability of entropy scaling relationships to describe kinetic properties of the extended simple point charge water model. We examine translational and rotational diffusivities, a characteristic relaxation time for rotational motion, and a collective relaxation time stemming from analysis of the coherent intermediate scattering function. We consider both the thermodynamic excess entropy and the contribution to the two-body excess entropy related to center-of-mass correlations as scaling variables. Calculations are performed over a broad range of conditions that span from the dense supercooled liquid regime to the critical region. We find that the thermodynamic excess entropy serves as a suitable metric for describing reduced transport properties for state conditions corresponding to temperatures above the onset of water's structurally anomalous region, defined by states points for which the excess entropy increases upon compression at constant temperature. In contrast, the aforementioned two-body contribution to the excess entropy cannot be used to quantitatively predict kinetic properties over the wide range of conditions explored here. For state points above the onset temperature of the structurally anomalous region, reduced transport property data collapse onto common curves when expressed as a function of the thermodynamic excess entropy. Below this temperature, data fall onto isochore-specific curves. Our results show a relatively strong correlation between the translational diffusivity and excess entropy and a noticeably weaker correlation between rotational mobility and excess entropy.

Phase equilibria of molecular fluids via hybrid Monte Carlo Wang–Landau simulations: Applications to benzene and n-alkanes

In recent years, powerful and accurate methods, based on a Wang-Landau sampling, have been developed to determine phase equilibria. However, while these methods have been extensively applied to study the phase behavior of model fluids, they have yet to be applied to molecular systems. In this work, we show how, by combining hybrid Monte Carlo simulations in the isothermal-isobaric ensemble with the Wang-Landau sampling method, we determine the vapor-liquid equilibria of various molecular fluids. More specifically, we present results obtained on rigid molecules, such as benzene, as well as on flexible chains of n-alkanes. The reliability of the method introduced in this work is assessed by demonstrating that our results are in excellent agreement with the results obtained in previous work on simple fluids, using either transition matrix or conventional Monte Carlo simulations with a Wang-Landau sampling, and on molecular fluids, using histogram reweighting or Gibbs ensemble Monte Carlo simulations.

Vapor−Liquid Phase Coexistence, Critical Properties, and Surface Tension of Confined Alkanes

Configurational-bias grand-canonical transition-matrix Monte Carlo simulations are conducted to investigate various thermophysical properties, such as phase coexistence, critical properties, density and orientation profiles of liquid and vapor phases, and vapor−liquid surface tension of methane, ethane, propane, n-butane, and n-octane in bulk and slit pores of graphite and mica surfaces. An exponential-6 (exp-6) model is used for the investigation of normal alkanes with a cutoff radius, 15 A. It is found that the surface tension of the bulk n-alkane, C1−C4, based on a truncated exp-6 model, agrees reasonably well with the experimental data. Critical properties are reported by means of the rectilinear diameter approach and least-squares technique. The shift in the critical temperature under confinements follows more than two linear regimes with an inverse in the slit width, as the slit width approaches the two-dimensional limit. This is contrary to what has been previously reported in the literature. The b...

Phase Behavior of Model Confined Fluids. Influence of Substrate−Fluid Interaction Strength

We examine the relationship between the macroscopic phase behavior of nanoconfined fluids and the nature of microscopic interactions between a confining substrate and fluid. Two model slit-pore systems are explored using grand canonical transition-matrix Monte Carlo simulation. One system consists of a square-well fluid interacting with a square-well substrate, and the other contains an embedded point charge model of lysozyme interacting with a mica surface. Fluid phase diagrams are constructed for a broad range of substrate conditions. Our results indicate that one observes a maximum in the critical temperature of the fluid phase envelope upon variation of substrate strength for a given slit width. Both systems studied exhibit such maxima at intermediate wall strength. The physical rationale for this observation suggests that this behavior should be generally expected. We introduce two metrics that enable one to predict conditions that produce maxima in critical temperature. The first is related to the contact angle a fluid develops at a single confining substrate. The second is based upon virial coefficient information and requires knowledge of the substrate-fluid and fluid-fluid interaction potentials only.

Fluid phase behavior of a model colloid-polymer mixture: Influence of polymer size and interaction strength

We examine how the fluid-fluid phase behavior of a model colloid-polymer mixture evolves with variation of polymer size and/or interaction strength. Polymer-polymer interactions are approximated through Gaussian-core potentials while colloid-colloid and colloid-polymer interactions are assumed purely hard. Grand canonical transition matrix Monte Carlo simulation is used to construct binary liquid-liquid phase diagrams for a wide range of polymer sizes and interaction strengths. Overall, our results indicate that independently decreasing the polymer size or increasing the strength of repulsion between polymers have similar effects; these variations lead to an increase in the critical colloid concentration and decrease in the critical polymer concentration. An examination of the phase diagrams collected in this study reduced by their critical reservoir polymer concentration and critical colloid concentration reveals a single common phase envelope that is quantitatively similar to that for common models used to describe simple fluids. Finally, evaluation of the osmotic second virial coefficient in the vicinity of the critical point suggests that this quantity cannot be used as a reliable predictor for the conditions under which colloid-polymer mixtures phase separate.