Generalized geometric cluster algorithm for fluid simulation

Evaluation of Osmotic Virial Coefficients via Restricted Gibbs Ensemble Simulations, with Support from Gas-Phase Mixture Coefficients

We present a method for computing osmotic virial coefficients in explicit solvent via simulation in a restricted Gibbs ensemble. Two equivalent phases are simulated at once, each in a separate box at constant volume and temperature and each in equilibrium with a solvent reservoir. For osmotic coefficient BN, a total of N solutes are individually exchanged back and forth between the boxes, and the average distribution of solute numbers between the boxes provides the key information needed to compute BN. Separately, expressions are developed for BN as a series in solvent reservoir density ρ₁, with the coefficients of the series expressed in terms of the usual gas-phase mixture coefficients B_ij. Normally, the B_ij are defined for an infinite volume, but we suggest that the observed dependence of B_ij on system size L can be used to estimate L dependence of the BN, allowing them to be computed accurately at L → ∞ while simulating much smaller system sizes than otherwise possible. The methods for N = 2 and 3 are demonstrated for two-component mixtures of size-asymmetric additive hard spheres. The proposed methods are demonstrated to have greater precision than established techniques, for a given amount of computational effort. The ρ₁ series for BN when applied by itself is (for this noncondensing model) found to be the most efficient in computing accurate osmotic coefficients for the solvent densities considered here.

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).

Efficient Irreversible Monte Carlo Samplers

We present here two irreversible Markov chain Monte Carlo algorithms for general discrete state systems. One of the algorithms is based on the random-scan Gibbs sampler for discrete states and the other on its improved version, the Metropolized-Gibbs sampler. The algorithms we present incorporate the lifting framework with skewed detailed balance condition and construct irreversible Markov chains that satisfy the balance condition. We have applied our algorithms to 1D 4-state Potts model. The integrated autocorrelation times for magnetization and energy density indicate a reduction of the dynamical scaling exponent from z ≈ 1 to z ≈ 1/2. In addition, we have generalized an irreversible Metropolis-Hastings algorithm with skewed detailed balance, initially introduced by Turitsyn et al. [ Physica D 2011, 240, 410] for the mean field Ising model, to be now readily applicable to classical spin systems in general; application to 1D 4-state Potts model indicate a square root reduction of the mixing time at high temperatures.

Early stages of phase selection in MOF formation observed in molecular Monte Carlo simulations

Metal–organic frameworks (MOF) comprising metal nodes bridged by organic linkers show great promise because of their guest-specific gas sorption, separation, drug-delivery, and catalytic properties. The selection of metal node, organic linker, and synthesis conditions in principle offers engineered control over both structure and function. For MOFs to realise their potential and to become more than just promising materials, a degree of predictability in the synthesis and a better understanding of the self-assembly or initial growth processes is of paramount importance. Using cobalt succinate, a MOF that exhibits a variety of phases depending on synthesis temperature and ligand to metal ratio, as proof of concept, we present a molecular Monte Carlo approach that allows us to simulate the early stage of MOF assembly. We introduce a new Contact Cluster Monte Carlo (CCMC) algorithm which uses a system of overlapping “virtual sites” to represent the coordination environment of the cobalt and both metal–metal and metal–ligand associations. Our simulations capture the experimentally observed synthesis phase distinction in cobalt succinate at 348 K. To the best of our knowledge this is the first case in which the formation of different MOF phases as a function of composition is captured by unbiased molecular simulations. The CCMC algorithm is equally applicable to any system in which short-range attractive interactions are a dominant feature, including hydrogen-bonding networks, metal–ligand coordination networks, or the assembly of particles with “sticky” patches, such as colloidal systems or the formation of protein complexes.

Clusters produced in simulations of MOF assembly at metal : ligand ratios of (a) 5 : 1 and (b) 1 : 1 with identical interaction potentials.

