The gas, liquid, and solid phases of dimerizing hard spheres and hard-sphere dumbbells Richard P. Seara) and George Jacksonb) Department of Chemistry, University of Sheffield, Sheffield, S3 7HF, United Kingdom

~Received 24 August 1994; accepted 4 October 1994! The complete phase diagram of a model associating molecule is determined, including the gas, liquid, and solid phases, the regions of coexistence between these three phases, and the location of the critical and triple points. The model molecule is a hard sphere with two very different attractive interactions, one a short ranged and directional attraction and the other a mean field. The first interaction only forms dimers as a molecule can only interact in this way with one other molecule. This saturable attraction mimics hydrogen and chemical bonding. The second interaction is an approximation for the dispersion forces between molecules. Thermodynamic functions for the liquid and gas phases of this model molecule are obtained from an existing theory for associating fluids but a new theory is developed for the solid phase. This is believed to be the first microscopic theory of a model associating molecule in the solid phase. In the low temperature limit no monomers are present; the system is then a fluid or solid of hard-sphere dumbbells. Simulation data are available in this limit and it is shown that in both the fluid and solid phases the theoretical predictions are close to those of simulation. The pressure equation of state for dumbbells is the most accurate theory available for the solid phase. An approximation for the free energy of a solid mixture of spheres and dumbbells is also presented. © 1995 American Institute of Physics.

I. INTRODUCTION

Hydrogen-bonding molecules, such as water, hydrogen fluoride, and the alcohols, exhibit phase behavior which is often very different from that of simpler molecules with only van der Waals attractive forces, such as the noble gases and alkanes. A striking example of this is ‘‘closed loop’’ immiscibility,1 where a mixture is completely miscible at low and high temperatures, but at intermediate temperatures, there is a region where the mixture phase separates. This behavior is only seen in hydrogen bonding systems. One approach to understanding this behavior is to develop a crude model of an associating molecule, which nevertheless reproduces the essential features of the real molecule. This crude model may then be studied in depth by theory and by computer simulation.2 Both these methods enable not only the study of the phase behavior of the model but also give information for the fluid or solid at the molecular level, e.g., for associating molecules it would be possible to determine how many hydrogen bonds have formed. In phase equilibrium studies for pure components and particularly for mixtures, theory is much more convenient than simulation. Even for a pure component the complete gas–liquid–solid phase diagram would require of the order of 50–100 coexistence points. Computer simulation would take a day or more per point whereas the theoretical phase diagrams shown in this contribution were generated and plotted in an afternoon. Association in fluids has been studied using an approximate theory, in both pure fluids3 and in fluid mixtures.1 The solid phase, however, has not been addressed due to the lack of a theory for association in the solid phase. a!

Author to whom correspondence should be sent. E-mail address: [email protected] E-mail address: [email protected]

b!

J. Chem. Phys. 102 (2), 8 January 1995

This lack of knowledge of the solid phase precluded locating the triple point and the study of the fluid–solid phase transition. The triple point defines the lower limit of gas–liquid coexistence, without it it is necessary to truncate the liquid– gas coexistence curve at an arbitrary temperature which is unsatisfactory. This contribution provides a theory for the solid state of a simple model of an associating fluid and uses it in conjunction with a theory of the fluid phase to generate the complete phase diagram. The model molecule is shown in Fig. 1~a!; it is a hard sphere with an attraction site3,4 mediating a strong, typically several times the thermal energy, but short-ranged and highly directional attraction. The short range and limited arc over which sites may overlap to form a bond results in the interaction saturating: The site on a sphere can overlap with at most one other site. Thus only dimers may form, whilst clusters of three or more spheres are precluded. A similar model with two sites has been studied in the fluid phase,3 but it is not studied here because this model forms chains of bonded spheres and it is not obvious how to deal with the solid phase of these chains. The phase diagrams calculated are for a modification of the model of Fig. 1~a!. This model does not exhibit gas–liquid coexistence as the only attractive forces are saturable, so van der Waals attractive forces are included at the mean-field level. The starting point for our theory of associating spheres is a mixture of dumbbells, hard-sphere dimers with the two spheres joined at contact @Fig. 1~b!#, and spheres. In Sec. II, a new theory for a pure solid of dumbbells is proposed and shown to predict the pressure of the solid phase accurately. In combination with an existing theory for a fluid of dumbbells the fluid–solid transition is predicted quite well. This theory is generalized to the mixture and then the free energy of the associating system is derived. Phase diagrams are then

0021-9606/95/102(2)/939/8/$6.00

© 1995 American Institute of Physics

939

940

R. P. Sear and G. Jackson: Dimerizing hard spheres and dumbbells

