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

Simulations in the grand canonical (GC) ensemble, where particles are exchanged between a virtual reservoir and positions throughout the simulated system, are useful in a variety of contexts. By establishing equilibrium between each position within the simulation box and a virtual reservoir of a known chemical potential μ, GC algorithms enable different regions of the simulation box to reach an equilibrated distribution of particles without relying on transport across the box. In cases where there are barriers to transport, the use of insertion and removal moves (or coupled exchange moves, as in the isomolar semigrand canonical ensemble¹) can speed up this equilibration among different microenvironments.² Furthermore, using the same virtual reservoir for different simulation environments—for instance, different phases or the interior of a porous solid versus its surrounding bulk medium³—enables the study of equilibrium partitioning of particles among these environments without incorporating multiple environments within the same simulation box. This may be more convenient and efficient for several possible reasons—for instance, to avoid unfavorable scaling of simulation cost with simulation size or to avoid structural perturbations and kinetic barriers arising at interfaces between environments.

Simulation methods that employ particle insertion and removal moves are often deemed impractical to perform in condensed phases. Within conventional Grand Canonical Monte Carlo (GCMC), the acceptance probability for adding a particle to a randomly selected position within a simulation box is given by min [V N/(N + 1)exp(−βΔU + β μ), 1], where V is the system volume, N the number of particles, β = 1/k_BT, and ΔU is the potential energy of interaction of the inserted particle with the system.⁴ Any positions that overlap existing particles will generally have high ΔU, so the acceptance probability is closely related to the probability of finding a cavity in the box large enough to accommodate the addition of a new particle. In a dense phase, this probability is typically small; it will decrease exponentially with the size of the particle to be added, so even optimizing the insertion success by actively searching for voids⁵ will eventually become unproductive for large particles.

One class of approaches to overcome this challenge uses gradual coupling and decoupling of the particle-particle interactions, either within a set of exchangeable systems via the simulated tempering approach⁶ or in multi-step growth and removal of a single “ghost” particle at a time.⁷ Another logical solution to this problem is to take away one or more existing particle to make room for the added one. In the isomolar semi-grand canonical ensemble techniques, which are extremely useful if the particle sizes and properties are not too dissimilar, each added particle is swapped for a single particle of a different type.^1,2

Configuration bias Monte Carlo⁸ (CBMC) is a proven method for sampling new configurations within crowded environments, most often applied to simulations of segmented chain structures with many internal degrees of freedom. The idea of this method is that a number k of test insertions are made for each segment of the structure as it is grown, allowing the selection of sites that are lower in potential energy and less likely to be rejected. In the context of a move that exchanges a large solute for multiple small solvent particles, the sequential insertion of individual solvent molecules into the cavity, as illustrated in Figure 1, takes the place of the insertion of segments. Variations of this approach have been used in canonical⁹ and Gibbs ensemble¹⁰ MC simulations of mixtures of large and small particles. The number of small particles to be exchanged for a large particle within a GCMC simulation is not predetermined, so we need to define a chemical potential μ_S for the small particles as well as the large particles to calculate the relative probabilities of configurations with different numbers of small particles. In exchanging one solute for one or more solvent molecules, one then can select the number of particles to be added or removed by considering the μ_S-dependent probability weights of each possible number of particles to exchange, as demonstrated by Lü and Kindt in GCMC simulations of self-assembled polymer chains.¹¹ This new approach is called solvent repacking Monte Carlo (SRMC).

As SRMC allows the generation of new solvent configurations with arbitrary numbers of particles within a defined region, it can be used for GCMC simulation even in the pure solvent. In solvent phases where the formation of cavities large enough to accommodate an additional solvent molecule is thermodynamically unfavorable, making the acceptance probability of conventional GCMC very low, solvent repacking offers the possibility of interchanging states with different numbers of solvents, bypassing an unfavorable intermediate state with a molecule-sized cavity.

Efficient application of CBMC as described above will become impossible when the number and density of solvent particles to repack in the solute’s cavity become high enough. Within a large cavity, the distribution of the first few trial positions generated will sample the Boltzmann distribution for a small number of solutes interacting with the cavity. Depending on the cavity size and the nature of the inter-particle interactions, they will either adopt a gas-like distribution in the cavity interior or a distribution characteristic of the vapor/solvent interface. These distributions, even when sampled perfectly (i.e., in the limit of infinite k), will in general differ from the distribution of solvent configurations in the bulk. Thus, particles that are initially placed poorly will preclude solvent positions generated later in the algorithm from packing efficiently in bulk-like arrangements, contributing to poor acceptance probabilities.