Designing Molecular Building Blocks for the Self-assembly of Complex Porous Networks

The creation of molecular or colloidal building blocks which can self-assemble into complex, ordered porous structures has been long sought-after, and so are the guiding principles behind this creation. The pursuit of this goal has led to the creation of novel classes of materials like metal organic frameworks (MOFs) and covalent organic frameworks (COFs). In theory, a tremendous number of structures can be formed by these materials due to the variety of geometries available to their building blocks. However, most realized crystal structures tend to be simple or homoporous and typically assemble from building blocks with high degrees of symmetry. Building blocks with low degrees of symmetry suitable for assembly into the more complex structures tend to assemble into polymorphous or disordered structures instead. In this work, we use Monte Carlo simulations of patchy vertex-like building blocks to show how the addition of chemical specificity via orthogonally reacting functional sites can allow vertex-like building blocks with even asymmetric geometries to self-assemble into ordered crystallites of various complex structures. In addition to demonstrating the utility of such a strategy in creating ordered, heteroporous structures, we also demonstrate that it can be used as a means for tuning specific features of the crystal structure, accomplishing such aims as the control of relative pore sizes. We also discuss heuristics for properly designing molecules so that they can assemble into target structures.

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.

Assembly of multi-flavored two-dimensional colloidal crystals

We systematically investigate the assembly of binary multi-flavored colloidal mixtures in two dimensions. In these mixtures all pairwise interactions between species may be tuned independently. This introduces an additional degree of freedom over more traditional binary mixtures with fixed mixing rules, which is anticipated to open new avenues for directed self-assembly. At present, colloidal self-assembly into non-trivial lattices tends to require either high pressures for isotropically interacting particles, or the introduction of directionally anisotropic interactions. Here we demonstrate tunable assembly into a plethora of structures which requires neither of these conditions. We develop a minimal model that defines a three-dimensional phase space containing one dimension for each pairwise interaction, then employ various computational techniques to map out regions of this phase space in which the system self-assembles into these different morphologies. We then present a mean-field model that is capable of reproducing these results for size-symmetric mixtures, which reveals how to target different structures by tuning pairwise interactions, solution stoichiometry, or both. Concerning particle size asymmetry, we find that domains in this model's phase space, corresponding to different morphologies, tend to undergo a continuous "rotation" whose magnitude is proportional to the size asymmetry. Such continuity enables one to estimate the relative stability of different lattices for arbitrary size asymmetries. Owing to its simplicity and accuracy, we expect this model to serve as a valuable design tool for engineering binary colloidal crystals from multi-flavored components.

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.

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.

Collective translational and rotational Monte Carlo cluster move for general pairwise interaction

Virtual move Monte Carlo is a cluster algorithm which was originally developed for strongly attractive colloidal, molecular, or atomistic systems in order to both approximate the collective dynamics and avoid sampling of unphysical kinetic traps. In this paper, we present the algorithm in the form, which selects the moving cluster through a wider class of virtual states and which is applicable to general pairwise interactions, including hard-core repulsion. The newly proposed way of selecting the cluster increases the acceptance probability by up to several orders of magnitude, especially for rotational moves. The results have their applications in simulations of systems interacting via anisotropic potentials both to enhance the sampling of the phase space and to approximate the dynamics.

Collective translational and rotational Monte Carlo moves for attractive particles

Virtual move Monte Carlo is a Monte Carlo (MC) cluster algorithm forming clusters via local energy gradients and approximating the collective kinetic or dynamic motion of attractive colloidal particles. We carefully describe, analyze, and test the algorithm. To formally validate the algorithm through highlighting its symmetries, we present alternative and compact ways of selecting and accepting clusters which illustrate the formal use of abstract concepts in the design of biased MC techniques: the superdetailed balance and the early rejection scheme. A brief and comprehensive summary of the algorithms is presented, which makes them accessible without needing to understand the details of the derivation.

Nanoparticle-controlled aggregation of colloidal tetrapods

