Lattice-switching Monte Carlo method for crystals of flexible molecules

We discuss implementation of the lattice-switching Monte Carlo method (LSMC) as a binary sampling between two synchronized Markov chains exploring separated minima in the potential energy landscape. When expressed in this fashion, the generalization to more complex crystals is straightforward. We extend the LSMC method to a flexible model of linear alkanes, incorporating bond length and angle constraints. Within this model, we accurately locate a transition between two polymorphs of n-butane with increasing density, and suggest this as a benchmark problem for other free-energy methods.


I. INTRODUCTION
The ability to predict which polymorph of a given organic crystal will crystallize or assemble under conditions of interest is still far from routine, despite significant progress in recent years [1][2][3][4][5][6][7].
When studying crystallization in a given model system, an important first step is to establish the available crystal structures and then map the equilibrium phase diagram.This provides essential context for design and interpretation of simulations which probe nucleation, polymorph selection, and crystal growth behavior.
Obtaining the solid phase diagram reduces to computing Gibbs free-energy differences g per atom or molecule between the available polymorphs.Typically one is interested in classical systems, as a fully quantum-mechanical treatment of crystal nucleation is unlikely to be tractable for some time.In the context of crystallization from the melt, one is also interested in solid-solid free-energy differences at temperatures close to melting, limiting the utility of the quasiharmonic approaches often used to compute free energies of "hard" crystalline solids.
Ideally one could compute g directly via thermodynamic integration along a reversible thermodynamic path between polymorphs.Such paths are difficult to obtain in general, but have been realized for simple transformations in atomic solids [8,9].Instead, one normally resorts to the computation of absolute free energies g for each solid phase, most commonly by connecting them via a fictitious path to an Einstein crystal [10].This is largely routine, particularly for atomic crystals [11], although other methods are available [12,13].However, obtaining a small g by subtracting two much larger numbers is less than ideal, particularly as the uncertainly in the absolute g of each phase can be laborious to quantify.
An appealing alternative is to use a so-called phase-switch algorithm, which samples from both phases in a single simulation without the need to construct a complex interphase path [14].g emerges naturally from the population of the samples obtained within each phase, and the associated uncertainty is extracted from the fluctuations in these populations at equilibrium.Within this class, the lattice-switching Monte Carlo (LSMC) algorithm [15][16][17] is appropriate for computing g between two ordered crystalline phases.
To our knowledge, the LSMC method has been applied to atomic crystals where the constituent bodies have only translational degrees of freedom [18], and also to crystals with rigid molecular units [19].In contrast, many models of crystallizing or self-assembling systems incorporate internal molecular degrees of freedom with multiple constraints.In this paper we demonstrate an extended LSMC scheme for the treatment of such systems.