FIG. 1. The two model molecules considered. ~a! is a hard sphere of diameter s represented by the larger circle, with an off-center spherical attraction site represented by the smaller shaded circle. ~b! is a hard-sphere dumbbell where the two spheres are fixed at a separation of their diameter s .

calculated for the associating spheres with dispersion forces treated at the level of a mean field. II. THE FREE ENERGY OF DUMBBELLS A. Dumbbells in the fluid phase

The fluid phase of dumbbells has been extensively studied by theory4 –7 and by simulation.7,8 The perturbation theory of Wertheim4 yields an equation for the dumbbell fluid which is essentially as accurate as simulation. The result of Zhou and Stell6 is equivalent. The Helmholtz free energy A for N dumbbells is thus

bA 5ln~ r L 6 ! 2112a ex,hs~ 2 r ! 2ln~ g hs~ s ;2 r !! , N

cell theory13 or from a fit to simulation data.14,15 Hard spheres freeze into a face-centered cubic or hexagonal close packed structure, the difference between the two is immaterial here as they have an identical first coordination shell and essentially indistinguishable free energies.15 The problem then is to find the change in free energy when the spheres are combined into dumbbells. We start from the free energy of N hard spheres and assume, as with cell theory, that each is constrained to lie in a cell of volume v 51/r ; these cells do not overlap so that the spheres are distinguishable. The exact shape of these cells is not critical, but an obvious choice is a dodecahedron centered on each lattice site; the sides of the dodecahedron are formed by planes bisecting the vectors from the central sphere to its twelve neighbors. Even at low, or zero density each sphere is still restricted to its cell.16 Constraining the spheres to lie within cells is expected to change the free energy of the solid phase by only a small amount even at the melting density of the solid. The free energy of N constrained spheres is related to the 17 partition function Z (1) N by

b A hs~ N ! 52ln Z ~N1 ! 1N ln L 3 with

~1!

where b 51/k B T, k B is Boltzmann’s constant and T is the temperature. The number density of dumbbells r 5N/V for V the volume. L is the de Broglie wavelength of a sphere. aex,hs is the excess ~over an ideal gas at the same temperature and density! free energy per sphere for a system of hard spheres divided by k B T, and g hs( s ;2 r ) is the contact value of the radial distribution function for hard spheres of diameter s ; both of these quantities are expressed as a function of the number density of spheres which is twice that for the dumbbells. The hard-sphere fluid is well described by the pressure equation of state of Carnahan and Starling9; a ex,hs and g hs consistent with this equation of state are used in Eq. ~1! when the phase diagrams are calculated. The pressure p here and in the following sections is calculated from the volume derivative of the free energy: p52( ] A/ ] V) T . The chemical potential m is then obtained from m 5A/N1 p/ r . Equation ~1! approximates the difference in free energy between a fluid of dumbbells and a fluid in which the dumbbells are broken up into their component hard spheres by N times the free energy change in forming one dumbbell in a fluid of spheres ~see Refs. 6 and 10 and references therein!. B. Dumbbells in the solid phase

The solid phase of dumbbells has been studied less extensively than the fluid phase, but computer simulation data11 and a theoretical equation of state12 are available. However, below we propose a new theoretical expression for the free energy of hard dumbbells in the solid phase. Just as for the fluid phase a perturbative approach is taken, and the difference in free energy calculated between the solid of dumbbells and a solid of spheres, each the same size as the spheres of the dumbbell, at the same volume fraction. The free energy of the solid phase of hard spheres may be obtained from

~2!

Z ~N1 ! 5

E

exp~ 2 b U ~ 1•••N ! 2 b U c ~ 1•••N !! d1•••dN, ~3!

where U(1•••N) is the potential energy function for the system; it is equal to zero if no two spheres overlap and infinity otherwise. U c (1•••N) is the potential function which constrains the spheres to their cells; it is equal to zero if sphere 1 is in cell 1, sphere 2 in cell 2 etc., and infinity otherwise. The cell theory free energy can be obtained from Eq. ~2! by first breaking Z (1) N up into a product of N identical one body partition functions, and then evaluating this partition function by fixing the twelve neighbors of each sphere at their positions in a perfect lattice: The lattice formed at close packing expanded to the desired density.13 It is helpful to look at the low density limit of the system of spheres constrained within N cells of volume v 51/r , this is lim 5 b A hs ~ N ! 2ln~ v ! N 1N ln L 3 5N ln r 1N ln L 3 .

r →0

~4!

In this limit, the system with cells imposed has the same pressure but higher free energy than the system without these constraints. The solid of dumbbells is believed to show solidlike order only in the positions of the centers of the hard-sphere sites, this was considered in the simulation studies of Vega et al.,11 and in more depth by Wojciechowski et al.18 for two dimensional hard dumbbells. For dumbbells of tangent hard spheres the packing is just as in a solid of free spheres; there is thus no driving force for long-range ordering of the bond vectors. The bond vector being that between the centers of the two spheres in a dumbbell. At close packing the spheres are in a face centered cubic or hexagonal close packed lattice, any configuration of bonds which satisfies the basic requirement of joining all the spheres into dimers is then equally likely. Motivated by this observation, we assume liq-

