Monte Carlo Simulation Methods for Computing Liquid–Vapor Saturation Properties of Model Systems

Flat-Histogram Monte Carlo Simulation of Water Adsorption in Metal-Organic Frameworks

Molecular simulations of water adsorption in porous materials often converge slowly due to sampling bottlenecks that follow from hydrogen bonding and, in many cases, the formation of water clusters. These effects may be exacerbated in metal-organic framework (MOF) adsorbents, due to the presence of pore spaces (cages) that promote the formation of discrete-size clusters and hydrophobic effects (if present), among other reasons. In Grand Canonical Monte Carlo (MC) simulations, these sampling challenges are typically manifested by low MC acceptance ratios, a tendency for the simulation to become stuck in a particular loading state (i.e., macrostates), and the persistence of specific clusters for long periods of the simulation. We present simulation strategies to address these sampling challenges, by applying flat-histogram MC (FHMC) methods and specialized MC move types to simulations of water adsorption. FHMC, in both Transition-matrix and Wang-Landau forms, drives the simulation to sample relevant macrostates by incorporating weights that are self-consistently adjusted throughout the simulation and generate the macrostate probability distribution (MPD). Specialized MC moves, based on aggregation-volume bias and configurational bias methods, separately address low acceptance ratios for basic MC trial moves and specifically target water molecules in clusters; in turn, the specialized MC moves improve the efficiency of generating new configurations which is ultimately reflected in improved statistics collected by FHMC. The combined strategies are applied to study the adsorption of water in CuBTC and ZIF-8 at 300 K, through examination of the MPD and the adsorption isotherm generated by histogram reweighting. A key result is the appearance of nontrivial oscillations in the MPD, which we show to be associated with water clusters in the adsorption system. Additionally, we show that the probabilities of certain clusters become similar in value near the boundaries of the isotherm hysteresis loop, indicating a strong connection between cluster formation/destruction and the thermodynamic limits of stability. For a hydrophobic MOF, the FHMC results show that the phase transition from low density to high density is suppressed to water pressure far above the bulk-fluid saturation pressure; this is consistent with results presented elsewhere. We also compare our FHMC simulation isotherm to one measured by a different technique but with ostensibly the same molecular interactions and comment on observed differences and the need for follow-up work. The simulation strategies presented here can be applied to the simulation of water in other MOFs using heuristic guidelines laid out in our text, which should facilitate the more consistent and efficient simulation of water adsorption in porous materials in future applications.

Comparison of Force Fields for the Prediction of Thermophysical Properties of Long Linear and Branched Alkanes

The prediction of thermophysical properties at extreme conditions is an important application of molecular simulations. The quality of these predictions primarily depends on the quality of the employed force field. In this work, a systematic comparison of classical transferable force fields for the prediction of different thermophysical properties of alkanes at extreme conditions, as they are encountered in tribological applications, was carried out using molecular dynamics simulations. Nine transferable force fields from three different classes were considered (all-atom, united-atom, and coarse-grained force fields). Three linear alkanes (n-decane, n-icosane, and n-triacontane) and two branched alkanes (1-decene trimer and squalane) were studied. Simulations were carried out in a pressure range between 0.1 and 400 MPa at 373.15 K. For each state point, density, viscosity, and self-diffusion coefficient were sampled, and the results were compared to experimental data. The Potoff force field yielded the best results.

Coupled Monte Carlo and Molecular Dynamics Simulations on Interfacial Properties of Antifouling Polymer Membranes

We use molecular simulation to study the wetting behavior of antifouling polymer-tethered membranes. We obtain the interfacial properties (e.g., contact angle) of water at various temperatures for five polymer membranes, including a base polysulfone (PSF) membrane and four other PSF membranes grafted with antifouling polymers (two poly(ethylene glycol) (PEG) tethers and two zwitterionic tethers). We implement a coupled Monte Carlo (MC)/molecular dynamics (MD) approach to determine the interface potentials of water on the membrane surfaces in an efficient manner. Within this method, short MC and MD simulations are performed in cycles to collect the surface excess free energy of a thin water film on polymer membrane surfaces. Simulation results show that the grafting of zwitterionic tethers provides a more significant enhancement in the hydrophilicity of the PSF membrane than that of the PEG tethers. Water completely wets the surface of zwitterionic polymer membranes.

