Effect of polydispersity on the relative stability of hard-sphere crystals
By extending the nonequilibrium potential refinement algorithm(NEPR) and lattice switch(LW) method to the semigrand ensemble, the semigrand potentials of the and structure of polydisperse hard sphere crystals are calculated with the bias sampling scheme. The result shows that the structure is more stable than the structure for polydisperse hard sphere crystals below the terminal polydispersity.
pacs:05.10.Ln, 68.35.Md, 82.70.Dd
A set of hard spheres under thermal agitation constitute a simple yet non trivial model of condensed matter and especially represents an idealization of a very important class of real colloid dispersions. The model has been extensively studied in the past decades. One of the important feature of the model is that the system undergoes a purely entropy-driven first-order phase transition from the fluid phase to a crystal phase at sufficiently high densityT. E. Wainwright (); W.W. Wood (); W.G. Hoover (). Simple estimations of the free energy of different crystal structures reveal that the possible structure could be face-centered cubic () or hexagonal close packed (). However, due to the very similarity in local environments of the two structures, the difference of free energy between them is extremely small and very hard to determine. The determination of the relative stability between the two structures from theoretical calculations has a long historyB. J. Alder (). A clear consensus was reached in the last decade that is the more stable phaseL. V. Woodcock (); P.G. Bolhuis1 (); A. D. Bruce (), however, there are still different views on the problem. The recent results of Pronk and Frenkel S. Pronk () indicate that a moderate deformation of a hard-sphere crystal may make the phase more stable than the phase. Kwak and Kofke S. K. Kwak () investigated the effect of monovacancies on the relative stability of and hard-sphere crystals.
The particle size of most artificial colloidal systems are polydisperse, the polydispersity is defined as the ratio of the standard deviation and the mean of the diameter distribution of particles, which is an intrinsic property of a colloidal system. The polydispersity may significantly affect the thermodynamic and dynamic properties of a hard sphere system, e.g., there exists a terminal polydispersity above which no cystallization can occurP. N. Pusey (); P.G. Bolhuis (); M. Fasolo (), the osmotic pressure of a polydisperse hard sphere crystal is higher than the one of a monodisperse system with the same volume fractionS. Phan (); S. Phan1 (), and there are local fractionations of particle sizes which has a strong retarding effect on nucleationS. Martin (); H. J. Schope (). However, the influence of size polydispersity on the relative stability of and hard-sphere crystals has not been addressed. In this paper we will compute the free energy difference between polydisperse and hard sphere crystals by Monte Carlo(MC) simulations, which suggests that the phase is still more stable than phase below the terminal polydispersity.
The simple simulating method in canonical ensemble is not suitable for a polydisperse system with continuous distribution of particle sizes. The reason is very simple, the simulation system is often too small to realize a given particle size distribution to the designed accuracy. Thus the grand canonical ensemble or semigrandJ. G. Briano () ensemble must be used. In these ensembles the number of particles of each size (characterized by the particle diameter ), (thus the size distribution ) is permitted to fluctuate, therefore, they can simulate a true polydisperse system in the average sense. Comparing with the grand canonical ensemble, the semigrand ensemble is more suitable because the total number of particles is fixed and the insertion or deletion of particles is not needed. The semigrand ensemble is especially more suited to simulate the dense fluid and crystalD. A. Kofke2 (); M. A. Bates (); P.G. Bolhuis (), which provides us a perfect framework to investigate the stability of polydisperse hard sphere crystals.
In the semigrand ensemble, the particle size distribution is not chosen a priori, which is obtained only after the simulation has been performedD. A. Kofke2 (); M. A. Bates (); P.G. Bolhuis (). This is because the imposed physical variables in the simulation are not the composition distribution but the chemical potential deference function (here is the diameter of an arbitrarily chosen reference component), which is a functional of the composition distribution . Consequently, in order to simulate a system with a prescribed distribution the inverse problem has to be solved. Recently, EscobedoF. Escobedo () and Wilding et alN. B. Wilding () have separately shown that the inverse problem can be solved by a histogram reweighting method. Alternatively, a more robust and convenient scheme, the so called nonequilibrium potential refinement algorithm(NEPR), was proposed by WildingN. B. Wilding1 () and works excellent in the grand canonical simulation. We will extend the algorithm to the semigrand ensemble, and use the extended method to determine the chemical potential deference function of an arbitrarily prescribed composition distribution in a semigrand ensemble. The resulting forms of are then used to study the stability of the polydisperse hard sphere crystals.
To avoid any confusion with the Helmholtz free energy, we will refer to the free energy of the semigrand ensemble as the semigrand free energy. In the semigrand ensemble the most stable phase has the lowest semigrand free energy. There are basically two routes to follow in the evaluation of the semigrand free energy. One is the thermodynamic integration routeD. Frenkel (); D. Frenkel1 () which determines the free energy of a system by integrating its derivatives along a parameter space path connecting the system of interest to a reference system( e.g., Einstein solid or ideal gas). The other is the lattice switch method proposed by Bruce et alA. D. Bruce (), from which the free energy difference between monodisperse and hard-sphere crystals can be calculated more directly than the thermodynamic integration method. The method is utilized in the canonicalA. D. Bruce () and isobaric-isothermal ensembleN. B. Wilding2 (). Therefore, We will extend the lattice switch approach to the semigrand ensemble in the present work and use it for the study of thermodynamic stability of the polydisperse hard sphere system.
The contents of the remain of this paper are as follows. In section II we formulate the statistical mechanics for a polydisperse system within the semigrand ensemble. In section III the methodology employed in the work is described and their validity is checked. The computational details and results are presented in section IV. Finally, we present our conclusions in section V.
Ii The semigrand canonical ensemble
The most convenient ensemble in the simulation of a polydisperse system is the so called semigrand canonical ensemble£¨SCE£© though other ensembles can also be used. In the SCE the total number of particles and the volume are fixed while the sizes of each particles can be changed. The average particle size distribution is determined by the chemical potential difference function . First, lets consider a system of hard-spheres in a volume , and the distribution of the diameter of the spheres is . Here we assume that the number of particles is large enough so that the distribution can be well defined. The Helmholtz free energy of the system is
The semigrand canonical free energy(SCFE) is defined through a Legendre transformJ. G. Briano ()
Here is the chemical potential of the reference particle (with diameter ). SCFE is a function of temperature , volume , total number of particles and a functional of . The partition function for SCE is
Here and are the diameter and the thermal wave length of the th particle, respectively. And is the canonical configuration integral
By setting as the excess chemical potential from idea gas, the semigrand canonical partition function can be written in a more symmetrical form
The SCFE of the system is related to the partition function by the following relation
Thus the stable state can be obtained in the semicanonical ensemble by the minimization of the semicanonical free energy, which is the criteria for the stability of the polydisperse hard sphere crystal.
In practical simulation calculations, the diameter of particles is discretised and the corresponding semigrand canonical partition function is
where and are the maximum and minimum values of the particle diameters, respectively. The above discussion can be extended straight forwardly to the polydispersity of other properties of the particles, such as charge dispersity, shape dispersity and mass dispersity etc.
Iii The methods
Under the semigrand canonical ensemble (SCE), the excess chemical potentials for different particles (different diameters) relative to the reference particle(diameter ), , are given. However, in experiments the fixed quantity is the distribution of particle diameters , in order to simulate the experimental controllable system with given particle size distribution, a proper excess chemical potential has to be chosen which can reproduce the required particle size distribution. This is in fact the solution of the functional equation . Wilding has proposed an effective and robust procedure, the nonequilibrium potential refinement (NEPR) algorithmN. B. Wilding1 (), to tackle this problem. The algorithm can be used to solve a wide range of the so-called inverse problemsD. Levesque () such as obtaining the inter particle interactions from experiment measured structure factors. It can also be used in our problem to find the excess chemical potential from the particle size distribution. The original NEPR algorithm was developed in the framework of the grand canonical ensemble, we extended it to the case of the SCE and used it in the calculation reported here. The following is the detailed description of the extension.
Consider a polydisperse system of hard spheres, the range of the particle diameter is , and the diameter of particles can take discrete values, , , , , . When is large enough, the diameter of the particles tends to a continuous variable which can resemble the real polydisperse system. The diameter distribution is , which is normalized in the following way
The excess chemical potential is solved by simulation in a recurrence way. First, initial guess of the excess chemical potential is assigned, then it is modified at every few Monte Carlo steps according to the instant diameter distribution, , which records the distribution of particles at the instant of the simulation, the simulation is terminated when the average of is the same as the required distribution within some tolerance. Then the is the solution of the problem. The initial value of the excess chemical potential is not a vital factor in the calculation process and may be assigned any reasonable values, for example, for all diameters. The detail implement is the following.
iii.1.1 The particle move
There are two kinds of particle moves in the simulation, the first kind is the random displacement of a randomly chosen particle, which is rejected or accepted depends on whether the new position overlaps to other particles or not, the second kind is the expansion or retraction of a randomly chosen particle, which is named as breathing move in literatureM. R. Stapleton (), the probability of acceptance of a breathing move which does not result in an overlap is
where and are diameters of the th particle before and after the test move. The move is rejected if it results in an overlap with other particles.
iii.1.2 The iteration
For a given particle size distribution, the excess chemical potential is calculated by a Monte Carlo iteration procedure. The central quantity in this procedure is the instantaneous particle size distribution , which is the histogram of the particle size distribution at the instant of the simulation and updated during the simulation. Another important quantity is the average particle size distribution , which is the average of the instantaneous particle size distribution in the simulation. The excess chemical potential is updated by the Wilding’s schemeN. B. Wilding1 () for every short intervals. The Wilding’s scheme in this iteration is given by
Here is the given particle size distribution, is a modification factor of the th iteration. For a given modification factor, the average size distribution is also recorded during the simulation. When the difference of the average size distribution and the given particel size distribution is less than a specified value
one loop of the iteration is finished. The modification factor is then reduced by a factor where is a small integer, , and the excess chemical potential of the last iteration is used as the initial input and start the next iteration. The iteration continues till the modification factor reaches a very small value, and the resulted excess chemical potential is then regarded as the solution of the problem. In practical calculations, the convergent criteria for , the threshold and the reduced factor are tuned to reach both high efficiency and accuracy.
The NEPR algorithm for SCE was tested by the simulation of a polydisperse hard sphere fluid and a polydisperse hard sphere crystal with structure. In the tests the particle size distribution is chosen to be the Schultz distribution, which is the most studied model distribution in polydisperse systems. The distribution is
where is the average diameter of the particles and controls the width of the distribution. In the Schultz distribution, the range of the diameter is , however, in a simulation calculation with finite number of particles, cuts off of the up and the lower limit of the distribution may be specified for convenience of computation. In the test studies the effect of the cutoff is not studied, the emphases is on the effect of polydispersity to the physical properties of the system. On the other hand, small particles may enter into the interstitial space of crystals and induce instabilities of crystal structure, this is beyond the subject of this study though it is an interesting subject of research.
In the test simulation, the threshold , the initial modification factor and the reduce factor is . The termination criteria is . For the hard sphere liquid, the volume fraction and the dispersity ; for the hard sphere crystal, the volume fraction and the dispersity . Figure 1 and figure 2 are the simulation results. Figure 1(a) is the calculated excess chemical potential for the truncated Schultz distribution as function of the diameter of particles in the liquid state, figure 1(b) is the comparison of the given truncated Schultz distribution and the distribution generated with calculated excess chemical potential, the agreement between the two is excellent. Figure 2(a) and figure 2(b) are the calculated excess chemical potential and the comparison of given and generated distributions in the crystal state, respectively, the agreement is also excellent as in the liquid case.
iii.2 Lattice switch
The free energy difference between the and the hard sphere crystals is extremely small, in order to obtain a reliable result for the difference, we need to find an accuracy method of calculation. There are different methods suggested in the past in the studies of monodisperse hard sphere crystals and many results were obtained for the problem. The lattice switch method(LW)A. D. Bruce (); A. D. Bruce1 () developed recently is a high precession method in the calculation of the free energy difference. The method has extended successfully to the calculation of the liquid-solid transition of monodisperse hard sphere systemsN. B. Wilding2 () and also used in the studies of soft sphere systemsA. N. Jackson (); J. R. Errington (). In this subsection, we extend it to the SCE.
A detailed presentation of the LW method for monodisperse hard sphere crystals in the canonical ensemble can be found in references A. D. Bruce (); A. D. Bruce1 (). Here we give a quick sketch of the method in the context of polydisperse hard sphere crystals. The system contains hard spheres in volume with periodic boundary conditions. The hard spheres can be in the or the structures respectively, and the spatial positions of the particles are specified by the position vectors or for the th particle in the structure or structure, respectively. The position vectors can be decomposed as
where the subscript may represent or . The and are lattice vectors of the idea structures, and are the displacements from the idea structures.
In principle, the displacement can be any vectors that only constrained by the geometry of the simulation box, however, in the crystal phase with dispersity smaller than the terminal dispersity, the displacements are naturally cutoff in the simulation time scale. We use to represent the phase space of the polydisperse system, here denotes the diameters of particles. Each structure () associates a set of displacements . In a typical simulation, in which a representative subset of the displacements for one structure is sampled, the transition from one structure to another can not happen because the transition probability between structures is extremely small. The spirit of the lattice switch method is that switch the ideal lattice vectors from one structure to another while keep the displacements frozen. The two sets and have a common intersection , which provides a gate to relate the two structures. All allowed(non overlap) configurations accessible by simulation which are associated with and structures can be divided into three subsets, (a). all the displacements allowed by the structure but not allowed by the structure, which we denote as , (b). all the displacements allowed by the structure but not allowed by the structure, which we denote as , (c). the displacements allowed by both the structure and the structure, denoted as .
The semigrand canonical ensemble partition function of the two structures can be written as
where is the potential energy of the system, which is or for the hard sphere system. The SCFE for the structure is
The SCFE difference of the two structures can be written as
Where , and represent the probabilities of three subsets , and respectively.
It is clear from the above discussion that the calculation of semigrand free energy difference from simulation may be proceed as follows: first, the excess chemical potentials for both structures are determined by the iteration procedure described in the last section. Then starting from any structure, particle size distribution and displacement, regard all of the displacement, particle size and the structure index as random variables, and make Monte Carlo moves. In principle, the system will move among the , and states, by recording the number of microstates corresponding to the three macrostates, the probabilities of each macrostate can be obtained and then follows the free energy difference. However, this prescription is not practical in real simulations, the system will trap in either or macrostate in the simulation period simply because the number of microstates of and are much larger than the number of microstates of . In order to overcome this difficulty, the bias sampling can be used. To achieve this goal, we first define an order parameter for the displacement field,
Here and represent the number of overlap pairs of the and structure for all samples of the displacements. The order parameter can take values of , , , , where corresponding to the macrostate . In the simulation process one of the and has to be zero since the domain of random walk is . The free energy difference can be represented by the macroscopic order parameter as
here is the probability that the order parameter takes value . Now we regard each value of the order parameter corresponding to a macroscopic state, biased sample alorithmB. A. Berg (); F. wang (); J. S. Wang () is to sample the system according to a probability so that the rate of visits to every macrostate is basically the same. If the sampling probability is the inverse of , then the rate of visits to each macrostate is exactly the same. Unfortunately, is the quantity we are looking for which is unknown before calculation. The problem was solved by several different methods in the context of density of states calculations, from which the multicanonical methodB. A. Berg () and Wang-Landau methodF. wang () are the most used. These methods, when used in the current problem, amount to starting the calculation with an initial guess of the probability , sampling the system with inverse of this probability, and modify according to the visit rates to the macrostates till the visit rates are constant for all macrostates within a given tolerance, then the resulted and visit rates together will give an accurate estimate of the free energy difference. One of the implementaion is to introduce a weight function , sampling the system with , calculating the probability distribution of the order parameter , , modifying the weight function till is close to constant and then obtain the required probability through
Based on these discussions, the acceptance probability of the sampling is summarized in the following expression
where and are the order parameters, displacement fields and diameters of the system corresponding to states before and after the test move, respectively.
Iv Computational details and Results
In this section we use the extended NEPR and lattice switch method described in the last section to study the stability of the polydisperse hard sphere crystal. The distribution of the hard spheres is chosen to be the Schultz distribution as was used by many researchers. The chemical potential difference obtained from this distribution with lattice can reproduce fairly accurate Schultz distribution for the specified lattice, this is because of the difference between the two structures is very small(see Fig. 3). In fact, if we specify a chemical potential difference, we may produce the same particle size distributions in both of the and the structures within the statistical errors. It is well known that there is a terminal poly-dispersity in polydisperse crystals above which the bulk crystal may not stably existed. The terminal polydispersity obtained from recent simulation is about P.G. Bolhuis (); D. A. Kofke (). In our simulation we set the maximum dispersity to be so that the stable crystal can be simulated. The simulation boxes are set up to suit the idea and crystal, periodic boundary conditions are used in the calculation, the initial configuration are idea and lattices, respectively. The Wang-Landau sampling method was used to obtain a crude estimation of the weight function of the order parameter, then the multicanonical algorithm is used to refine the result. The calculated results are shown in Table 1, for the polydispersity used in the simulation, the lattice has the smaller free energy and more stable than the crystal. In the case of , there are possibilities that the smallest sphere may jump from the cage of its lattice positions so that a defect may be created, to avoid this situation we used a larger volume fraction as shown in the last column of Table 1. The effect of finite size effect was studied in the case of medium polydispersity by enlarging the simulation boxes with fixed volume fraction, it was found that the value of the free energy is affected by the box size slightly but the relative stability is unchanged.
In order to have a clear picture of the relative stability of the structures, we plotted the probability distribution of the macrostates for two different polydispersities in figures 4–5. For convenience of comparison, we plotted the distribution of and in the same half plane, in the case the absolute value of the order parameter is used. From the figures we see that there is a maximum of probability for each structure, and the probability maximum of structure is larger then that of the structure which means that structure is more favorable. Figure 6 shows that the ratio of the probability of “gate” states () to the probability maximum is about which means that the “gate” states will never be reached if a simple sampling scheme is used. Considering the extreme difference of the probability between the “gate” state and the maximum state, the refinement of macrostates and bias sampling are the necessary scheme to obtain meaningful results.
The particle size distribution used in the calculation is the truncated schultz distribution, this is the widely used distribution in theoretical studies. The results based on this distribution may not be extended to all polydisperse systems, however, we expect that it represents a class of polydisperse systems. The distribution dependence of structures of polydisperse systems requires detailed computations of various polydisperse systems. It should also be noted that the accuracy of our calculation is limited by the computation resources, the number of particles is pretty small, typical Monte Carlo steps are million which is still too small to determine accurately the relationship between the difference of semigrand canonical free energy and the polydispersity. On the other hand, the conclusion that the is more favorable than the structure for the polydisperse hard sphere crystals with truncated Schultz distribution of diameters of spheres is clear and reliable.
V conclusion and discussion
Polydispersity of colloid system is common in artificial colloids. We studied here the effect of polydispersity on the stability of structures of colloid crystals, and found that structure is more stable than the structure, which is the same as the monodisperse case with the same calculations. To study the problem we have extended the NEPR algorithm and the lattice switch method to the semigrand ensemble. The extension provides a powerful tool in the studies of other thermodynamical problems of polydisperse systems. A direct application is the determination of the phase diagrams of the polydisperse systems which may replace or at least complement the current Gibbs-Duhem integration methodP.G. Bolhuis (); D. A. Kofke (). The monodisperse system has already studied by the original LW methodN. B. Wilding2 (). The above extension can also be extended to the problems of soft sphere system simulations by some extra techniquesA. N. Jackson (); J. R. Errington ().
Acknowledgements.The work is supported by the National Natural Science Foundation of China under grant No.10334020 and in part by the National Minister of Education Program for Changjiang Scholars and Innovative Research Team in University.
- (1) B. J. Alder and T. E. Wainwright, J. Chem. Phys. 27, 1208 (1957)
- (2) W.W. Wood and J. D. Jacobson, J. Chem. Phys. 27, 1207 (1957).
- (3) W.G. Hoover and F. H. Ree, J. Chem. Phys. 49, 3609 (1968).
- (4) B. J. Alder, B. P. Carter, and D. A. Young, Phys. Rev. 183, 831 (1969).
- (5) L. V. Woodcock, Nature (London) 384, 141 (1997).
- (6) P.G. Bolhuis, D. Frenkel, Siun-Chuon Mau and David A. Huse, Nature 388, 235(1997).
- (7) A. D. Bruce, N. B. Wilding, and G. J. Ackland, Phys. Rev. Lett. 79, 3002 (1997).
- (8) S. Pronk and D. Frenkel, Phys. Rev. Lett. 90, 255501 (2003).
- (9) S. K. Kwak and D. A. Kofke, J. Chem. Phys. 122, 176101 (2005)
- (10) P. N. Pusey and W. van Megen, Nature (London) 320, 340 (1986).
- (11) P.G. Bolhuis and D. A. Kofke, Phys. Rev. E 54, 634 (1996).
- (12) M. Fasolo and P. Sollich, Phys. Rev. Lett. 91, 068301 (2003).
- (13) S. Phan, W. B. Russel, Z. Cheng, J. Zhu, P. M. Chaikin, J. H. Dunsmuir, and R. H. Ottewill, Phys. Rev. E 54, 6633 (1996).
- (14) S. Phan, W. B. Russel, J. Zhu and P. M. Chaikin, J. Chem. Phys. 108, 9789 (1998).
- (15) S. Martin, G. Bryant, and W. van Megen, Phys. Rev. E 67, 061405 (2003).
- (16) H. J. Schope, G. Bryant and W. van Megen, Phys. Rev. Lett. 96, 175701 (2006).
- (17) J. G. Briano and E. D. Glandt, J. Chem. Phys. 80, 3336 (1984).
- (18) D. A. Kofke and E. D. Glandt, J. Chem. Phys. 87, 4881 (1987).
- (19) M. A. Bates and D. Frenkel, J. Chem. Phys. 109,6193 (1998).
- (20) F. Escobedo, J. Chem. Phys. 115, 5642 (2001); 115, 5653 (2001).
- (21) N. B. Wilding and P. Sollich, J. Chem. Phys. 116, 7116 (2002).
- (22) N. B. Wilding, J. Chem. Phys. 119, 12163 (2003).
- (23) D. Frenkel and A. J. C. Ladd, J. Chem. Phys. 81, 3188 (1984)
- (24) D. Frenkel and B. Smit, Understanding Molecular Simulation (Academic, San Diego,1996).
- (25) N. B. Wilding and A. D. Bruce, Phys. Rev. Lett. 85, 5138 (2000).
- (26) D. Levesque, J. J. Weis, and L. Reatto, Phys. Rev. Lett. 54, 451 (1985).
- (27) M. R. Stapleton and D. J. Tildesley, J. Chem. Phys. 92, 4456 (1990)
- (28) A. D. Bruce, A. N. Jackson, G. J. Ackland, and N. B. Wilding, Phys. Rev. E 61, 906 (2000)
- (29) A. N. Jackson, A. D. Bruce, and G. J. Ackland, Phys. Rev. E 65, 036710 (2002)
- (30) J. R. Errington, J. Chem. Phys. 120,3130 (2004)
- (31) B. A. Berg and T. Neuhaus, Phys. Rev. Lett. 68, 9 (1992)
- (32) F. Wang, D. P. Landau, Phys. Rev. Lett. 86, 2050 (2001); Phys. Rev. E 64, 056101 (2001).
- (33) J. S. Wang and R. H. Swendsen, J. Stat. Phys. 106, 245 (2002)
- (34) D. A. Kofke and P.G. Bolhuis, Phys. Rev. E 59, 618 (1999)
FIG. 1:(a), The solved excess chemical potential difference of hard
spheres as function of particle diameter in the fluid
state. (b), The line is the plot of the Schultz function, and the
dots are the particle diameter distribution obtained from simulation
by using the plotted in (a).
FIG. 2: The same as figure 1, for solids.
FIG. 3: The calculated excess chemical potentials for the truncated
Schultz distribution of particle diameters from NEPR method. The
line is for the structure and the points are for the
structure, the difference is smaller then the statistical errors.
FIG. 4: Left: the excess chemical potential as function of the
particle diameter, insert is the particle size distribution for
crystal. Right:the distribution of macrostates (order
parameters) for (open circles) and (filled square)
crystals. The polydispersity , volume fraction
FIG. 5: Same as figure 4, , volume fraction
FIG. 6: The logarithm of the probability distribution of macrostates of hard sphere crystals, the solid line is for the crystal and the dashed line is for the crystal.