J. Chem. Phys., Vol. 102, No. 2, 8 January 1995

R. P. Sear and G. Jackson: Dimerizing hard spheres and dumbbells

uidlike ordering of the bond vectors. The assumption is that the dimer can ‘‘rotate,’’ i.e., sample all orientations if the excluded volume interactions are switched off. The free energy is then that of a spherical particle, given by Eq. ~4!. A pair of free spheres has two unhindered rotational degrees of freedom. On formation of a dimer from these spheres, one of these degrees of freedom is lost completely as in the fluid phase and the other is unaffected if free rotation of the dimer is assumed. This is all in the limit s 50; with finite s there are excluded volume interactions and the volume available to the spheres of a dumbbell will decrease denying certain orientations to the dumbbell. These excluded orientations are not included in the approximate dumbbell partition function derived below because they correspond to configurations in which a sphere of a dimer overlaps with another sphere, and configurations in which two spheres overlap are excluded when the hard-sphere partition function is calculated. The limits r 50 and s 50 are equivalent; in both cases the result is a system of point particles in the dodecahedron cages. For N dumbbells the volume available for translational motion is that available to a sphere, i.e., v /2—the volume per sphere is half that per dimer which is v . The solid is of N dumbbells and hence 2N spheres. Excluded volume effects lead to a term which is the potential of mean force, i.e., minus the log of the radial distribution function, for two spheres in contact just as in a fluid. As for the fluid it derives from the increase in volume available to the remaining spheres when two spheres are stuck together; at contact, the excluded volumes of the two spheres overlap reducing the total volume denied to the rest of the fluid by the pair of spheres. The free energy for N dumbbells in the solid phase is thus

bA 52ln~ v /2! 1ln L 6 12a ex,hs~ 2 r ! 2ln~ g hs~ s ;2 r !! , N bA 5ln~ 2 r L 6 ! 12a ex,hs~ 2 r ! 2ln g hs~ s ;2 r ! , N

~5!

where a ex,hs is the free energy per sphere divided by k B T in excess to that for s 50, i.e.,

b A hs~ 2N ! 52N ln~ 2 r L 3 ! 12Na ex,hs~ 2 r ! .

~6!

It is instructive to compare Eqs. ~1! and ~5! for the fluid and solid phases, respectively. In particular, note the factor of two inside the logarithm of the solid free energy, in the absence of excluded volume interactions, the translational motion of a dumbbell is restricted to a cell the size of the cells which hold the spheres. Now, g(r;2 r ) for r5 s is related to the pressure by the virial equation:2

b p hs 5114 h g ~ s ; r ! , r

FIG. 2. The predictions of theory and simulation for the reduced pressure isotherm of hard dumbbells. The circles are the fluid phase simulation data of Ref. 7, the squares and diamonds are the fluid and solid phase results respectively of Ref. 11, and the dashed line with the two crosses is the fluid–solid transition predicted by Vega et al.12 The dotted line with two crosses is the location of the transition predicted by the theory of Paras et al.12 The solid and dot–dashed curves are obtained from the theory presented here, i.e., from the free energies Eqs. ~1! and ~5! for the fluid and solid phases, respectively. The solid curve uses cell theory for the reference solid of hard spheres while the dot–dashed curve uses the fit of Hall; both use the Carnahan and Starling equation for the fluid phase.

S

b p5 b p hs~ 2 r ! 2 r 11

where h 5 ps 3 r /3 is the volume or packing fraction. The pressure of a solid of hard spheres p hs may be obtained from cell theory13 or from a fit to simulation data.14 The pressure p of the dumbbells, obtained by differentiating the free energy Eq. ~5! is

r g ~ s ;2 r ! hs

DS

D

] g hs~ s ;2 r ! . ~8! ]r

The pressure in the fluid phase is given by the same expression.

C. Comparison with computer simulation data

The location of the fluid–solid transition has been determined by equating the pressures and chemical potentials in the two phases, see Fig. 2. These two equalities form two simultaneous nonlinear equations in two unknowns, i.e., the fluid and solid densities at coexistence. This pair of equations is solved using a generalization of the method of false position.19,20 The free energy, pressure, and contact radial distribution function used for the fluid of hard spheres are consistent with the Carnahan Starling equation of state.9 In the solid phase, both the cell theory13 and a fit to simulation data14 are used for the hard sphere properties. Using cell theory the fluid–solid transition is predicted at

h f 50.497 99, h s 50.563 07, p * 510.958, m * 528.384, while using the fit to simulation data