Entropy scaling close to criticality: From simple to metallic systems

Entropy has recently drawn considerable interest both as a marker to detect the onset of phase transitions and as a reaction coordinate, or collective variable, to span phase transition pathways. We focus here on the behavior of entropy along the vapor-liquid phase coexistence and identify how the difference in entropy between the two coexisting phases vary in ideal and metallic systems along the coexistence curve. Using flat-histogram simulations, we determine the thermodynamic conditions of coexistence, critical parameters, including the critical entropy, and entropies along the binodal. We then apply our analysis to a series of systems that increasingly depart from ideality and adopt a metal-like character, through the gradual onset of the Friedel oscillation in an effective pair potential, and for a series of transition metals modeled with a many-body embedded-atoms force field. Projections of the phase boundary on the entropy-pressure and entropy-temperature planes exhibit two qualitatively different behaviors. While all systems modeled with an effective pair potential lead to an ideal-like behavior, the onset of many-body effects results in a departure from ideality and a markedly greater exponent for the variation of the entropy of vaporization with temperature away from the critical temperature.

Entropy in Molecular Fluids: Interplay between Interaction Complexity and Criticality

Using flat-histogram simulations, we calculate the entropy of molecular fluids along the vapor-liquid phase boundary. Our simulation approach is based on the evaluation of the canonical and grand-canonical partition functions, which, in turn, provide access to entropy through the statistical mechanics formalism. The results allow us to determine the critical entropy of molecular fluids and to uncover that the transition occurs symmetrically from an entropic standpoint. This can best be seen through the patterns exhibited by the thermodynamic variables temperature and pressure when plotted against the entropy of the coexisting phases. This behavior is found to hold for apolar, quadrupolar, and dipolar fluids. Finally, we identify functional forms that characterize the relation between thermodynamic variables and entropy along the coexistence curve up to the critical point.

Parallel Prefetching for Canonical Ensemble Monte Carlo Simulations

In order to enable large-scale molecular simulations, algorithms must efficiently utilize multicore processors that continue to increase in total core count over time with relatively stagnant clock speeds. Although parallelized molecular dynamics (MD) software has taken advantage of this trend in computer hardware, single-particle perturbations with Monte Carlo (MC) are more difficult to parallelize than system-wide updates in MD using domain decomposition. Instead, prefetching reconstructs the serial Markov chain after computing multiple MC trials in parallel. Canonical ensemble MC simulations of a Lennard-Jones fluid with prefetching resulted in up to a factor of 1.7 speedup using 2 threads, and a factor of 3 speedup using 4 threads. Strategies for maximizing efficiency of prefetching simulations are discussed, including the potentially counterintuitive benefit of reduced acceptance probabilities. Determination of the optimal acceptance probability for a parallel simulation is simplified by theoretical prediction from serial simulation data. Finally, complete open-source code for parallel prefetch simulations was made available in the Free Energy and Advance Sampling Simulation Toolkit (FEASST).

Improving the efficiency of Monte Carlo simulations of ions using expanded grand canonical ensembles

While ionic liquids have promising applications as industrial solvents, predicting their fluid phase properties and coexistence remains a challenge. Grand canonical Monte Carlo simulation is an effective method for such predictions, but equilibration is hampered by the apparent requirement to insert and delete neutral sets of ions simultaneously in order to maintain charge neutrality. For relatively high densities and low temperatures, previously developed methods have been shown to be essential in improving equilibration by gradual insertion and deletion of these neutral sets of ions. We introduce an expanded ensemble approach which may be used in conjunction with these existing methods to further improve efficiency. Individual ions are inserted or deleted in one Monte Carlo trial rather than simultaneous insertion/deletion of neutral sets. We show how charge neutrality is maintained and show rigorous quantitative agreement between the conventional and the proposed expanded ensemble approaches, but with up to an order of magnitude increase in efficiency at high densities. The expanded ensemble approach is also more straightforward to implement than simultaneous insertion/deletion of neutral sets, and its implementation is demonstrated within open source software.