II. LATTICE-SWITCHING MONTE CARLO
In this section we review the standard LSMC method, and discuss an implementation which expresses the lattice switch as a binary selection between two synchronized Markov chains in a manner which naturally extends to molecular crystals or indeed any other system.
We first define generalized absolute coordinates q α j denoting the j th degree of freedom within a configuration which can be identified as "belonging to" phase α.In the context of solid-solid phase transitions, this implies the configuration is within the basin of attraction surrounding a local minimum in the potential energy landscape, corresponding to a periodically repeating crystal structure.We also define fixed reference coordinates Q α j , normally taken to be q α j at exactly the local energy minimum, i.e., the perfect or ideal crystal.The phase index α can take values of 0 or 1 corresponding to the two crystal structures between which one requires a free-energy difference.
In what follows, we denote the currently selected phase index as α = A and refer to the corresponding lattice as "active."The alternative lattice index (α = B = 1 − δ A,1 ) is termed "passive."Beginning from either of these two choices for A, the LSMC scheme consists of the usual MC moves in q A j , sampling from the currently active lattice.An additional move switches the selected phase index while preserving all relative coordinates δq j = q A j − Q A j .Provided the acceptance probability for the switch move is chosen appropriately, the scheme generates samples across both phases (A = 0 and A = 1) in the correct ratio for the ensemble of interest, leading directly to a free-energy difference via Boltzmann inversion.
For simple atomic crystals, simulated at constant volume, the 3N q α j run over the components of fractional coordinates s i for each atom i. N is the number of atoms within the simulated unit cell.Absolute coordinates r i are computed via the 3 × 3 matrix h which contains the three unit cell vectors arranged into columns such that r i = hs i in the usual fashion.In this context we denote the ideal position of atom i in phase α as R α i .Some attention must be given to how the atom indices are mapped between the two structures.In their original study on close-packed hard-sphere systems, Bruce et al. [15,16] identified an optimal mapping as corresponding to the translation of entire atomic planes, altering the stacking sequence from fcc (face centered cubic) to hcp (hexagonal close packed) and vice versa.The switch move hence has a clear geometric interpretation and can be expressed as a translation vector for each atom , which is zero for atoms i lying in planes which are not translated by the switch.A natural implementation of LSMC with this mapping need only store a single set of coordinates for the currently active lattice.The computation of pair distances for the passive lattice is easily accomplished by adding R i + R j to each separation vector r ij = r i − r j .Hence the scheme is formulated in terms of a single Markov chain which can explore both crystal phases.
For more complex crystals, the geometric operation which transforms one lattice to another with minimal change of local atomic environments may involve the rotation of molecules, as well as the stretching, bending, and rotation of molecular bonds.The nontrivial task of designing an optimal mapping must be undertaken for each pair of crystals under consideration.A mapping without a geometric interpretation will likely be suboptimal in terms of maximizing the similarity of atomic environments between the two phases.Indeed Bruce et al. [15] considered random mappings and demonstrated this to be true.
We adopt the following approach in our implementation.At all times we store the coordinates of both the active and passive crystals.Beginning from the reference coordinates (q α j = Q α j ) each MC move is used to update q α j for both α = 0 and α = 1 simultaneously.The move is accepted or rejected based on the energetics of the currently active lattice.The LSMC scheme in this form can be considered as two synchronized Markov chains sampling from different regions of the energy landscape.The switch operation reduces to swapping the identity of the active and passive lattices, and can be attempted at zero additional computational expense.
We stress that this logical (rather than geometrical) interpretation of the lattice switch is mathematically equivalent to the original formulation, but suboptimal in terms of index mapping.However, its extension to more complex crystals is much more straightforward provided generalized coordinates are appropriately chosen (see Sec. III).The mapping of atomic indices between phases can be chosen arbitrarily, however, it is sensible to maximize the similarity of local environments where possible to maximize sampling efficiency.Some consideration must be given to maintaining the synchronization of the two chains over simulations of many billion MC moves in light of the finite accuracy with which real numbers can be represented in the computer.We periodically enforce this synchronization by computing the change in the generalized coordinates from the ideal reference configuration in the active lattice.The coordinates in the passive lattice are overwritten with the result of applying this same change to the corresponding reference configuration This need only be performed periodically (typically every 100 000 MC moves) and as such constitutes a vanishingly small overhead.

A. Non-Boltzmann sampling
The efficacy of LSMC is hampered by a very low acceptance of switch moves when using the Boltzmann acceptance criterion to sample physically meaningful ensembles.In the context of hard-sphere models, the switch can only be accepted when zero hard-sphere overlaps are present in both the active and passive lattices.Bruce et al. [15,16] overcame this problem by introducing a discrete order parameter M, which we define in the present notation as where m({q α j }) denotes the number of overlaps generated from the configuration in phase α.Configurations for which M = 0 define "gateway states" from which the switch move can be accepted.The sampling of these states is enhanced by introducing a weighting function η (M) with units of k B T .
Trial MC moves which introduce overlaps into the currently active lattice can be rejected without further computation.The remaining moves are accepted with probability where the primes denote quantities evaluated after the trial move.Evaluating the acceptance probability hence necessitates computing the number of overlaps in both lattices after every trial move which does not introduce overlaps in the active lattice.This is easily accomplished with shared memory parallelism over the two synchronized Markov chains.The trial volume or lattice vector moves are biased in a similar fashion.The function η is chosen to achieve uniform sampling in M. In our work we adopt the recursion scheme commonly employed in the Wang-Landau [20] method to compute this function.
The overlap parameter M is discrete, with the kth value denoted as M k .Upon visiting a configuration {δq j } with M k , we increment the corresponding weight function η k by δη = ln f and accumulate a histogram h i → h k + 1.When this histogram is "flat," defined as each entry h k lying within 5% of the mean, f is reduced according to f → √ f and the histogram is reset to zero.The process repeats until f reaches a vanishingly small value.Typically we begin with f = 1.0005 and halt the simulation when 1 − f becomes smaller than 10 −7 .
Once we have converged an accurate set of weights {η k }, we simulate using these without further modification to accumulate a converged sampling across M. Unbiased histograms are then recovered from these data via histogram reweighting [21].