h f 50.502 75, h s 50.556 32, p * 511.459, m * 529.386, which may be compared to the simulation prediction of11

h f 50.490, ~7!

941

h s 50.539,

p * 510.2, m * 527.5.

Here, h f and h s are the packing fractions of the coexisting fluid and solid phases, respectively. The pressure and chemical potential at the transition are also given; in terms of the and m * 5 bm reduced quantities p * 5bp s3p/3 2 2 3 ln(L /s). There seems to be some fortuitous cancellation of errors with the cell theory as the transition is nearer to the prediction of simulation than with Hall’s fit despite the solid

J. Chem. Phys., Vol. 102, No. 2, 8 January 1995

942

R. P. Sear and G. Jackson: Dimerizing hard spheres and dumbbells

isotherm drifting off that from simulation at densities near melting. The cell theory pressure equation of state for hard spheres drifts off from the true pressure to a similar extent near melting.12 It is gratifying to see the excellent agreement between the isotherm obtained using Hall’s fit and the simulation data. The theoretical chemical potential in the solid phase is now compared with simulation data. At the density at which the solid coexists with the fluid h 50.539; with the fit to simulation data for the hard-sphere solid we have m*527.339 and with cell theory m*525.095. At a higher density h50.6173: simulation m*543.5, fit to simulation data for reference m*541.41 and cell theory reference m*539.73. With the accurate equation of state of Hall, the errors for the solid phase chemical potential are around 1–5 %. In the coexisting fluid phase for h50.490, with Carnahan–Starling for the hard-sphere fluid, m*527.480.

III. ASSOCIATING SPHERES A. Associating spheres in the fluid phase

Just as with dumbbells the fluid phase of dimerising spheres has been considered extensively. Wertheim,21 and Zhou and Stell22 have derived the same expression for a dimerizing fluid:

bA 1 5a hs~ r ! 1ln X1 ~ 12X ! , N 2

~9!

where a hs is the free energy per sphere of a hard-sphere fluid divided by k B T. X is the fraction of monomers, i.e., spheres which are not bonded. This fraction X is determined by a mass-action equation which differs very slightly between Refs. 21 and 22. From Refs. 3 and 21 we have 15X1X 2 r g ~ s ; r ! KF,

~10!

where K is the bonding volume;23 it is the phase space, i.e., the relative positions and orientations, for which the sites on a pair of spheres overlap without the two hard cores overlapping. The site–site potential is a square well; if two sites are within a cutoff distance there is an energy ehb ~assumed positive! released. F5exp(behb!21 is the Mayer f function of this site–site potential. However, the solid phase of associating spheres will be treated within the approach of Olaussen and Stell;24,25 for consistency we derive the equations for associating spheres in the fluid phase within their approach. This starts from the free energy of a mixture of spheres and dumbbells:

b A5N 0 ~ ln~ r 0 L 3 ! 21 ! 1N d ~ ln~ r d L 3d ! 21 ! 1Na ex,hs~ r ! 2N d ln g ~ s ; r ! ,

~11!

for N 0 spheres and N d dumbbells. The number densities r 0 5N 0 /V and r d 5N d /V, and N5N 0 12N d is the total number of spheres both free and in dumbbells. L d is the de Broglie wavelength of a dumbbell which as such includes a factor from integration over any internal degree of freedom the dimer might possess. a ex,hs is the free energy of hard spheres in excess of the ideal gas free energy, divided by

k B T. The chemical potentials of the spheres m 0 and of the dumbbells m d may be obtained by differentiating the free energy

bm 0 5

S D ]b A ]N0

5ln~ r 0 L 3 ! 1 bm ex,hs~ r ! V,T,N d

2

bm d 5

S D ]b A ]Nd

S

]g~ s;r ! Nd g~ s;r ! ]N0

D

,

~12!

V,T,N d

5ln~ r d L 3d ! 12 bm ex,hs~ r ! 2ln g ~ s ; r ! V,T,N 0

2

S

]g~ s;r ! Nd g~ s;r ! ]Nd

D

,

~13!

V,T,N 0

where mex,hs is the excess chemical potential of the hardsphere fluid. The composition of a mixture, the fraction of dumbbells, is fixed whereas in an associating system the fraction of dimers varies with temperature and density. At a particular temperature and density the fraction of dimers will be that which minimizes the free energy. This fraction is, of course, a narrow distribution about the most probable value but in the thermodynamic limit the contributions to the free energy of states with compositions differing materially from this most probable value will be negligible. Discarding the contribution of these states, the maximum term method, the free energy is that of a mixture with the equilibrium composition. This composition is determined by the equality of chemical potentials

m 0 5 21 m d ,

~14!

the chemical potentials of a free monomer and a monomer which forms half of a dimer must be equal. Inserting Eqs. ~12! and ~13! into Eq. ~14!, we note that the derivatives of g( s ; r ) cancel, and rearrange to form