Effect of Carboxylic Acid on the Wetting Properties of a Model Water-Octane-Silica System

Monte Carlo simulations are employed to determine the effects of acetic acid on the wetting properties of a model water-octane-silica system. We first compute the bulk liquid-vapor saturation properties of pure acetic acid and subsequently explore the bulk liquid-liquid saturation properties of the ternary water-octane-acid system. We introduce an expanded ensemble approach to compute the coexistence properties of the ternary system. An interface potential approach is then used to capture the evolution of the wetting properties of the water-octane-silica system upon the addition of acetic acid. We track the change in the octane-water interfacial tension and the contact angle of water droplet on a silica substrate in a mother octane fluid over a range of acetic acid activities. The structure of the fluid, including the partitioning of acetic acid within the interfacial system, is also considered at several state points. We observe that acetic acid has a strong tendency to adsorb at the octane-water interface, resulting in a reduction in the octane-water interfacial tension. The response of the contact angle is more sensitive to the temperature and the hydrophilicity of the silica substrate.

Flat-Histogram Monte Carlo as an Efficient Tool To Evaluate Adsorption Processes Involving Rigid and Deformable Molecules

Monte Carlo simulations are the foundational technique for predicting thermodynamic properties of open systems where the process of interest involves the exchange of particles. Thus, they have been used extensively to computationally evaluate the adsorption properties of nanoporous materials and are critical for the in silico identification of promising materials for a variety of gas storage and chemical separation applications. In this work we demonstrate that a well-known biasing technique, known as "flat-histogram" sampling, can be combined with temperature extrapolation of the free energy landscape to efficiently provide significantly more useful thermodynamic information than standard open ensemble MC simulations. Namely, we can accurately compute the isosteric heat of adsorption and number of particles adsorbed for various adsorbates over an extremely wide range of temperatures and pressures from a set of simulations at just one temperature. We extend this derivation of the temperature extrapolation to adsorbates with intramolecular degrees of freedom when Rosenbluth sampling is employed. Consequently, the working capacity and isosteric heat can be computed for any given combined temperature/pressure swing adsorption process for a large range of operating conditions with both rigid and deformable adsorbates. Continuous thermodynamic properties can be computed with this technique at very moderate computational cost, thereby providing a strong case for its application to the in silico identification of promising nanoporous adsorbents.

Monte Carlo simulation of cylinders with short-range attractions

Cylindrical or rod-like particles are promising materials for the applications of fillers in nanocomposite materials and additives to control rheological properties of colloidal suspensions. Recent advances in particle synthesis allows for cylinders to be manufactured with short-ranged attractions to study the gelation as a function of packing fraction, aspect ratio and attraction strength. In order to aid in the analysis of small-angle scattering experiments of rod-like particles, computer simulation methods were used to model these particles with specialized Monte Carlo algorithms and tabular superquadric potentials. The attractive interaction between neighboring rods increases with the amount of locally-accessible surface area, thus leading to patchy-like interactions. We characterize the clustering and percolation of cylinders as the attractive interaction increases from the homogenous fluid at relatively low attraction strength, for a variety of aspect ratios and packing fractions. Comparisons with the experimental scattering results are also presented, which are in agreement.

Using reweighting and free energy surface interpolation to predict solid-solid phase diagrams