B. Hard-sphere system
Before discussing the extension to model molecular crystals, we first demonstrate that our implementation of LSMC in terms of synchronised Markov chains and a purely logical phase switch is consistent with earlier published results.Specifically we study the entropy difference between fcc and hcp packings of N = 216 hard spheres at ρ = 0.7778 (constant volume) and P = 14.58 (constant pressure) where accurate data are available for validation.All results are presented in reduced units.
Following the above procedure for the generation of multicanonical weights, we obtained the function plotted in Fig. 1.This was then used for an extended simulation of 10 9 MC sweeps, where one MC sweep corresponds to N trial displacements of a randomly selected particle.The mappings between particles in each phase where chosen arbitrarily.
Figure 2 plots the current estimate of s as a function of MC sweep number during this simulation.Convergence is obtained to within a tolerance of 10 −3 (reduced units) after only a few million MC sweeps.An accuracy of 10 −5 (marginally better than previously reported LSMC accuracies) is obtained by one billion sweeps, which currently requires only a single core for less than a week on the current generation of desktop hardware.The final result for s is within the range of uncertainly reported in previous calculations.
Further validation tests have been performed on this system at constant pressure, reproducing literature results for the Δs FIG. 2. Convergence of the entropy difference per particle s (circles) between fcc and hcp phases of 216 hard spheres at a reduced density of ρ = 0.7778 sampled using two synchronized Markov chains and a logical phase switch.Dashed lines indicate the limits set by the uncertainly of calculations reported by previous authors [16] using a geometric phase switch.TABLE I. Validation of our synchronized Markov chain implementation of LSMC against previous work on the fcc-hcp free-energy difference [16].The effect of including anisotropic volume moves in the MC sampling is demonstrated to be negligible at the pressure studied.All values of s and g have been multiplied by 10 5 for comparison to Fig. 2. Gibbs free-energy difference g to similar accuracy.It is interesting to note that all previous calculations of g in this system have used constant pressure MC with isotropic volume moves.We have performed calculations using both isotropic and Parrinello-Rahman anisotropic MC to determine if any small differences in anisotropy of particle displacements between fcc and hcp stackings exert a measurable influence on their relative stability.The results of all our validation tests are presented in Table I.The only significant discrepancy is for 216 spheres at P = 14.58 [22].

III. EXTENSION TO FLEXIBLE SYSTEMS
In this section we demonstrate the extension of LSMC to molecular crystals incorporating internal degrees of freedom subject to constraints.In the framework of a purely logical switch, this reduces to a careful choice of relative coordinates {δq} to synchronize between the two lattices.
One choice would involve synchronizing the fractional lattice coordinates of all atoms within the two molecular crystals.This approach may be suitable in the case of two crystals comprised entirely of fully flexible molecules which adopt a uniform conformation.In the more general case considered here, we apply the following principles.
First, all bond length and angle constraints within each molecule must be simultaneously satisfied in both the active and passive lattices.A change in fractional lattice coordinates can preserve such a constraint in one lattice while violating it in another.Coordinates which are independent of the constrained quantity are essential if gateway states which satisfy simultaneously satisfy all constraints in both lattices are to be sampled in a finite simulation.
Second, the coordinates used to specify the position and orientation of a molecule must be independent of its internal degrees of freedom.This greatly simplifies the process of synchronization between two lattices where molecules are in different conformations.Molecular center of mass, or end-to-end vectors, are therefore not suitable as generalized coordinates.
Third, special consideration must be given to rotational degrees of freedom.If a linear molecule j in lattice 1 is aligned parallel to the x axis, then a π/2 rotation about an axis parallel to x will produce little change in the overlap parameter.However, should molecule j in lattice 0 be aligned perpendicular to x, the same rotation will produce a dramatic change in orientation, introducing many new overlaps.It is therefore important to synchronize changes in orientation from the reference configurations {Q α }, and not absolute orientation of molecules relative to a common axis to limit the range of the overlap parameter which must be explored.