Tetrapods are among the most promising building blocks for nanoscale self-assembly, offering various desirable features. Whereas these particles can be fabricated with remarkable precision, comparatively less is known about their aggregation behavior. Employing a novel, powerful simulation method, we demonstrate that charged nanoparticles offer considerable control over the assembly of tip-functionalized tetrapods. Extending these findings to tetrapods confined to a gas/liquid interface, we show that regular structures can be achieved even without functionalization.

Depletion potentials in highly size-asymmetric binary hard-sphere mixtures: comparison of simulation results with theory

We report a detailed study, using state-of-the-art simulation and theoretical methods, of the effective (depletion) potential between a pair of big hard spheres immersed in a reservoir of much smaller hard spheres, the size disparity being measured by the ratio of diameters q ≡ σ(s)/σ(b). Small particles are treated grand canonically, their influence being parameterized in terms of their packing fraction in the reservoir η(s)(r). Two Monte Carlo simulation schemes--the geometrical cluster algorithm, and staged particle insertion--are deployed to obtain accurate depletion potentials for a number of combinations of q ≤ 0.1 and η(s)(r). After applying corrections for simulation finite-size effects, the depletion potentials are compared with the prediction of new density functional theory (DFT) calculations based on the insertion trick using the Rosenfeld functional and several subsequent modifications. While agreement between the DFT and simulation is generally good, significant discrepancies are evident at the largest reservoir packing fraction accessible to our simulation methods, namely, η(s)(r) = 0.35. These discrepancies are, however, small compared to those between simulation and the much poorer predictions of the Derjaguin approximation at this η(s)(r). The recently proposed morphometric approximation performs better than Derjaguin but is somewhat poorer than DFT for the size ratios and small-sphere packing fractions that we consider. The effective potentials from simulation, DFT, and the morphometric approximation were used to compute the second virial coefficient B(2) as a function of η(s)(r). Comparison of the results enables an assessment of the extent to which DFT can be expected to correctly predict the propensity toward fluid-fluid phase separation in additive binary hard-sphere mixtures with q ≤ 0.1. In all, the new simulation results provide a fully quantitative benchmark for assessing the relative accuracy of theoretical approaches for calculating depletion potentials in highly size-asymmetric mixtures.

A cluster algorithm for Monte Carlo simulation at constant pressure

We propose an efficient algorithm to sample the volume in Monte Carlo simulations in the isobaric-isothermal ensemble. The method is designed to be applied in the simulation of hard-core models at high density. The algorithm is based in the generation of clusters of particles. At the volume change step, the distances between pairs of particles belonging to the same cluster do not change. This is done by rescaling the positions of the center of mass of each cluster instead of the position of each individual particle. We have tested the performance of the algorithm by simulating fluid and solid phases of hard spheres, finding that in both cases the algorithm is much more efficient than the standard procedure. Moreover, the efficiency of the method measured in terms of correlation "time" does not depend on the system size in contrast with the standard method, in which the sampling becomes rapidly inefficient as the system size increases. We have used the procedure to compute with high precision the equation of state of the face-centered-cubic phase of the hard sphere system for different system sizes. Using these results we have estimated the equation of state at the thermodynamic limit. The results are compared to different equations of state proposed in literature.

Methods for Monte Carlo simulations of biomacromolecules

The state-of-the-art for Monte Carlo (MC) simulations of biomacromolecules is reviewed. Available methodologies for sampling conformational equilibria and associations of biomacromolecules in the canonical ensemble, given a continuum description of the solvent environment, are reviewed. Detailed sections are provided dealing with the choice of degrees of freedom, the efficiencies of MC algorithms and algorithmic peculiarities, as well as the optimization of simple movesets. The issue of introducing correlations into elementary MC moves, and the applicability of such methods to simulations of biomacromolecules is discussed. A brief discussion of multicanonical methods and an overview of recent simulation work highlighting the potential of MC methods are also provided. It is argued that MC simulations, while underutilized biomacromolecular simulation community, hold promise for simulations of complex systems and phenomena that span multiple length scales, especially when used in conjunction with implicit solvation models or other coarse graining strategies.