Many physical properties of small organic molecules are dependent on the current crystal packing, or polymorph, of the material, including bioavailability of pharmaceuticals, optical properties of dyes, and charge transport properties of semiconductors. Predicting the most stable crystalline form at a given temperature and pressure requires determining the crystalline form with the lowest relative Gibbs free energy. Effective computational prediction of the most stable polymorph could save significant time and effort in the design of novel molecular crystalline solids or predict their behavior under new conditions. In this study, we introduce a new approach using multistate reweighting to address the problem of determining solid-solid phase diagrams and apply this approach to the phase diagram of solid benzene. For this approach, we perform sampling at a selection of temperature and pressure states in the region of interest. We use multistate reweighting methods to determine the reduced free energy differences between T and P states within a given polymorph and validate this phase diagram using several measures. The relative stability of the polymorphs at the sampled states can be successively interpolated from these points to create the phase diagram by combining these reduced free energy differences with a reference Gibbs free energy difference between polymorphs. The method also allows for straightforward estimation of uncertainties in the phase boundary. We also find that when properly implemented, multistate reweighting for phase diagram determination scales better with the size of the system than previously estimated.

Temperature extrapolation of multicomponent grand canonical free energy landscapes

We derive a method for extrapolating the grand canonical free energy landscape of a multicomponent fluid system from one temperature to another. Previously, we introduced this statistical mechanical framework for the case where kinetic energy contributions to the classical partition function were neglected for simplicity [N. A. Mahynski et al., J. Chem. Phys. 146, 074101 (2017)]. Here, we generalize the derivation to admit these contributions in order to explicitly illustrate the differences that result. Specifically, we show how factoring out kinetic energy effects a priori, in order to consider only the configurational partition function, leads to simpler mathematical expressions that tend to produce more accurate extrapolations than when these effects are included. We demonstrate this by comparing and contrasting these two approaches for the simple cases of an ideal gas and a non-ideal, square-well fluid.

Multivariable extrapolation of grand canonical free energy landscapes

We derive an approach for extrapolating the free energy landscape of multicomponent systems in the grand canonical ensemble, obtained from flat-histogram Monte Carlo simulations, from one set of temperature and chemical potentials to another. This is accomplished by expanding the landscape in a Taylor series at each value of the order parameter which defines its macrostate phase space. The coefficients in each Taylor polynomial are known exactly from fluctuation formulas, which may be computed by measuring the appropriate moments of extensive variables that fluctuate in this ensemble. Here we derive the expressions necessary to define these coefficients up to arbitrary order. In principle, this enables a single flat-histogram simulation to provide complete thermodynamic information over a broad range of temperatures and chemical potentials. Using this, we also show how to combine a small number of simulations, each performed at different conditions, in a thermodynamically consistent fashion to accurately compute properties at arbitrary temperatures and chemical potentials. This method may significantly increase the computational efficiency of biased grand canonical Monte Carlo simulations, especially for multicomponent mixtures. Although approximate, this approach is amenable to high-throughput and data-intensive investigations where it is preferable to have a large quantity of reasonably accurate simulation data, rather than a smaller amount with a higher accuracy.

Hybrid Monte Carlo with LAMMPS

We describe a strategy for performing canonical and isothermal-isobaric ensemble hybrid Monte Carlo (HMC) simulations with the widely-used Large-scale Atomic/Molecular Massively Parallel Simulator (LAMMPS) molecular dynamics (MD) software package. The overall workflow for the HMC simulations is handled using an external Python driver script, which invokes LAMMPS’ library interface to perform numerically intensive tasks such as MD integration. We document several rigorous consistency checks that have been used to validate our HMC implementation. We also demonstrate that our approach can be readily extended to implement biased HMC sampling schemes for computing free energies. Codes and input files from the documented examples are available on the web.

Quasi-Two-Dimensional Phase Transition of Methane Adsorbed in Cylindrical Silica Mesopores