A. Linked hard-sphere alkanes
To demonstrate these principles, we simulate a specific example system.In a pair of papers [23,24], Monson et al. have introduced a family of hard-sphere united atom models for normal alkanes.Molecular chains consist of hard sphere beads, diameter d, linked by bonds of fixed length 0.4 d.All bond angles are constrained at 109.47 • .We have implemented the most recent model described by Cao and Monson in which the intrachain overlap is permitted between first, second, and third nearest neighbors.All other bead pairs interact via a hard-sphere overlap potential.We report all quantities in reduced units where d = 1.Note that this differs from the unit system used in Ref. [23] where the unit of length is taken as the diameter of a sphere with equal volume to that of the molecule.
The flexibility of these chains arises from their ability to explore torsional degrees of freedom.The corresponding dihedral angle is restricted to lie within ±17.4 • of zero (trans configuration) or ±10 • of either +120 • or −120 • (gauche configuration).MC moves leading to dihedral angles outside of this range are rejected.
Our MC simulations consist of whole molecule translations and rotations, plus two moves which sample from the conformations available to each of the N chain molecules, torsion, and configurational bias Monte Carlo (CBMC) moves [25].Torsion moves are made by selecting a random dihedral along the chain and rotating through a random angle φ, between −φ max and +φ max .With a probability of 50%, the magnitude of φ is increased by 120 • to allow sampling of transitions between trans and gauche configurations.We have also CBMC moves, but find these to be largely unnecessary for the simulations of short-chain solid phases reported there.When simulating at constant pressure, the anisotropic nature of the solid necessitates the introduction of Parrinello-Rahman [26] moves which perturb a randomly selected subdiagonal entry in the matrix of cell vectors h.

B. Solid butane phases
To illustrate the applicability of LSMC to systems of this kind, we compute free-energy differences between two solid phases which form when the chains consist of four beads, mimicking butane molecules.Previous studies of this model have located only solid-fluid equilibria [27].We are not aware of any previous study of solid-solid phase equilibria in this system.
The structure of solid butane has previously been explored via a combination of neutron diffraction and molecular dynamics simulation [28,29].At high temperatures a plastic crystal (phase I) is formed in which molecules are orientationally disordered about their long axes, adopting a continuous range of angles.This differs from the rotator phases observed for longer alkanes in that these long axes are not aligned, dividing into two mutually orthogonal groups.
At lower temperatures, butane forms two ordered structures, one of which (phase II) is metastable with respect to the other (phase III) under atmospheric pressure.The structure considered in the solid-fluid equilibria explored by Malanoski and Monson [23] corresponds to phase III, in which the long axes of all molecules are aligned.
We have mapped these structures onto the hard-sphere model of Cao and Monson.Some minor refinement of the resulting coordinates was necessary to remove overlaps, preserve the bond and angle constraints, and to ensure dihedral angles lay within the allowable range.All three structures remain mechanically stable within the present model, however, the plastic crystal nature of phase I is lost and the distribution of rotation angles about the long molecular axis is rather narrow.The structures of phases I and III within the present model are shown in Fig. 3.
The calculations reported below are all performed at constant volume.For anisotropic crystal structures such as these, the dimensions of the unit cell are not uniquely defined by the density.We focus on unit cells that correspond to the average cell shape under hydrostatic pressure.To create these, we first compute pressure versus density curves for both phases using anisotropic [26] constant pressure MC simulations and interpolate on these data to identify the pressure of each phase at the density of interest.We then simulate, at this pressure, to identify the average components of each cell vector.This averaged cell is then scaled isotropically to correct for any small residual difference between the averaged and target density, and populated with chain molecules.
This procedure necessarily introduces a source of error in calculations of Helmholtz free-energy difference at constant density, due to the statistical uncertainty in the computed cell parameters.We return to this issue and quantify the magnitude of this error below.

C. Generalized coordinates
Suitable choices for relative coordinates {δq α } in this system are as follows.The superscripts denote the index of a bead within a chain and subscripts denote the chain index from 1 to N. We define δs i to be the fractional lattice coordinates of the first bead on the ith chain, measured relative to its position in the reference configuration.The quaternion w i defines the rotation of the bond vector r 2 i − r 1 i relative to the x axis, and δw i the change in this quaternion needed to rotate this bond into its current orientation from that in the reference configuration.The change in internal dihedral angle of chain i from its value in the reference configuration is denoted δφ i .
Each of the coordinates δs i , δw i , δφ i , for all i, define the 6N + N relative coordinates {δq} synchronized between the two Markov chains simultaneously exploring phase I and phase III.These coordinates are linearly independent, and each MC move modifies only a single coordinate type.The exception to this last property would be CBMC moves, which despite being largely superfluous in the present context, may be of use for MC switch calculations in other systems.This would require the change in chain conformation generated by a proposed CBMC move to be recalculated as changes in the generalised relative coordinates.The expense of this calculation is negligible compared to the CBMC move itself.