The use of a biasing function that departs from the Boltzmann weighting could in principle improve the probability that the CBMC algorithm will generate bulk-like packing arrangements within a cavity, if that biasing function favors the selection of positions for solvent insertion that will better accommodate the solvents that are not yet inserted. The use of auxiliary functions that differ from the actual potential energy function to generate configurations within CBMC is well established.¹² For the final distribution to reflect the landscape of the actual function rather than the auxiliary function, a correction that depends on the difference of the two functions as calculated only for the final result must be included in the acceptance probability. In previous applications, an auxiliary function was used because it was less computationally expensive than the true potential; the full potential function then needs only to be calculated once and not k times. In the current approach, the aim instead is to use auxiliary functions to generate better configurations and thus to improve acceptance probabilities. A natural choice for this auxiliary biasing function is the radial distribution function g(r) of the pure solvent, whose natural logarithm gives the potential of mean force between solvent particles in units of k_BT.¹³ For the hard disk systems explored in this study, the structure in this function is entirely the reflection of the degree to which different interparticle distances permit efficient packing of the system as a whole, since the true potential energy function is independent of r beyond its value at contact.

This report will contain the derivations of acceptance probabilities for solvent repacking MC moves, along with details of implementing these methods. Section II will be followed by Sec. III exploring the efficiency and accuracy of SRMC in the context of two-dimensional monodisperse and bidisperse hard disk systems. Hard disk systems are used here for the simplicity of their interparticle potential, which highlights the importance of the auxiliary potential for selecting particle positions, and because they are widely studied model systems whose thermodynamic, structural, and dynamic properties remain topics of active research.^14–20 The capabilities of SRMC are demonstrated by showing that it can reproduce the features of the hard disk ordering transition—a first-order transition between fluid and hexatic phases followed by a second-order hexatic to solid transition—as recently determined by canonical ensemble simulations.^16,17 Its effectiveness for mixtures is shown in simulations of binary mixtures of disks of different sizes, both in the fluid and ordered phases. The partitioning of large disks between ordered and disordered phases, and the nature of the translational order in the ordered phase, is determined at a single point on the phase boundary for binary mixtures with disk diameter ratios d = 1.4, 2.5, and 3.0.

To derive expressions for the acceptance probabilities of the SRMC trial moves, we will start with moves that maintain N in a pure solvent, proceed to moves that allow fluctuations of N at constant chemical potential in a pure solvent, and conclude with moves that exchange a variable number of solvents for a single solute at fixed chemical potentials of both species.

The first type of move is an attempted rearrangement of the subset of solvent molecules found within a randomly selected subvolume (“cavity”) within the system, keeping the number of molecules constant. This is effectively an application of the CBMC method where order of insertion of the segments is arbitrary; it is not likely to be particularly efficient, but is included for explanatory purposes. The simplest method to select and define a cavity is to select a point in the system at random and identify all solvent molecules located within a radius r_cav of that point. Defining this number of molecules as n_cav, generation of a new configuration for the system consists of n_cav cycles. In each cycle i, a number k of positions are generated that randomly and uniformly sample the volume of the cavity. For each position j generated in cycle i, the potential energy u_i,j is calculated as the sum of the energies of interaction for a particle at that position with all particles outside the cavity and with particles added to the cavity in earlier cycles of the same repacking attempt. After k positions are generated, one of them is chosen to place the ith particle for the current trial move; the probability of choosing a position j′ is proportional to its Boltzmann factor,

In the case of the hard disk system, the Boltzmann factors are either 0 (for a position that overlaps other disks) or 1 (for a position that does not), so all allowable positions are equally probable. If P_i,j = 0 for all trial positions at some i ≤ n_cav, the move fails. Otherwise, a new set of n_cav positions has been generated. To determine the acceptance probability for keeping the new configuration of these particles, the Rosenbluth weight W of the new configuration must first be found

Next, the Rosenbluth weight of the current configuration of particles in the cavity must be calculated. First, k − 1 positions within the empty cavity are generated at random and their potential energies u_i,j determined. Next, one of the existing n_cav particle positions is chosen to be the kth position; its potential energy is calculated, it is added to the cavity, and another set of k − 1 positions are generated. The process is repeated for n_cav cycles, yielding values of u_i,j to calculate W for the current configuration. The probabilities α_R→R′ of generating a trial move between configurations R and R′ and its reverse α_R′→R are unequal, but their ratio can be calculated as

In general, for a biased Monte Carlo method in the canonical ensemble, the acceptance probability must include this ratio as

Plugging in Eq. (1) yields the result familiar from the CBMC method as commonly used

The actual difference between the potential energies of the new and current configurations, U_R′ − U_R, does not enter directly into Eq. (5) because their Boltzmann weights are already built in to the probability of choosing a position in the numerator of Eq. (1).

The method as described is prone to high failure rates in dense systems, where random selection of positions for insertion is unlikely to produce an efficiently packed cavity. The acceptance probability might be improved by using an auxiliary potential U′ instead of the actual potential in selecting one out of the k positions generated

This causes a shift in the probability of a given trial structure being generated, so the acceptance probability must be adapted to compensate

where W′ is defined as in Eq. (2) with the substitution of W′ and $u i , j ′$ for W and u_i,j, as previously derived.¹² 

As the cavity fills up with solvent, the need for an auxiliary potential drops. The choice for the position of the last solvent molecule added to the cavity should optimally be based solely on the true potential, because there is no longer any need to “leave room” for further insertions. The auxiliary potential can be switched over from the effective potential to the true potential as i is increased. Conversely, the number of trial positions k needed to find a good position will tend to increase as the cavity fills up and also may be made dependent on i.