Using Monte Carlo and molecular dynamics simulations, we examine the adsorption of methane in cylindrical silica mesopores in an effort to understand a possible phase transition of adsorbed methane in MCM-41 and SBA-15 silica that was previously identified by an unexpected increase in the adsorbed fluid density following capillary condensation, as measured by small-angle neutron scattering (SANS) [Chiang, W-S., et al., Langmuir 2016, 32, 8849]. Our initial simulation results identify a roughly 10 % increase in the density of the liquidlike adsorbed phase for either an isotherm with increasing pressure or an isobar with decreasing temperature and that this densification is associated with a local maximum in the isosteric enthalpy of adsorption. Subsequent analysis of the simulated fluid, via computation of bond-orientational order parameters of specific annular layers of the adsorbed fluid, showed that the layers undergo an ordering transition from a disordered, amorphous state to one with two-dimensional hexagonal structure. Furthermore, this two-dimensional restructuring of the fluid occurs at the same thermodynamic state points as the aforementioned densification and local maximum in the isosteric enthalpy of adsorption. We thus conclude that the densification of the fluid is the result of structural reorganization, which is signaled by the maximum in the isosteric enthalpy. Owing to the qualitative similarity of the structural transitions in the simulated and experimental methane fluids, we propose this hexagonal reorganization as a plausible explanation of the densification observed in SANS measurements. Lastly, we speculate how this structural transition may impact the transport properties of the adsorbed fluid.

Efficient Application of Continuous Fractional Component Monte Carlo in the Reaction Ensemble

A new formulation of the Reaction Ensemble Monte Carlo technique (RxMC) combined with the Continuous Fractional Component Monte Carlo method is presented. This method is denoted by serial Rx/CFC. The key ingredient is that fractional molecules of either reactants or reaction products are present and that chemical reactions always involve fractional molecules. Serial Rx/CFC has the following advantages compared to other approaches: (1) One directly obtains chemical potentials of all reactants and reaction products. Obtained chemical potentials can be used directly as an independent check to ensure that chemical equilibrium is achieved. (2) Independent biasing is applied to the fractional molecules of reactants and reaction products. Therefore, the efficiency of the algorithm is significantly increased, compared to the other approaches. (3) Changes in the maximum scaling parameter of intermolecular interactions can be chosen differently for reactants and reaction products. (4) The number of fractional molecules is reduced. As a proof of principle, our method is tested for Lennard-Jones systems at various pressures and for various chemical reactions. Excellent agreement was found both for average densities and equilibrium mixture compositions computed using serial Rx/CFC, RxMC/CFCMC previously introduced by Rosch and Maginn (Journal of Chemical Theory and Computation, 2011, 7, 269-279), and the conventional RxMC approach. The serial Rx/CFC approach is also tested for the reaction of ammonia synthesis at various temperatures and pressures. Excellent agreement was found between results obtained from serial Rx/CFC, experimental results from literature, and thermodynamic modeling using the Peng-Robinson equation of state. The efficiency of reaction trial moves is improved by a factor of 2 to 3 (depending on the system) compared to the RxMC/CFCMC formulation by Rosch and Maginn.

Relationship between Pore-size Distribution and Flexibility of Adsorbent Materials: Statistical Mechanics and Future Material Characterization Techniques

Measurement of the pore-size distribution (PSD) via gas adsorption and the so-called "kernel method" is a widely used characterization technique for rigid adsorbents. Yet, standard techniques and analytical equipment are not appropriate to characterize the emerging class of flexible adsorbents that deform in response to the stress imparted by an adsorbate gas, as the PSD is a characteristic of the material that varies with the gas pressure and any other external stresses. Here, we derive the PSD for a flexible adsorbent using statistical mechanics in the osmotic ensemble to draw analogy to the kernel method for rigid materials. The resultant PSD is a function of the ensemble constraints including all imposed stresses and, most importantly, the deformation free energy of the adsorbent material. Consequently, a pressure-dependent PSD is a descriptor of the deformation characteristics of an adsorbent and may be the basis of future material characterization techniques. We discuss how, given a technique for resolving pressure-dependent PSDs, the present statistical mechanical theory could enable a new generation of analytical tools that measure and characterize certain intrinsic material properties of flexible adsorbents via otherwise simple adsorption experiments.