r d 5 r 20

S D L6

L 3d

g~ s;r !.

~15!

The de Broglie wavelength of the dimer L d is needed to specify the relative stabilities of the monomer and dimer in the low density limit.24 In this limit g( s ; r )51, with that substitution in Eq. ~15!, we have

rd L6 5 5K 0 . r 20 L 3d

~16!

K 0 is the equilibrium constant for the equilibrium between monomers and dimers, it is defined by the left most fraction of Eq. ~16! in the ideal gas limit. Although K 0 is all that is required to specify the strength of association in the system, we would like to relate it to our model. So, K 0 is expressed in terms of the parameters of our model as K 0 5Kexp~ be hb!.

~17!

Inserting K 0 in Eq. ~15! and replacing r d by ( r 2 r 0 )/2, the mass action equation is derived as

r 5 r 0 12 r 20 K 0 g ~ s ; r ! . With the help of Eq. ~15!, Eq. ~11! becomes

J. Chem. Phys., Vol. 102, No. 2, 8 January 1995

~18!

R. P. Sear and G. Jackson: Dimerizing hard spheres and dumbbells

bA 1 5ln~ r 0 L 3 ! 2 ~ 11X ! 1a ex,hs~r!, N 2

~19!

which is equivalent to Eq. ~9!. The pressure is then

b p5 b p hs~ r ! 1 r 2

S DS D ]X ]r

T

1 1 2 , X 2

~20!

where p hs is the pressure of the fluid of hard spheres. The derivative of X may be obtained by differentiating Eq. ~18!. For associating spheres the reduced pressure p * is defined slightly differently, p * 5 b p s 3 p /6.

B. Associating spheres in the solid phase

The problem of associating molecules in the solid phase does not appear to have been considered from the point of view of a microscopic statistical mechanical theory. The utility of such a theory is obvious; without it the triple point cannot be predicted. Although a theory for the fluid phase of an associating molecule can be used to calculate the gas– liquid coexistence curve, without the triple point this curve must be cutoff at some arbitrary temperature. Another possible application is gas–solid equilibria. Just as for the fluid phase the free energy of a mixture of monomers and dimers is required, it is

b A5N 0 ln~ r 0 L 3 ! 1N d ln~ 2 r d L 3d ! 1Na ex,hs~ r ! 2N d ln g ~ s ; r ! .

~21!

bm 0 5ln~ r 0 L 3 ! 111 bm ex,hs~ r !

S

Nd ]g~ s;r ! g~ s;r ! ]N0

D

~22!

, V,T,N d

bm d 5ln~ 2 r d L 6 /K 0 ! 1112 bm ex,hs~ r ! 2ln g ~ s ; r ! 2

S

]g~ s;r ! Nd g~ s;r ! ]Nd

D

~23!

, V,T,N 0

where mex,hs is the chemical potential of a hard-sphere solid in excess to that in the s 50 limit. Using the equilibrium condition Eq. ~14! with Eqs. ~22! and ~23! 2rd

r 20 K 0

5eg ~ s ; r ! .

~24!

The equilibrium constant for a system of point particles constrained within cells is thus e/2 times that for an ideal gas. This assumes liquidlike ordering of the bond vectors. Substituting for r d ,

r 5 r 0 1 r 20 Keg ~ s ; r ! .

~25!

Using Eq. ~24! in the free energy Eq. ~21!

~26!

This expression may be compared with that for the fluid phase Eq. ~19!. The free energies for the fluid and solid phases as a function of r , X, and T differ only in that the solid phase expression is that for the fluid phase plus one, but note that two different definitions are used for the excess free energy of the hard sphere system in the fluid and solid phases. The resulting pressure equation for the solid phase is the same as for the fluid phase Eq. ~20!. The numbers of dimers will differ, however, in the fluid and solid phases at the same density and temperature due to the slightly different mass action Eqs. ~18! and ~25!. C. Phase diagrams

The associating spheres exhibit only one fluid phase whereas in reality associating molecules exhibit a gas and liquid phase. This phase separation is caused by dispersion forces which are missing in our model. These can be included in an approximate way with a mean-field term.2 The mean-field free energy per molecule is simply proportional to the density, i.e.,

b A mf 52 be mfh . N

~27!

The mean-field contribution to the pressure is

The procedure is the same as in the preceding section. Again L d is substituted using Eq. ~16! which remains unchanged. However, there are differences: a ex,hs in Eq. ~21! is the free energy in excess to that with s 50 but with the spheres still constrained within their cells. In this limit, the number density of dimers is not given by Eq. ~16!, as will be seen below. The two chemical potentials are

2

bA 1 5ln~ r 0 L 3 ! 1 ~ 12X ! 1a ex,hs. N 2