The true utility of applying the configuration bias concept to the solvent repacking problem becomes apparent in grand canonical ensemble simulations, where n_cav is not restricted to have the same value in the initial and final states. The ratio of weights of states with different N is associated with an additional factor of (f/Λ^d)^ΔN, where fugacity f = exp(βμ) and d is the dimension of the system. (Note that we assign a value of 1 to the de Broglie thermal wavelength Λ in all calculations. Throughout the text, chemical potentials or fugacities will be given depending on which is more convenient in context.) If we were to assume an equal probability of move attempts to raise or lower the number of particles in a cavity by 1, we need to incorporate this additional factor into Eq. (4). To use Eq. (4), we further need the ratio of the probability of generating and selecting the current state (with n_cav particles in the cavity) in a removal-type move to the probability of generating the configuration with the additional particle,

The factorial terms account for the number of sequences in which a given set of positions can be selected, and the 1/v_cavⁿ terms account for the probability of selecting a set of n particular positions within the cavity volume v_cav.

Therefore, substituting the definition of the Rosenbluth weight (Eq. (2)) into Eq. (4) yields

The approach of predetermining whether a given solvent repacking move will end with an increased, decreased, or equal number of particles in the cavity in distinct types of move attempt can be improved with a move where all these possibilities are built in as options, using the “polydisperse insertion” strategy developed for GCMC simulations of equilibrium polymers.¹¹ Already in the constant NVT method described above, configurations containing 0, 1, …, n_cav solvents occupying the cavity are generated and can serve as possible final states. In the GCMC method, the procedure for filling the cavity is continued until some predetermined condition is met (e.g., a maximum value for n_cav or a minimum W_new; for hard disks, W_new will eventually reach zero as the cavity runs out of space for more disks), whatever may be the original number of particles in the cavity.

Once n_max new positions for solvents in the cavity have been generated, a choice may be made to use a subset i′ ∈ {0, n_max} of them,

where the weight ω_i of a state containing a specific number i of particles in the cavity incorporates the factors that go along with its weight in the grand canonical ensemble as well as the probability of generating it,

To accurately account for the probability that a given repacking move will be reversed, new particle positions must also be generated for the current configuration of solvent within the cavity until the maximum particle condition is met. The overall acceptance probability of

then makes the method consistent with Eq. (9). If this acceptance criterion is not met, a second possibility to change the state of the system presents itself: the weights ω_i,current for the cavity containing its current configuration with particles added or removed having been calculated anyhow, a number i′ of particles for the cavity to contain is chosen with probability P(i′) given by Eq. (11). The most probable state to choose is typically i′ = n_cav, which results in no change, but a finite chance remains of adding or removing one or more particles.

In a binary mixture of small and large disks, the SRMC algorithm can be used as described above to maintain a constant chemical potential of the small disk component; the presence of a large disk in the cavity does not affect the algorithm except as it contributes to the potential energies u and the auxiliary potential u′ for interaction with the small disks being added or removed. The exchange of a large disk for small disks, maintaining constant chemical potential of the large disk component and allowing local equilibration of large disk concentration without requiring diffusion of large disks, is the last type of move we must consider. For this study, each insertion or removal move involving large disks will only involve a single large disk. Insertion is attempted at a randomly chosen position in the system; insertion attempts where the large disk overlaps another large disk fail immediately. All small disks overlapping the large disk are slated for removal; in other words, the cavity is a circle whose diameter is the average of the large and small disks’ diameters, centered at the point of insertion (see Fig. 1(a)). (In a system where solvent configuration is strongly coupled to the presence of solute, the cavity could be expanded to include a region of strongly perturbed solvent—a solute and its associated solvent shell would then be exchanged for a number of bulk solvent particles.)

The procedure to attempt removal of a large disk is to repack the cavity with solvent, with the number of solvents added chosen just as in the all-solvent move. Therefore, to calculate the acceptance probability of the large disk insertion move, the probability weights ω_i for repacking the cavity up to i = n_max must be calculated. (If the number of solutes to be replaced, n_cav, is less than n_max, then the insertion move automatically fails, because it cannot be reversed by a removal move.) The probability of attempting a large disk insertion at a specific position is simply 1/V, while the probability of replacing that large disk with a particular arrangement of solvent is inversely proportional to the number of large disks present and proportional to factors already discussed in Eqs. (8) and (9),

From this ratio, the acceptance probabilities for addition and removal moves can be obtained as

and

where u_L is the potential energy of interaction between the large disk and the rest of the system and f_L is the large disk fugacity. As in solvent-repacking for the single-component system, a failed attempt to insert a large disk can be followed by an attempt to change the number of small disks within the cavity, using the weights ω_i already calculated.