Classical and quantum many-body effects on the critical properties and thermodynamic regularities of silicon

Using molecular simulation, we determine the critical properties of Si as well as the loci for several remarkable thermodynamic contours spanning the supercritical region of the phase diagram. We consider a classical three-body potential as well as a quantum (tight-binding) many-body model, and determine the loci for the ideality contours, including the Zeno line and the H line of ideal enthalpy. The two strategies (classical or quantum) lead to strongly asymmetric binodals and to critical properties in good agreement with each other. The Zeno and H lines are found to remain linear over a wide temperature interval, despite the changes in electronic structure undergone by the fluid along these contours. We also show that the classical and quantum model yield markedly different results for the parameters defining the H line, the exponents for the power-laws underlying the line of minima for the isothermal enthalpy and for the density required to achieve ideal behavior, most notably for the enthalpy.

Depletion-driven crystallization of cubic colloids sedimented on a surface

Cubic colloids, sedimented on a surface and immersed in a solution of depletant molecules, were modeled with a family of shapes which smoothly varies from squares to circles. Using Wang-Landau simulations with expanded ensembles, we observe the formation of rhombic lattices, square lattices, hexagonal lattices, and a fluid phase. This systematic investigation includes locating transitions between all combinations of the three lattice structures upon changing the shape and transitions between the fluid and crystal upon changing the depletant concentration. The rhombic lattice deforms smoothly between square-like and hexagonal-like angles, depending on both the shape and the depletant concentration. Our results on the effect of the depletant concentration, depletant size, and colloid shape to influence the stability of the fluid and the lattice structures may help guide experimental studies with recently synthesized cubic colloids.

Multicomponent adsorption in mesoporous flexible materials with flat-histogram Monte Carlo methods

We demonstrate an extensible flat-histogram Monte Carlo simulation methodology for studying the adsorption of multicomponent fluids in flexible porous solids. This methodology allows us to easily obtain the complete free energy landscape for the confined fluid-solid system in equilibrium with a bulk fluid of any arbitrary composition. We use this approach to study the adsorption of a prototypical coarse-grained binary fluid in "Hookean" solids, where the free energy of the solid may be described as a simple spring. However, our approach is fully extensible to solids with arbitrarily complex free energy profiles. We demonstrate that by tuning the fluid-solid interaction ranges, the inhomogeneous fluid structure inside the pore can give rise to enhanced selective capture of a larger species through cooperative adsorption with a smaller one. The maximum enhancement in selectivity is observed at low to intermediate pressures and is especially pronounced when the larger species is very dilute in the bulk. This suggest a mechanism by which the selective capture of a minor component from a bulk fluid may be enhanced.

Using Grand Canonical Monte Carlo Simulations to Understand the Role of Interfacial Fluctuations on Solvation at the Water-Vapor Interface

The present work investigates the effect of interfacial fluctuations (predominantly capillary wave-like fluctuations) on the solvation free energy (Δμ) of a monatomic solute at the water-vapor interface. We introduce a grand-canonical-ensemble-based simulation approach that quantifies the contribution of interfacial fluctuations to Δμ. This approach is used to understand how the above contribution depends on the strength of dispersive and electrostatic solute-water interactions at the temperature of 400 K. At this temperature, we observe that interfacial fluctuations do play a role in the variation of Δμ with the strength of the electrostatic solute-water interaction. We also use grand canonical simulations to further investigate how interfacial fluctuations affect the propensity of the solute toward the water-vapor interface. To this end, we track a quantity called the interface potential (surface excess free energy) with the number of water molecules. With increasing number of water molecules, the liquid-vapor interface moves across a solute, which is kept at a fixed position in the simulation. Hence, the dependence of the interface potential on the number of waters models the process of moving the solute through the water-vapor interface. We analyze the change of the interface potential with the number of water molecules to explain that solute-induced changes in the interfacial fluctuations, like the pinning of capillary-wave-like undulations, do not play any role in the propensity of solutes toward water-vapor interfaces. The above analysis also shows that the dampening of interfacial fluctuations accompanies the adsorption of any solute at the liquid-vapor interface, irrespective of the chemical nature of the solute and solvent. However, such a correlation does not imply that dampening of fluctuations causes adsorption.