943

* 52 be mfh 2 . p mf

~28!

Including these in the free energies and pressures of both the fluid and solid phases will result in the characteristic phase diagram of a simple substance.2,26 The mean-field parameter emf is used to define the temperature scale, i.e., the reduced temperature is T * 5k B T/ e mf . The diameter of a sphere s defines our unit of length. Now, there is only one model parameter which defines the relative strengths of the association and the mean-field attractions: The ratio ehb /emf . The phase diagram of a simple substance consists of three regions: gas–liquid, fluid–solid, and gas–solid coexistence. In generating a phase diagram these three parts are calculated separately. Analytical equations of state such as those presented here possess van der Waals loops in the two phase region for gas–liquid coexistence. We take advantage of this fact in locating the critical point by searching for the temperature at which these loops vanish. Next the gas–liquid and fluid–solid coexistence curves are determined from the critical temperature or a selected upper temperature, respectively to a temperature believed to be below the triple point. Similarly, the gas–solid curve is determined from low temperatures to above the expected triple point. Each of these three curves now contains not only part of the equilibrium phase diagram but also a metastable continuation. The triple point temperature may be found from any of the three points at which the six density curves cross and then the metastable continuations removed leaving the phase diagram. Coexistence is found in the same way as used to locate the fluid–solid transition of dumbbells and by a method which brackets the two densities either by using the spinodal or a guess. Then one of the densities is held fixed and the

J. Chem. Phys., Vol. 102, No. 2, 8 January 1995

944

R. P. Sear and G. Jackson: Dimerizing hard spheres and dumbbells

FIG. 3. The temperature-density projections of three phase diagrams of associating hard spheres. The letters s, l, and g indicate regions where the stable phase is the solid, liquid, or gas, respectively; the regions without a letter are the two-phase regions. The solid curves are for e hb /emf50, i.e., without association, the dashed curves are for e hb /emf50.5, and the dotted curves are for ehb /emf51. The triple and critical points are tabulated in Table I. The bonding volume K was fixed at 4.2731024 s for all three systems. This volume corresponds to a site displaced 0.4s from the center of the sphere with a range of 0.28s .

density of the second phase is varied until the chemical potentials are equal; the density of the second phase is then held constant and the density of the first varied until the two pressures are equal. The procedure is repeated until the densities have converged; this method was found to be robust. See Ref. 27 for a discussion of the solution of nonlinear equations. Fluid–solid coexistence has been determined without the mean-field term ~not shown!; as the strength of the hydrogen bonding increases the fluid and solid densities at coexistence increase and the pressure decreases smoothly and monotonically. The densities increase towards the theoretical values for dumbbells which are slightly larger than the densities for hard spheres; this is a deficiency of the theory as the simulation densities for dumbbells11 are a little below those for hard spheres.16 This increase is presumably also found in a mixture of spheres and dumbbells as the fraction of dumbbells increases; verifying this would be straightforward using the mixture free energies Eqs. ~11! and ~21!. Phase diagrams for ehb /emf5 0, 0.5, and 1 are plotted in Fig. 3, and the corresponding values of X along the coexistence curves are plotted in Fig. 4. Table I contains the temperatures, densities, and pressures at the critical and triple points. The limit of complete association, dumbbells, has been considered elsewhere with a different theory.28 Association, unsurprisingly, increases the critical temperature but note that association affects the triple point before the critical point in the sense that even with the weaker association the triple point is moved up quite substantially. This results in the liquid range, defined by the ratio of the critical and triple point temperatures first decreasing then increasing as the strength of the association increases. Perhaps the most remarkable feature of Table I is the reduction in the pressure at the triple point. Association in the very dilute gas phase costs much more translational entropy than in the dense liquid or solid. Therefore, there is less association in the gas phase ~see Fig. 4! than in the liquid or solid phases destabilizing the gas phase with respect to these denser phases. Another unusual feature is seen for ehb /emf51; the density of the fluid

FIG. 4. The fraction of monomers at coexistence for ehb /emf50.5 ~dashed curve! and ehb /emf51 ~dotted curve!. The bonding volume K is the same as in Fig. 3, and so the dashed curves corresponds to the dashed curves of Fig. 3, and the dotted curves to the dotted curves of Fig. 3. For both values of ehb /emf , above the triple point, the curves are from left to right: solid, fluid, liquid, and gas phases.

phase in coexistence with the solid phase goes through a maximum at around T * 50.12. IV. CONCLUSION