For simulations of single-component systems of hard disks, solvent repacking move attempts and conventional single-particle translation attempts were made in a 1:4 ratio. For each solvent repacking move, the cavity radius was selected from a uniform distribution over the range [1,1.4] in units of the disk diameter σ, so that n_cav was typically in the range of ∼3-8 disks at densities near the freezing point. In simulations of mixtures, the total ratio of solvent repacking moves to conventional translational moves was 1:2; 50% of solvent repacking moves were attempts to change solvent number and configuration alone, 25% were attempts to add a solute, and 25% were attempts to remove a solute. The biasing function used in all simulations was a pairwise potential derived from the radial distribution function g₀(r) of hard disks at area fraction η = 0.69, slightly below the freezing transition, which is shown in Fig. 2. The auxiliary bias potential was scaled by a factor γ_i for the ith disk to be added,

Values for γ_i and for the number of positions generated k_i are given in Table I; these were chosen after some preliminary tests to optimize acceptance probabilities for the addition and removal of d = 3 large disks.

The typical acceptance probability for trial moves was often very small (<1 × 10⁻¹⁰). A 32-bit algorithm generates a value of exactly 0 within an average of 2³² = 4.3 × 10⁹ function calls and can thus lead to the spurious acceptance of a proposed move during a long simulation. To avoid this problem without increasing the computational cost of random number generation, uniform random numbers were generated by a 32-bit Mersenne Twistern algorithm on the interval [0,1]; those that were lower than 2.33 × 10⁻⁹ (i.e., 10 times 2⁻³²) were discarded and replaced by a newly generated random number multiplied by 2.33 × 10⁻⁹.

Simulations were parallelized using a domain-decomposition scheme. The system (area = 350 σ × 350 σ) was divided into an 8 × 8 or a 12 × 12 grid with a randomly selected origin. Simple translational move attempts and solvent repacking attempts were performed within each grid square on separate processors. Due to the short-range nature of the hard disk potential, changes to the structure in a grid do not affect the sampling of configurations in a neighboring grid except within a border of width σ. To reduce the necessity for communication among processors, all disk positions were fixed within these borders; likewise insertions, removals, and translational moves from the grid square interior to the border were prohibited. After a predetermined number of solvent repacking move attempts (500-1000) and simple translational move attempts (1000-2000) were performed within each grid square, a new origin for the grid was chosen and disks were redistributed among the processors. The introduction of large disks necessitated the use of correspondingly wider borders. Note that the number of large disks within the local grid square and the area of the grid square can be used for N_L and V where they appear in Eqs. (13)–(15).

After each set of simple translational moves, histograms of pair distances for like and unlike types of disk were updated. Pressure (in units of k_BT/σ²) is calculated as

where the values for the radial distribution functions g_SS, etc., at contact were obtained by fitting the pair distance histograms to a third-order polynomial out to 1.05 times the contact length.

The complex hexagonal order parameter S₆ was calculated for each small disk by searching for all small-disk neighbors within a distance of 1.5 σ and taking the average of exp(6iθ) over the neighbors, where θ is the angle formed by the difference vector between the disk and its neighbor and the x axis. The average over all small-disk neighbors was used, even for particles with fewer than 6 neighbors, so that small disks that neighbor large disks would be treated on an even footing with those with all small-disk neighbors.

A histogram of displacement vectors between small disks that lie within ±0.1 σ of a multiple of a unit vector u was constructed for a single configuration and divided by $( ρ S 2 L 2 / a bin )$⁠, where a_bin = 0.2 σ × 0.02 σ is the area of the bin (projected distances along u were binned with a resolution of 0.02 σ) to yield g(Δy). For the structures analyzed, with S₆ ∼ − 1.0, the vector u used was near (0,1); the precise angle to use was obtained by varying the angle of u in increments of 1/1400 radians away from π/2 and choosing the value that gave the longest correlations in g(Δy). The resulting function is directly comparable with g(Δx) as presented by Bernard and Krauth,¹⁶ except for a factor of 2 in the x-axis arising from a difference in the definitions of σ here and in that work.

The first part of this section includes demonstrations that the SR-GCMC method, and particularly the use of the auxiliary potential, enables simulations in the grand canonical ensemble where other methods pose challenges. Thereafter, validation that the method gives results that are consistent with published work on monodisperse hard disks and are thermodynamically self-consistent in the case of hard disk mixtures will be given. Finally, the applications to the phase behavior of mixtures containing small fractions of large disks will be described, emphasizing the compositions and structures of the ordered and disordered phases at coexistence.

Whereas the use of CBMC to choose small particle positions to fill cavities left by larger particles has some precedents,^9,10 the use of an auxiliary potential to improve packing in this context is a previously unexplored development. To demonstrate the benefit of using an auxiliary potential, we show results from simulations of a hard disk mixture with a diameter ratio d of 6.5. The fugacity f_S of the small disks is set to 1.0 × 10⁴, which for a pure small-disk system produces a fluid with area fraction 0.634. The large disk fugacity was selected at a value (f_L = 1.0 × 10⁹⁷) that gives an area coverage of large disks near 0.2 in the small disk medium. Under these conditions, addition or removal of a single large disk is coupled with removal or repacking of approximately 30 small disks. As evident from Fig. 3, when an auxiliary potential is used, the rate of number fluctuations of the large disks is an order of magnitude greater than in its absence for the same number k of trial attempts. Without the auxiliary potential, even increasing the number of move attempts k by a factor of 10 did not appreciably improve the success probability. For the systems described below, the benefit of the auxiliary potential was more typically near a factor of 2 or 3.