Monte Carlo study of the triangular lattice gas with first- and second-neighbor exclusions

We formulate a Swendsen-Wang-like version of the geometric cluster algorithm. As an application, we study the hard-core lattice gas on the triangular lattice with the first- and second-neighbor exclusions. The data were first analyzed by finite-size scaling without including logarithmic corrections. We determine the critical chemical potential as mu_{c}=1.75682(2) and the critical particle density as rho_{c}=0.180(4) . From the Binder ratio Q and susceptibility chi , the thermal and magnetic exponents are estimated as y_{t}=1.51(1) approximately 32 and y_{h}=1.8748(8) approximately 158 , respectively, while the analyses of energylike quantities yield y_{t} ranging from 1.440(5) to 1.470(5). Nevertheless, the data for energylike quantities are also well described by theoretically predicted scaling formulas with logarithmic corrections and with exponent y_{t}=32 . These results are very similar to the earlier study for the four-state Potts model on the square lattice [J. Stat. Phys. 88, 567 (1997)], and strongly support the general belief that the model is in the four-state Potts universality class. The dynamic scaling behavior of the Metropolis simulation and the combined method of the Metropolis and the geometric cluster algorithm is also studied; the former has a dynamic exponent z_{int} approximately 2.21 and the latter has z_{int} approximately 1.60 .

Avoiding unphysical kinetic traps in Monte Carlo simulations of strongly attractive particles

We introduce a "virtual-move" Monte Carlo algorithm for systems of pairwise-interacting particles. This algorithm facilitates the simulation of particles possessing attractions of short range and arbitrary strength and geometry, an important realization being self-assembling particles endowed with strong, short-ranged, and angularly specific ("patchy") attractions. Standard Monte Carlo techniques employ sequential updates of particles and can suffer from low acceptance rates when attractions are strong. In this event, collective motion can be strongly suppressed. Our algorithm avoids this problem by proposing simultaneous moves of collections (clusters) of particles according to gradients of interaction energies. One particle first executes a "virtual" trial move. We determine which of its neighbors move in a similar fashion by calculating individual bond energies before and after the proposed move. We iterate this procedure and update simultaneously the positions of all affected particles. Particles move according to an approximation of realistic dynamics without requiring the explicit computation of forces and without the step size restrictions required when integrating equations of motion. We employ a size- and shape-dependent damping of cluster movements, motivated by collective hydrodynamic effects neglected in simple implementations of Brownian dynamics. We discuss the virtual-move algorithm in the context of other Monte Carlo cluster-move schemes and demonstrate its utility by applying it to a model of biological self-assembly.

Simulation of phase transitions in highly asymmetric fluid mixtures

We present a novel method for the accurate numerical determination of the phase behavior of fluid mixtures having large particle-size asymmetries. By incorporating the recently developed geometric cluster algorithm within a restricted Gibbs ensemble, we are able to probe directly the density and concentration fluctuations that drive phase transitions, but that are inaccessible to conventional simulation algorithms. We develop a finite-size scaling theory that relates these density fluctuations to those of the grand-canonical ensemble, thereby enabling accurate location of critical points and coexistence curves of multicomponent fluids. Several illustrative examples are presented.

Effective interactions in mixtures of silica microspheres and polystyrene nanoparticles

We investigate the effect of small concentrations of highly charged nanoparticles on the stability of uncharged colloidal microspheres using large-scale simulations. Employing pair potentials that accurately represent mixtures of silica microspheres and polystyrene nanoparticles as studied experimentally, we are able to demonstrate that nanoparticle-induced stabilization can arise from a relatively weak van der Waals attraction between the colloids and nanoparticles. This demonstrates that the nanoparticle haloing mechanism for colloidal stabilization is of considerable generality and potentially can be applied to large classes of systems. The range of optimal nanoparticle concentrations can be tuned by controlling the attraction between colloids and nanoparticles.