D. Results
We first compute the free-energy difference between phases I and III at a constant density of ρ = 0.5.We use the same overlap parameter as for the hard-sphere system, excluding overlaps between spheres which are directly bonded or linked by only a single atom.Overlaps between the first and last beads on a chain are included.An additional overlap is counted for each dihedral angle in the passive phase which lies outside of the allowable range.
Our mapping between atom positions in phases I and III is essentially arbitrary, beyond the criterion that the sequence of beads within a chain is preserved by the lattice switch.In terms of minimizing the range of overlaps which must be explored, this mapping is certainly suboptimal compared to the (unknown) ideal mapping.We mitigate this by storing both h k and η k on a nonuniform grid, with a higher density of points near M = 0 where η (M) is rapidly varying.This requires a minor modification to the weight generation procedure, incrementing both h k and η k in inverse proportion to the width of bin k.To avoid ambiguity in what follows, f refers to the modification factor employed for bins of unit width.
To accelerate the weight generation further, our simulations employ p independent pairs of synchronized Markov chains executing concurrently on a parallel computer.Each chain pair updates a separate copy of the histogram and weighting arrays h k and η k , periodically consolidated via global communication.The efficacy of this approach reduces as the interval between the communication of these arrays increases, evolving the system with increasingly out-of-date bias.In practice the computational cost of communicating data between processors necessitates a compromise, using a finite communication interval of 100 MC sweeps.
Figure 4 illustrates the accelerated convergence the scheme with increasing p.We typically choose p between 8 and 24 for generation of all weights in the butane system, representing a compromise between speed and efficiency.
A suitably converged set of weights η k is obtained after 4 × 10 6 MC sweeps, at which time the modification factor is reduced to ln f = 10 −7 .Convergence of the free-energy difference with respect to the number of MC sweeps subsequently performed using these weights is presented in Table II.The uncertainly of the computed A should be taken in the context of previous calculations by Malanoski and Monson [23] who computed A/N k B T = −0.4(3) in favor of phase III at this density by means of thermodynamic integration from an Einstein crystal reference state [30].
With the improved accuracy of our calculations, we are able accurately resolve the free-energy difference between these two phases as a function of density.Calculations of the freeenergy difference were performed at densities of ρ = 0.50, 0.51, and 0.52, each using at least 2 × 10 7 MC sweeps with 24 walkers and 192 butane molecules.The Helmholtz free-energy difference between phases I and III reverses over this range.Linear regression suggests a transition at ρ = 0.518.A final LSMC simulation was performed at this density, yielding the double-peaked histogram P (M) shown in Fig. 5.The corresponding free-energy difference is A/N k B T = 0.005.This difference from zero is comparable to the finite-size error established in Table II.Combining these two sources of uncertainty with the gradient of the extrapolated fit (inset of Fig. 5), we report the density at which phase III becomes thermodynamically preferable to phase I as ρ = 0.518 (1).The inset shows the variation in free energy with ρ and a fit (dashed line) to the free energies computed at ρ = 0.5, 0.51, and 0.52 (black circles).The free energy difference computed at ρ = 0.518 is shown as a square.Error bars on computed free energies are comparable to the symbol size.
We now return to the uncertainly in our result due to the use of a particular set of cell parameters to realize a target density.To estimate this, we have repeated the calculation of A at ρ = 0.5 using unit cell parameters which differ from the averages used to compute data in Table II by between one and two standard deviations, while maintaining exactly the same density.The resulting free-energy difference is A/N kBT = 0.165(1) suggesting an additional source of error no bigger than 0.007, and certainly much smaller than this when using averaged cell parameters within one standard error of the mean.This does necessitate an upward revision of the uncertainly in the transition density.