We aim to verify that grand canonical ensemble simulations using the solvent-repacking method reproduce the phase behavior documented recently using canonical simulations.^16,17 Previously, GCMC simulations using simulated tempering⁶ have been reported, but only small system sizes (up to 20 × 20 σ²) were practical due to the scaling of this method; as demonstrated previously¹⁶ and below, this box dimension is smaller than the orientational correlation length of the fluid phase at coexistence, so will induce severe finite-size effects. Here, GCMC simulations on a monodisperse hard disk system in a square box of dimensions 350 × 350 σ² were performed over a narrow range of fugacities f targeted to span the transition, ranging from 3.5 × 10⁵ to 4 × 10⁵, starting either from a hexagonally ordered state or a disordered state.

As apparent from Fig. 4, at a fugacity of f = 3.7 × 10⁵ (μ = 12.82 k_BT), a system initiated in an isotropic state at an area fraction near η = 0.7 increases its area fraction significantly over the course of the trajectory, whereas at a slightly lower fugacity of 3.65 × 10⁵, the area fraction remains stable with an average of η = 0.699. Figs. 5(a)-5(c) show that the jump in area fraction is accompanied by the adoption of a uniform hexagonal orientational order throughout the system. At a fugacity f = 3.6 × 10⁵ the reverse process takes place quickly, while a system initiated in an ordered state and equilibrated at f = 3.65 × 10⁵ (μ = 12.81 k_BT) loses density gradually; at the end of the trajectory, it had undergone a small amount of melting as shown in Fig. 5(f). From these results we may conclude that the fugacity at the order/disorder transition is between 3.65 × 10⁵ and 3.7 × 10⁵, placing a lower bound to the two-phase region at η = 0.699 (corresponding to the isotropic phase at f = 3.65 × 10⁵) and the upper bound at η = 0.718 (corresponding to the ordered phase at f = 3.7 × 10⁵). This range encompasses the two-phase region of 0.700 < η < 0.716 determined by Bernard and Krauth¹⁶ from canonical ensemble calculations. Pressures at points directly above and below the phase transition are similarly consistent with those obtained in Ref. 16, as shown in Fig. 6. These results indicate that the solvent repacking MC approach can be used to equilibrate a hard-disk system within the grand canonical ensemble at densities up to and beyond the ordering transition and to yield information about the phase diagram.

The canonical ensemble simulations of Ref. 16 take advantage of efficient, rejection-free “event chain MC” cluster moves,²¹ which promote sampling of density fluctuations across long distances relative to local moves. Where we expect the grand canonical ensemble to have a real advantage, even over these cluster moves, is for sampling composition fluctuations in a mixed system without requiring the diffusion of particles of different sizes from one region to another. Before discussing applications of the method to phase transitions in binary hard-disk mixtures, we present results of tests for the internal consistency of the method as applied to mixtures of hard disks in the fluid phase. In this test, the fugacity of small disks was kept constant, while the fugacity of large disks was increased. According to the Gibbs-Duhem relation, if volume (box dimensions), temperature, and the chemical potential of one component are held constant, the pressure of the system should vary in a predictable way with increases in chemical potential,

Fixing μ_S by setting the small disk fugacity at f_S = 3.6 × 10⁵, just within the fluid region of the phase diagram for a single-component system (area fraction = 0.698), simulations were performed for several values of the fugacity f_L of large disks with diameter d = 3 times σ. Simulations were initiated with an equilibrated box of pure small disks at the appropriate fugacity, then equilibrated over 2.56 × 10⁹ SRMC moves to introduce the large disks at the appropriate f_L, followed by a production period of 9.6 × 10⁹ SRMC moves. The dependence of ρ_L = N_L/V on f_L could be fit reasonably well with the empirical form

with A = 1.86 × 10⁻³⁵ and B = 2.35 × 10⁻³⁴, as shown in Fig. 7(a). Combining Eq. (18) with Eq. (19) leads to the dependence of pressure on f_L,

Pressures calculated from simulation data using Eq. (17) are in reasonable agreement with the predictions of Eq. (20), as shown in Fig. 7(b). We conclude that the method and its implementation are sound, although pressure calculations in the mixed-disk system are subject to some uncertainty.