Phase transitions in highly asymmetric binary hard-sphere fluids: Fluid-fluid binodal from a two-component mixture theory

Fluid-fluid binodals of binary hard-sphere mixtures are computed from the recently proposed fundamental measure functional-mean spherical approximation closure of the two-component Ornstein-Zernike equation. The results, especially in the dense fluid region that was not accessible by previous theoretical methods, are compared with the corresponding ones for the one-component fluid of big spheres with effective potential obtained from the same closure. The general trends are those expected for hard-sphere potentials but small difference are detectable. The overall agreement found validates the equivalence of the two descriptions for size ratios R = 8.5 or greater.

Influence of polydispersity on the critical parameters of an effective-potential model for asymmetric hard-sphere mixtures

We report a Monte Carlo simulation study of the properties of highly asymmetric binary hard-sphere mixtures. This system is treated within an effective fluid approximation in which the large particles interact through a depletion potential [R. Roth, Phys. Rev. E 62 5360 (2000)] designed to capture the effects of a virtual sea of small particles. We generalize this depletion potential to include the effects of explicit size dispersity in the large particles and consider the case in which the particle diameters are distributed according to a Schulz form having a degree of polydispersity 14%. The resulting alteration (with respect to the monodisperse limit) of the metastable fluid-fluid critical point parameters is determined for two values of the ratio of the diameters of the small and large particles: q(triple bond)sigma(s)/(-)sigma(b)=0.1 and q=0.05. We find that the inclusion of polydispersity moves the critical point to lower reservoir volume fractions of the small particles and high volume fractions of the large ones. The estimated critical point parameters are found to be in good agreement with those predicted by a generalized corresponding states argument which provides a link to the known critical adhesion parameter of the adhesive hard-sphere model. Finite-size scaling estimates of the cluster percolation line in the one phase fluid region indicate that inclusion of polydispersity moves the critical point deeper into the percolating regime. This suggests that phase separation is more likely to be preempted by dynamical arrest in polydisperse systems.

Colloidal stabilization via nanoparticle halo formation

We present a detailed numerical study of effective interactions between micrometer-sized silica spheres, induced by highly charged zirconia nanoparticles. It is demonstrated that the effective interactions are consistent with a recently discovered mechanism for colloidal stabilization. In accordance with the experimental observations, small nanoparticle concentrations induce an effective repulsion that counteracts the intrinsic van der Waals attraction between the colloids and thus stabilizes the suspension. At higher nanoparticle concentrations an attractive potential is recovered, resulting in reentrant gelation. Monte Carlo simulations of this highly size-asymmetric mixture are made possible by means of a geometric cluster Monte Carlo algorithm. A comparison is made to results obtained from the Ornstein-Zernike equations with the hypernetted-chain closure.

Structure of highly asymmetric hard-sphere mixtures: An efficient closure of the Ornstein-Zernike equations

A simple modification of the reference hypernetted chain (RHNC) closure of the multicomponent Ornstein-Zernike equations with bridge functions taken from Rosenfeld's hard-sphere bridge functional is proposed. Its main effect is to remedy the major limitation of the RHNC closure in the case of highly asymmetric mixtures--the wide domain of packing fractions in which it has no solution. The modified closure is also much faster, while being of similar complexity. This is achieved with a limited loss of accuracy, mainly for the contact value of the big sphere correlation functions. Comparison with simulation shows that inside the RHNC no-solution domain, it provides a good description of the structure, while being clearly superior to all the other closures used so far to study highly asymmetric mixtures. The generic nature of this closure and its good accuracy combined with a reduced no-solution domain open up the possibility to study the phase diagram of complex fluids beyond the hard-sphere model.