Evaluation of the grand-canonical partition function using expanded Wang-Landau simulations. IV. Performance of many-body force fields and tight-binding schemes for the fluid phases of silicon

We extend Expanded Wang-Landau (EWL) simulations beyond classical systems and develop the EWL method for systems modeled with a tight-binding Hamiltonian. We then apply the method to determine the partition function and thus all thermodynamic properties, including the Gibbs free energy and entropy, of the fluid phases of Si. We compare the results from quantum many-body (QMB) tight binding models, which explicitly calculate the overlap between the atomic orbitals of neighboring atoms, to those obtained with classical many-body (CMB) force fields, which allow to recover the tetrahedral organization in condensed phases of Si through, e.g., a repulsive 3-body term that favors the ideal tetrahedral angle. Along the vapor-liquid coexistence, between 3000 K and 6000 K, the densities for the two coexisting phases are found to vary significantly (by 5 orders of magnitude for the vapor and by up to 25% for the liquid) and to provide a stringent test of the models. Transitions from vapor to liquid are predicted to occur for chemical potentials that are 10%-15% higher for CMB models than for QMB models, and a ranking of the force fields is provided by comparing the predictions for the vapor pressure to the experimental data. QMB models also reveal the formation of a gap in the electronic density of states of the coexisting liquid at high temperatures. Subjecting Si to a nanoscopic confinement has a dramatic effect on the phase diagram with, e.g. at 6000 K, a decrease in liquid densities by about 50% for both CMB and QMB models and an increase in vapor densities between 90% (CMB) and 170% (QMB). The results presented here provide a full picture of the impact of the strategy (CMB or QMB) chosen to model many-body effects on the thermodynamic properties of the fluid phases of Si.

Impact of Friedel oscillations on vapor-liquid equilibria and supercritical properties in two and three dimensions

We determine the impact of the Friedel oscillations on the phase behavior, critical properties, and thermodynamic contours in films [two dimensions (2D)] and bulk phases [three dimensions (3D)]. Using expanded Wang-Landau simulations, we calculate the grand-canonical partition function and, in turn, the thermodynamic properties of systems modeled with a linear combination of the Lennard-Jones and Dzugutov potentials, weighted by a parameter X (0<X<1). Varying X allows us to control the height of the first Friedel oscillation and to provide a complete characterization of the effect of the metal-like character in the potential on the thermodynamic properties over a wide range of conditions. For 3D systems, we are able to show that the critical parameters exhibit a linear dependence on X and that the loci for the thermodynamic state points, for which the system shows the same compressibility factor or enthalpy as an ideal gas, are two straight lines spanning the subcritical and supercritical regions of the phase diagram for all X values. Reducing the dimensionality to 2D results in a loss of impact of the Friedel oscillation on the critical properties, as evidenced by the virtually constant critical density across the range of X values. Furthermore, our results establish that the straightness of the two ideality lines is retained in 2D and is independent from the height of the first Friedel oscillation in the potential.

Scaling Laws and Critical Properties for fcc and hcp Metals

The determination of the critical parameters of metals has remained particularly challenging both experimentally, because of the very large temperatures involved, and theoretically, because of the many-body interactions that take place in metals. Moreover, experiments have shown that these systems exhibit an unusually strong asymmetry of their binodal. Recent theoretical work has led to new similarity laws, based on the calculation of the Zeno line and of the underlying Boyle parameters, which provided results for the critical properties of atomic and molecular systems in excellent agreement with experiments. Using the recently developed expanded Wang-Landau (EWL) simulation method, we evaluate the grand-canonical partition function, over a wide range of conditions, for 11 fcc and hcp metals (Ag, Al, Au, Be, Cu, Ir, Ni, Pb, Pd, Pt, and Rh), modeled with a many-body interaction potential. This allows us to calculate the binodal, Zeno line, and Boyle parameters and, in turn, obtain the critical properties for these systems. We also propose two scaling laws for the enthalpy and entropy of vaporization, and identify critical exponents of 0.4 and 1.22 for these two laws, respectively.