After validation of the methods, the effects of introducing low concentrations of large disks on the order/disorder transition was investigated. Although several groups have performed careful canonical ensemble Monte Carlo simulations of mixed-disk systems with moderate or large size ratios,^22,23 these have remained well below the transition (η ≤ 0.60). Phase diagrams for several mixed-disk systems have been investigated using approximations to the equation of state based on an approximate equation of state for the fluid phase and a cell model for the ordered phase,²⁴ and some simulations have been reported of the effects of polydispersity on the freezing transition of equimolar mixtures with small size disparities, up to a size ratio d = 1.4 where freezing was not observed.²⁰ The phase diagrams of non-additive and additive hard disk mixtures have been studied by Guáqueta²⁵ using cluster algorithms^26,27 that are highly efficient for very asymmetric mixtures but only up to a small-disk area fraction of about η = 0.62.²⁶ An important conclusion of Guáqueta is that, in spite of the depletion forces that provide an effective attraction between large disks, there is not expected to be a stable fluid-fluid phase separation in additive hard-disk mixtures, as it will be superseded by a fluid-solid phase separation. Hard-disk mixtures have been the subject of considerable interest as a model system for glassy dynamics;^18,19 the simulations reported here do not, however, reach densities near the glass transition.

Here, we will confine our attention to systems in which large disk number density ρ_L ≪ ρ_S, in which case the large disk might be considered as a solute within the small disk solvent. Efficient simulations outside of this regime should be possible with some modifications to the methods, but the current algorithm is not efficient at equilibrating high area fractions of large disks.

Mixtures at three ratios of large to small disk diameter (d = 1.4, 2.5, and 3.0) were investigated. The small disk fugacity f_S was maintained at 6.0 × 10⁵ in all cases, representing a chemical potential μ_S = 13.3 k_BT, which is 0.5 k_BT above its value at coexistence for the pure small-disk system (assuming that coexistence lies near f_S = 3.65 × 10⁵). The average area fraction for the pure small-disk system at this fugacity was found to be η = 0.729. The pressure of the ordered phase of the pure small-disk system obtained from simulation, 9.63 ± 0.01 k_BT/σ², matches the value predicted by the first-order expansion from the Gibbs-Duhem equation,

when the values from Ref. 16 are used for the coexistence pressure p₀ and the density ρ_ordered,0 of the ordered phase at coexistence.

For a two-component system at fixed μ_S, there should be at most one unique value of μ_L where two phases coexist, at which μ_S, μ_L, and pressure all coincide in two phases with different ρ_S and ρ_L. In practice, if there is a coexistence region, then we expect a range of μ_L over which one phase or the other is metastable over the simulation time scale. Here, mixed-disk simulations were performed over a range of μ_L; at each value, simulations were initiated either from a pure small-disk system equilibrated in the ordered phase or from an empty box. For each value of d investigated, both ordered and fluid phases were stable at one or two of the values of μ_L chosen for simulations. At higher μ_L, the system melted even when initiated in the ordered state. At lower μ_L, systems initiated in a disordered state showed gradually decreasing ρ_L and persistent ordered domains. Once the size of the largest ordered domain approached the box dimensions, these simulations were discontinued and the fluid phase was assumed to be unstable at that μ_L. Figure 8 shows the large disk contents for all trajectories that remained stable in the starting phase, with the scales of μ_L on the x-axis offset to give rough correspondence in the compositions of the isotropic phases. In principle, the exact point of coexistence can be found by applying the constraint that pressure along with μ_S and μ_L must be equal across the coexistence region; unfortunately, uncertainties in pressure calculations were comparable to the difference between pressures over the observed stable range of the fluid phase expected from the Gibbs-Duhem equation (Δp ≈ ρ_L Δμ_L < 0.003 k_BT/σ²).

Table II gives concentrations of large and small disks at values of μ_L where both phases were stable over the simulation timeframe. As expected, the large disks are depleted in the ordered phase relative to the fluid phase. The ratio K_L = ρ_L,fluid/ρ_L,ordered at fixed chemical potential near coexistence varied strongly among disks of different sizes. The partitioning was strongest for d = 2.5, with K_L ∼ 20; next for d = 1.4, with K_L slightly below 10, and weakest for d = 3, with the large disks enriched in the fluid phase by a factor less than 5. The relatively high solubility of d = 3 disks in the ordered phase can be explained by its ability to substitute for an arrangement of seven small disks (one central disk plus six nearest neighbors) in the hexagonal array with little disruption to the packing. The twelve second-shell neighbors lie at distances close to the distance of closest approach of the large and small disks (2 σ) marked by the outer red arc in Fig. 9(a). With relatively little disruption to the extended lattice—some shifting nearer and some shifting farther, they can “solvate” this large disk. The positions of the inner and middle red arcs in Fig. 9(a) suggest that large disks in the d = 1.4 mixture can only substitute for a single small disk in the lattice if all six nearest neighbors shift far to the outer reaches of their distributions, while the d = 2.5 large disks leave too much empty space if they substitute for a single disk and its first-neighbor shell. Therefore, in these mixtures, the large disks are more likely to be found in regions of packing disruption, which in turn affects the spatial and angular correlations between large disk positions. Only the d = 3 mixture features peaks in the distribution of displacement vectors between large disks that mimic those of the underlying small-disk lattice, as shown in Fig. 9(b). It may help, here, that the circle of closest approach for a pair of large disks (outer blue arc in Fig. 9(a)) intersects twelve more peaks in the third neighbor shell of the hexagonal lattice. In contrast, as evident in Figs. 9(c) and 9(d), the relative positions of disks of diameter 1.4 and 2.5 are distributed more isotropically, both in their angular distributions and in their absolute distances.