IV. CONCLUSION
We have demonstrated that the highly accurate LSMC method can be applied to the study of organic crystals in which the molecules possess both rotational and internal conformational degrees of freedom subject to constraints.By abandoning any geometric interpretation of the lattice transformation move, the application of LSMC to any system is straightforward provided suitable generalised coordinates can be devised.
In the case of the linked hard-sphere model of Ref. [24], we have used LSMC to locate a (previously unreported) transition between phases I and III on increasing density, and propose this as a benchmark problem for free-energy methods in the context of model molecular solids.When combining this information with previously published data [27], one can assert that this transition lies at higher pressures than the equilibrium between phase I and the fluid, but lower than the pressure at which phase III will coexist with the metastable fluid.This is hence an ideal model for studying polymorph selection during crystallization under compression.
While analogous to a phase transition in butane we stress that the present model does not correctly capture the plastic nature of the phase I crystal, lacking rotational freedom about the end-to-end axes of the molecules.If this freedom could be captured within a modified model, the expected result would be an increase in the phase I to phase III transition pressure to offset the greater entropic stability of the plastic crystal phase.Sampling rare transitions between molecular orientations which occur only in one lattice presents an additional challenge to LSMC, which we defer to future work.We also note that the transition from phase III to phase I in actual butane can be purely temperature induced, whereas the present model is athermal.Nonetheless, the trend of the favoring the less open structure at higher density or pressure is correctly reproduced as would be expected.
We have not reported the computation of Gibbs free-energy differences in the butane system.Attempts to do so have identified a particular implementation issue arising when performing LSMC using two cells of hard particles which significantly differ in volume and shape at constant pressure.To satisfy a detailed balance, the probability of accepting a lattice switch between two crystals of different volume must be In the case of the difference cell shapes, a synchronized change in simulation cell vectors leads to different changes in volume of the two simulation boxes, and may result in initially similar volumes diverging as one samples the currently active crystal.Identifying (and biasing to sample) gateway states purely by means of Eq. (2) will hence fail, as the set of states visited at M = 0 becomes dominated by those possessing volumes for which Eq. ( 4) is negligible.For soft particles, where Eq. ( 2) is replaced by the difference in lattice energies at the current relative coordinates, this problem is simply circumvented by subtracting the exponent of Eq. ( 4) from the overlap parameter such that values of zero once again correspond only to configurations where the lattice switch has a high probability of acceptance (see, for example, Ref. [31]).In the case of hard particles as studied here, this would result in sufficiently large and enthalpically favorable volume differences offsetting a non-zero number of overlaps in the passive phase, leading to nonphysical sampling upon accepting a switch.Similar arguments preclude the use of separate branches for the weight function depending on which phase is currently active, as suggested for soft crystals of substantially differing volume by Jackson et al. [18].As our interest in extending LSMC is motivated primarily by soft coarse-grained models for polymorphic molecular solids [32], we do not pursue alternative solutions to this issue here.
We acknowledge that the length of MC simulations reported in this paper would be intractably expensive if applied to detailed atomistic models with explicit electrostatics (possibly incorporating multipole interactions).The LSMC method is hence best viewed as a complement to existing methods based on thermodynamic integration or lattice dynamics, for use where high accuracy is desirable.For example, in coarsegrained models of polymorphic molecular solids, highly accurate free-energy calculations are needed at reference interaction parameters, before refining with thermodynamic perturbation theory to match experimentally obtained calorimetric data on the relative stability for polymorphs.

FIG. 1 .
FIG. 1. Multicanonical weight function η (M) (circles, right axis) and unbiased probability histogram P (M) (triangles, left axis) for 216 hard spheres as a function of the fcc-hcp overlap parameter M at a reduced density of ρ = 0.7778.The flat histogram (squares, left axis) is that obtained directly during the multicanonical sampling.Calculations were performed using two synchronized Markov chains and a logical phase switch.

FIG. 3 .
FIG. 3. (Color online) Structures of butane phase I (left) and phase III (right) in the linked hard-sphere model of Malanoski and Monson.The supercells shown each contain 32 molecules.Beads within the same molecule share color, and coloring of molecules is arbitrary.

FIG. 4 .
FIG. 4. (Color online) Reduction of the weight increment δη = ln f while generating a weight function over the first 3.5 million MC sweeps of an LSMC simulation of the hard-sphere butane model.Data are shown for various values of p, the number of concurrent pairs of chains periodically synchronizing the histogram h k and the tabulated weight η k .The initial reduction factor is f = 0.05.Simulations were conducted with 192 chains at βP = 50.0.

FIG. 5 .
FIG.5.Histogram of M near the phase I to III transition at ρ = 0.518 for 192 chains.The inset shows the variation in free energy with ρ and a fit (dashed line) to the free energies computed at ρ = 0.5, 0.51, and 0.52 (black circles).The free energy difference computed at ρ = 0.518 is shown as a square.Error bars on computed free energies are comparable to the symbol size.

TABLE II .
Convergence with run length and system size of the Helmholtz free-energy difference between phases I and III of the hard-sphere butane model at a reduced density of ρ = 0.5.Negative values favor phase I.