Self-assembly of trimer colloids: effect of shape and interaction range

Trimers with one attractive bead and two repulsive beads, similar to recently synthesized trimer patchy colloids, were simulated with flat-histogram Monte Carlo methods to obtain the stable self-assembled structures for different shapes and interaction potentials. Extended corresponding states principle was successfully applied to self-assembling systems in order to approximately collapse the results for models with the same shape, but different interaction range. This helps us directly compare simulation results with previous experiment, and good agreement was found between the two. In addition, a variety of self-assembled structures were observed by varying the trimer geometry, including spherical clusters, elongated clusters, monolayers, and spherical shells. In conclusion, our results help to compare simulations and experiments, via extended corresponding states, and we predict the formation of self-assembled structures for trimer shapes that have not been experimentally synthesized.

Relation between pore size and the compressibility of a confined fluid

When a fluid is confined to a nanopore, its thermodynamic properties differ from the properties of a bulk fluid, so measuring such properties of the confined fluid can provide information about the pore sizes. Here, we report a simple relation between the pore size and isothermal compressibility of argon confined in such pores. Compressibility is calculated from the fluctuations of the number of particles in the grand canonical ensemble using two different simulation techniques: conventional grand-canonical Monte Carlo and grand-canonical ensemble transition-matrix Monte Carlo. Our results provide a theoretical framework for extracting the information on the pore sizes of fluid-saturated samples by measuring the compressibility from ultrasonic experiments.

Computational study of trimer self-assembly and fluid phase behavior

The fluid phase diagram of trimer particles composed of one central attractive bead and two repulsive beads was determined as a function of simple geometric parameters using flat-histogram Monte Carlo methods. A variety of self-assembled structures were obtained including spherical micelle-like clusters, elongated clusters, and densely packed cylinders, depending on both the state conditions and shape of the trimer. Advanced simulation techniques were employed to determine transitions between self-assembled structures and macroscopic phases using thermodynamic and structural definitions. Simple changes in particle geometry yield dramatic changes in phase behavior, ranging from macroscopic fluid phase separation to molecular-scale self-assembly. In special cases, both self-assembled, elongated clusters and bulk fluid phase separation occur simultaneously. Our work suggests that tuning particle shape and interactions can yield superstructures with controlled architecture.

Optimizing Noble Gas-Water Interactions via Monte Carlo Simulations

In this work we present optimized noble gas-water Lennard-Jones 6-12 pair potentials for each noble gas. Given the significantly different atomic nature of water and the noble gases, the standard Lorentz-Berthelot mixing rules produce inaccurate unlike molecular interactions between these two species. Consequently, we find simulated Henry's coefficients deviate significantly from their experimental counterparts for the investigated thermodynamic range (293-353 K at 1 and 10 atm), due to a poor unlike potential well term (εij). Where εij is too high or low, so too is the strength of the resultant noble gas-water interaction. This observed inadequacy in using the Lorentz-Berthelot mixing rules is countered in this work by scaling εij for helium, neon, argon, and krypton by factors of 0.91, 0.8, 1.1, and 1.05, respectively, to reach a much improved agreement with experimental Henry's coefficients. Due to the highly sensitive nature of the xenon εij term, coupled with the reasonable agreement of the initial values, no scaling factor is applied for this noble gas. These resulting optimized pair potentials also accurately predict partitioning within a CO2-H2O binary phase system as well as diffusion coefficients in ambient water. This further supports the quality of these interaction potentials. Consequently, they can now form a well-grounded basis for the future molecular modeling of multiphase geological systems.