Equations of state have been proposed for the solid phase of hard-sphere dumbbells, for mixtures of these dumbbells with hard spheres, and for associating hard spheres. The equation of state for the dumbbells has been shown to be accurate giving some confidence in the other two. The validity of the dumbbell as a molecular model is doubtful, when two atoms chemically bond to form a diatomic such as nitrogen the result is a molecule which is best represented by a pair of spheres which overlap substantially. A dumbbell in which the spheres are closer than their diameter is not free to rotate. It follows that the phase behavior is very different to that of the dumbbell studied here.11 However, a theory for a solid in which the molecules cannot rotate is possible within the current framework. Indeed, given the use of cell theory ideas it is perhaps more natural. The dumbbell solid is unusual in that while creating a cavity, i.e., an empty space in the lattice of spheres, is as always very unfavorable, rotation of two dumbbells so that the lattice site occupied by a sphere on one dumbbell is now occupied by a sphere from the other dumbbell and vice versa TABLE I. The triple and critical points of dimerizing spheres.a

ehb/emf

a

T c*

hc

p c*

T t*

h

p t* 24

0

0.094 35

0.131

0.0468

0.04370

2.82310

0.5

0.097 23

0.134

0.0481

0.04860

1.1631024

1

0.118 5

0.129

0.0440

0.04905

1.831026

2.831024 0.472 0.575 1.231024 0.482 0.583 1.831026 0.483 0.583

T c* , h c , p c* , T t* , and p * t are the critical temperature, critical volume fraction, critical pressure, triple point temperature, and triple point pressure, respectively. h indicates the packing fractions at the triple point.

J. Chem. Phys., Vol. 102, No. 2, 8 January 1995

R. P. Sear and G. Jackson: Dimerizing hard spheres and dumbbells

results in a configuration which is no more unfavorable than the original one. However, solid phases in which this rotational diffusion did not occur have been studied by simulation.11 It is straightforward to consider the effect of ordering the bond vectors within our approximate theory. For an ordered solid, for example one of the lattices considered by Vega et al.,11 the free energy Eq. ~5! should be modified to yield

bA 5ln~ 2 r L 6 ! 12a ex,hs~ 2 r ! 2ln~ g hs~ s ;2 r ! /6! . N

~29!

Since we restrict the bonding so that a sphere bonds only to a specific other sphere, we should divide the distribution function by the coordination number, here 12. However, rotation of a dumbell by 180° results in a configuration indistinguishable from the original one. Thus, we subtract ln 2 from the free energy per dumbbell. The result is the addition of ln 6 to the free energy; this leaves the pressure unchanged. The free energy Eq. ~29! with the Hall fit for the hard-sphere free energy14,15 predicts the fluid–solid phase transition at

h f 50.564 70, hs 50.628 63, p * 520.784, m * 546.769, which should be compared with the simulation result for the CP1 structure of Vega et al.11

h f 50.554,

h s 50.617,

p * 518.6,

m * 543.5.

The theory predicts densities and so a pressure which is a little too high, just as for the disordered dumbbell solid. The extent of the overestimation is similar in the two cases. The theory for the thermodynamic functions of both associating spheres, and dumbbells is applicable to dimensions other than three. For example, given accurate expressions for the fluid and solidlike phases of hard discs in two dimensions, the theory presented here yields the thermodynamic functions of associating discs and disc dimers. This has been done29 for the two dimensional dumbbells but the results are disappointing. The excellent agreement with simulation seen for the pressure in both fluid and solid phases in three dimensions is absent in two dimensions. This is a little surprising as Wertheim’s formalism has been shown to yield exact results in one dimension.30 The nature of the approximations made in the theory for the fluid phase is well understood, see Ref. 10 and references therein. The poorer understanding of the solid phases in two and three dimensions makes understanding possible deficiencies of the theory in the solid phase more difficult. This lack of understanding, however, applies to contributions to the entropy which are nearly density independent; the calculation of the pressure of the solid phase from simulation is quite straightforward. Given that the theoretical pressure in two dimensions is significantly in error for both the fluid and the solid phases it seems quite likely that the error arises from the excluded volume term. The assumption that the structure of the fluid, i.e., g( s ; r ), is unaffected by forming dumbbells is the only approximation in the fluid phase. Neglecting any change in structure is perhaps less reasonable in two dimensions than in three due to the smaller number of surrounding particles. In the two dimensional solid a disc of a dumbbell is bonded to one of its six neighbors while in

945

three dimensions the sphere is bonded to one of twelve. If this argument is correct theories of the type studied here will be highly accurate in four or more dimensions. The limitation of this study to association of spherical particles into dimers excludes most substances of experimental interest. Clearly, theories of nonspherical molecules and of molecules which bond more than once would be of greater interest. For nonspherical molecules at least the situation looks hopeful, apart from the dumbbells considered here, there has been work on more general nonspherical molecules in solids, see Ref. 12 and references therein. This could be used as a reference system in place of the hard spheres considered here. For a molecule such as water the difficulties are greater; its solid structure is very different from that of hard spheres. This is due to a problem which does not arise with our simple model: packing effects and association compete, they cannot be satisfied simultaneously. Even for a simple two site model which poses few difficulties in the fluid phase, if the angle between the two sites is fixed it may not be possible for both sites to bond simultaneously in a face centered cubic lattice. For dumbbells, we have an accurate equation of state. Can this be extended to nmers with n greater than two? For an ordered solid, writing down the free energy is straightforward for linear, rigid nmers at least:

bA 5ln~ n r L 3n ! 1na ex,hs~ n r ! N 2 ~ n21 ! ln~ g hs~ s ;n r ! /12! 1ln 2.

~30!

How rapidly the accuracy of this expression decreases as n increases is an open question. For nonlinear rigid and flexible chains there will be combinatorial factors associated with the number of ways the molecules may be fitted into the lattice. These factors will not affect the pressure but will move the position of the fluid–solid phase transition. Note that the nonlinear rigid chains must have angles between the bonds holding the spheres together which are consistent with a face centered cubic lattice.

ACKNOWLEDGMENTS

R.P.S. would like to thank the EPSRC and British Petroleum for the award of a CASE studentship. We acknowledge support from the Computational Initiative of the SERC ~Grant No. GR/H58810-C91! for computer hardware.

G. Jackson, Mol. Phys. 72, 1365 ~1991!. J.- P. Hansen and I. R. McDonald, Theory of Simple Liquids ~Academic, London, 1986!. 3 G. Jackson, W. G. Chapman, and K. E. Gubbins, Mol. Phys. 65, 1 ~1988!. 4 M. S. Wertheim, J. Chem. Phys. 85, 2929 ~1986!. 5 W. G. Chapman, G. Jackson, and K.E. Gubbins, Mol. Phys. 65, 1057 ~1988!. 6 Y. Zhou and G. Stell, J. Chem. Phys. 96, 1507 ~1992!; G. Stell and Y. Zhou, ibid., 91, 3618 ~1989!. 7 A. L. Archer and G. Jackson, Mol. Phys. 73, 881 ~1991!. 8 D. J. Tildesley and W. B. Streett, Mol. Phys. 41, 85 ~1980!. 9 N. E. Carnahan and K.E. Starling, J. Chem. Phys. 51, 635 ~1969!. 10 R. P. Sear and G. Jackson, Mol. Phys. 81, 801 ~1994!. 1 2

J. Chem. Phys., Vol. 102, No. 2, 8 January 1995

946 11

R. P. Sear and G. Jackson: Dimerizing hard spheres and dumbbells

C. Vega, E. P. A. Paras, and P. A. Monson, J. Chem. Phys. 96, 9060 ~1992!, 97, 8543 ~1992!. 12 E. P. A. Paras, C. Vega, and P. A. Monson, Mol. Phys. 77, 803 ~1992!. 13 R. J. Buehler, R. H. Wentorf, J. O. Hirschfelder, and C. Curtiss, J. Chem. Phys. 19, 61 ~1951!. 14 K. R. Hall, J. Chem. Phys. 57, 2252 ~1972!. 15 G. Jackson and F. van Swol, Mol. Phys. 65, 161 ~1988!. 16 W. G. Hoover and F. H. Ree, J. Chem. Phys. 49, 3609 ~1968!. 17 J. G. Kirkwood, J. Chem. Phys. 18, 380 ~1951!. 18 K. W. Wojciechowski, A. C. Bran´ka, and D. Frenkel, Physica A 196, 519 ~1993!. 19 F. S. Acton, Numerical Methods That Work ~Mathematical Association of America, Washington, 1990!.

R. P. Sear and G. Jackson, Mol. Phys. 82, 473 ~1994!. M. S. Wertheim, J. Stat. Phys. 35, 19, 135 ~1984!. 22 Y. Zhou and G. Stell, J. Chem. Phys. 96, 1504 ~1992!. 23 R. P. Sear and G. Jackson, Phys. Rev. E 50, 386 ~1994!. 24 K. Olaussen and G. Stell, J. Stat. Phys. 62, 221 ~1991!. 25 R. P. Sear and G. Jackson, Mol. Phys. ~in press!. 26 H. C. Longuet-Higgins and B. Widom, Mol. Phys. 8, 549 ~1964!. 27 W. H. Press, S. A. Teukolsky, W. T. Vetterling, and B. P. Flannery, Numerical Recipes, 2nd edition ~Cambridge University, Cambridge, 1992!. 28 E. P. A. Paras, C. Vega, and P. A. Monson, Mol. Phys. 79, 1063 ~1993!. 29 R. P. Sear and G. Jackson ~unpublished!. 30 E. Kierlik and M. L. Rosinberg, J. Stat. Phys. 68, 1037 ~1992!. 20 21

J. Chem. Phys., Vol. 102, No. 2, 8 January 1995