Looking at the orientationally averaged radial distribution functions shown in Fig. 10, the d = 3 case gives much longer coherence but shorter enrichment, while d = 2.5 shows the greatest and longest-ranged enhancement of local concentration of large disks near other large disks. Although the effective interactions are rather long in range, they decayed over a 10-15 σ length scale that is well below the box dimension. All of these evidence suggest that large disks of diameters 1.4 σ and 2.5 σ induce (or find) regions of local disorder, which then attract more large disks. Individual disks of size 3 σ better able to occupy lattice positions without disrupting the underlying structure and can even pack in pairs while preserving the lattice symmetry to a large degree. The effective attractions among large disks, mediated by their effects on the small-disk lattice, lead to the nucleation and growth of fluid domains once a saturation threshold is reached; for the three large disk sizes studied here, stronger attractions are correlated with lower thresholds.

In contrast to the varied levels of large disks found in the ordered phase at coexistence, the concentration of large disks found in the fluid phase was much less dependent on disk size ratio d. As in the familiar phenomenon of freezing point depression, the addition of large disks to the system is expected to extend the region of stability of the disordered phase to a degree that is, to first order, a function of the mole fraction of solutes but not their specific properties. Given coexistence at a pressure p₀ and small disk chemical potential μ₀ of phases A and B with number densities ρ_0,A and ρ_0,B, the pressure of phase A or B can be written for an ideal-dilute solution of large disks of number density ρ_L at a different μ_S as

Equalizing pressure in two phases containing an ideal-dilute solute near the transition for the pure solvent yields

The difference in chemical potentials between the system at f_S = 6.0 × 10⁵ and the pure disk system at the point of phase coexistence is between 0.48 k_BT and 0.51 k_BT, given the finding above that μ₀ = k_BT ln f_S, where f_S is between 3.6 × 10⁵ and 3.7 × 10⁵. The difference in densities at coexistence, from Bernard and Krauth,¹⁶ is 0.0204 σ⁻². Thus, we expect the difference in large disk concentrations between the phases to be 0.010 σ⁻². In contrast, as shown in Table II, the observed range of values found for Δρ_L (at μ_L where both phases are stable during the simulations) varied from 0.011 to 0.19 σ⁻².

The greatest contribution to this discrepancy is likely to be the non-ideality of the large disks’ interactions in the fluid phase. This non-ideality is illustrated in Fig. 11, which compares the dependence of large (d = 3) disk concentration on fugacity f_L under three conditions. In the absence of small disks, concentration falls below fugacity above f_L = 0.005 (where area fraction is 0.035) due to excluded volume effects. When small disks are introduced at a constant fugacity of f_S = 360 000, repulsions between large disks appear to be canceled out by depletion effects from the small disks, and the ideal-dilute proportionality ρ_L ∝ f_L holds up to ρ_L = 0.014 σ⁻² (large disk area fraction 0.1). At f_S = 600 000, we cannot extrapolate to ρ_L = 0 because the fluid phase is not stable below ρ_L = 0.014 σ⁻²; however, even if we assume ideal-dilute behavior for the large disks at this concentration (which is unlikely), the subsequent rise in ρ_L with f_L is faster than the ideal-dilute model would predict.

Consistent with this thermodynamic evidence that the addition of large disks is cooperative, radial distribution functions for large disks at all d show an enrichment of neighbors within a length of 10 σ or more, as seen in Fig. 12. This enrichment is much reduced or absent in a fluid phase system with f_S = 3.6 × 10⁵ and f_L = 1.0 × 10³³ compared to the system with nearly identical ρ_L at f_S = 6.0 × 10⁵ and f_L = 5.0 × 10³⁴. The qualitative difference in mixing behavior is also evident in snapshots presented in Fig. 13. With the small disks color-coded to indicate the complex value of S₆ in different regions, it is apparent that the mixture with higher f_S (Fig. 13(b)) consists of larger ordered domains, with regions enriched in large disks separating these domains. At lower f_S (Fig. 13(a)), these domains are significantly smaller. For a quantitative measure of domains sizes, the radial correlations of hexagonal packing order within several monodisperse and mixed-disk systems are shown in Fig. 14. In the monodisperse system, the orientational correlation lengths are very sensitive to f_S near coexistence. In the presence of the large disks, the orientational correlation length is shorter than for the pure small disks at coexistence. (This is evident qualitatively in the comparison of snapshots Fig. 5(a) versus Fig. 13(b).) The large disks might thus be playing a surfactant-like role, stabilizing the boundaries between ordered domains. The description of 2-d fluids near the ordering transition as a microphase separated mosaic of ordered domains has been suggested by Patashinski et al.;²⁸ the presence of the large disks at domain boundaries could have interesting implications for the mechanics and dynamics of these systems.

In constant NVT simulations of pure hard disks, Bernard and Krauth¹⁶ found that the correlations between particle positions in the ordered phase at η = 0.718 decayed exponentially with distance over a length scale of 50 particle diameters, while at η = 0.720, these correlations decayed as a power law. From these observations, they concluded that over a range of area fractions from the lower limit of stability of the ordered phase at η = 0.716 to around η = 0.720, the hard disk system is in the hexatic liquid crystalline phase, with only short-ranged positional order. Further simulations showed that systems with stiff, continuous repulsive potentials followed a similar process (continuous solid-hexatic, first-order hexatic-liquid), while softer potentials showed continuous hexatic-liquid transitions.²⁹ Here, we compare the translational order in the pure hard disk systems near coexistence to that in mixtures, keeping in mind that finite size effects and slow relaxation times may influence the results;¹⁷ changes in the positional order are associated with the slow processes of creation, migration, and annihilation of dislocation defects. In Fig. 15(a), the monodisperse hard disk systems show the expected transition from short-ranged translational order at η = 0.718 to quasi-long ranged order at η = 0.721. The same analysis performed on the mixed-disk ordered phases, shown in Fig. 15(b), shows that the presence of the large disks reduces the long-ranged order relative to the pure system at η = 0.728 (at the same μ_S), but that the positional correlation lengths in these impure structures are much greater than seen in pure hard disk hexatics at η = 0.717 or 0.718 and appear closer to 2-d solid structures. Given the uncertainty in the precise point of the coexistence between the isotropic and ordered mixtures and the difficulties in converging these structures, the possibility that a stable hexatic mixture with a smaller positional correlation length exists cannot be ruled out, but these results are suggestive that in these mixtures, the (already narrow) region of stability of the hexatic phase is diminished or entirely swallowed up by the two-phase coexistence region.

The hexatic phase is characterized by short-ranged positional order and quasi-long range orientational order, while the 2-d solid displays quasi-long-ranged positional order and long-ranged orientational order. Since the theoretical basis for the transition from a 2-d solid to a hexatic phase is grounded in the thermodynamics of edge dislocations (see, e.g., Ref. 30 for a concise summary), a method to help visualize these dislocations and variations in orientational order is a useful tool for qualitative insight. In a plot of disk positions as points on a plane with a distorted aspect ratio as in Fig. 16, the structure (in systems where the mean value of S₆ is close to −1) appears as a row of approximately vertical stripes. In relatively ordered regions of the system, each stripe followed from y = 0 to y = 350 ends up shifted along x by an integer number of rows. At η = 0.717 (Fig. 16(a)), these stripes change direction, both along the y axis and along the x axis; a dislocation, appearing as a branch point in a stripe, is highlighted with a red ellipse. In contrast, at η = 0.728, the stripes maintain their order with only slight distortions. The introduction of large disks in the ordered phase at f_s = 6.0 × 10⁵ (Fig. 16(c)) is accompanied by an increase in the orientational disorder, but not to the degree observed at η = 0.717. It is clear that most large disks are not associated with dislocations.

Grand canonical simulations of hard disks and additive hard disk mixtures have been performed covering both sides of the order-disorder transition. A SRMC approach based on the configuration bias algorithm, which uses auxiliary potential to improve success in packing, was used to equilibrate the systems and establish the coexistence conditions. For monodisperse hard disk systems, results were in qualitative and quantitative agreement with recent work^16,17 on the nature of the transition. The transition was observed at a chemical potential of 12.81 ± 0.01 k_BT. In dilute solutions of large disks in a small disk “solvent”, the partition coefficient for large disks across the ordered and disordered phases at coexistence could be measured. At a diameter ratio of d = 3, the density of large disks in the ordered phase was over 20% that in the fluid phase. Partitioning toward the solid phase was less favorable at d = 1.4 than for d = 3, and still less favorable for d = 2.5. By analyzing angle-resolved positional correlations of large disks embedded in the ordered phases, this result could be rationalized based on the greater ability of disks with d = 3 to insert substitutionally within the hexagonal lattice. As large disks were introduced at constant μ_S, the range of stability of the hexatic phase appears to be partially or fully suppressed, superseded by the transition to the isotropic phase. The isotropic phases of mixtures at coexistence with the ordered phase showed significant non-ideality, with domains of mostly small disks sharing common local hexagonal order separated from each other by disordered border regions enriched in large disks. Variations on the current algorithms could be employed to map out broader regions of the phase diagram in these mixed-disk systems. Finally, the results show that the SRMC approach can make grand canonical Monte Carlo simulations practical in a system with some features in common to other solution-phase systems of interest: dense packing of the “solvent” and considerable size difference between “solute” and “solvent”. Whether it can be successfully applied in molecular systems of chemical and biophysical interest remains to be explored.

This material is based upon work supported by the National Science Foundation under Grant No. CHE-1213904. This work used the Extreme Science and Engineering Discovery Environment (XSEDE), which is supported by National Science Foundation under Grant No. ACI-1053575, and the resources of the Cherry Emerson Center for Scientific Computation.