5212
J . Phys. Chem. 1990, 94, 5212-5220
FEATURE ARTICLE Freezing of Mixtures: The Density Functional Theory Steven W. Rick and A. D. J. Haymet* Department of Chemistry, University of Utah, Salt Lake City, Utah 841 12 (Received: September 27, 1989; In Final Form: March 1, 1990)
We present the density functional theory for the freezing of several different kinds of mixtures. The phase diagrams of Lennard-Jones binary systems are presented, including models of real binary mixtures of argon, krypton, and methane. The phase diagrams of binary mixtures of (a) hard spheres and (b) charged hard spheres are also presented. The freezing of polydisperse mixtures is also discussed.
1. Introduction
The density functional (DF) theory of classical statistical mechanics' has been developed by Ramakrishnan and Yussouffz and Haymet and Oxtoby3 to describe successfully the freezing transition. The goal of this theory is to predict the thermodynamic conditions under which a liquid will freeze, and to predict the symmetry and density of the crystal with which the liquid coexists. Central to the theory, for both one-component and multicomponent systems, is a thermodynamic functional Taylor series expansion of the free energy of the crystal about the free energy of the liquid. Implicit in this expansion is the assumption of a similarity between the local environments of the crystal and the liquid phases along the coexistence line. Quantitatively, the similarities in the local environments can be measured in terms of a coordination number, the number of nearest neighbors around a particle, which is often the same for both phases. For multicomponent systems, the local environment depends not only upon the location and number, but also the species, of nearest neighbors. Starting in 1983,4 the DF freezing theory has been applied numerically to many one-component systems, such as hard ~pheres,~-'O the Lennard-Jones l i q ~ i d , ~ ,the " . ~hard-core ~ Yukawa fluid,I3 water.I4 and hard ellipsoids and dumbbell^,'^-'^ and even
( 1 ) Evans, R. Paper presented at Les Houches Summer Lectures, Session XLVIII, 1989. Evans, R. Adu. Phys. 1979, 28, 143. See also: Mermin, D. Phys. Rea. A 1965, 137, 1441. (2) Ramakrishnan, T. V.; Yussouff, M. Phys. Reo. B 1979, 19, 2775. (3) Haymet, A . D. J.; Oxtoby, D. W. J . Chem. Phys. 1981, 74, 2559. (4) Haymet, A. D.J . J . Chem. Phys. 1983, 78, 464. (5) Laird, B. B.; McCoy, J. D.; Haymet, A . D. J. J . Chem. Phys. 1987, 87, 5449. (6) Laird, B. B.; McCoy, J . D.; Haymet, A. D. J. J . Chem. Phys. 1988, 88, 390. (7) Baus, M.; Colot, J. L. Mol. Phys. 1985, 55, 653. (8) Tarazona, P. Phys. Rea. A 1985, 31, 2672; erratum 1985, 32, 3148. (9) Jones,,G.; Mohanty, U. Mol. Phys. 1985, 54, 1241. (IO) Curtin. W.; Ashcroft, W. Phys. Rec. A 1985, 32, 2909. ( 1 I ) Marshall. C.: Laird, B. B.; Haymet, A. D. J. Chem. Phys. Left. 1985, 122, 320. (12) Curtin, W . ; Ashcroft, W. Phys. Rea. Lett. 1986, 56, 2775: erratum 1986, 57, 1192. (13) Kloczkowski, A.; Samborski, A . J . Chem. Phys. 1988, 88, 5834. (14) Ding, K.; Chandler, D.; Smithline, S. J.; Haymet, A. D. J. Phys. Reo. Lett. 1987, 59, 1698. (15) Smithline, S . J.: Rick, S. W.; Haymet, A. D. J. J . Chem. Phys. 1988, 88, 2004. (16) Chandler, D.; McCoy, J . D.; Singer, S . J. J . Chem. Phys. 1987,85, 5971. 5977. (17) McCoy, J. D.; Singer. S. J.: Chandler, D. J . Chem. Phys. 1987, 87, 4853.
0022-3654/90/2094-52 12$02.50/0
reformulated for auantum l i a u i d ~ . ~This ~ ? ~article ~ discusses in detail the generali'zation of th'is theory to a variety of multicomponent systems, including binary mixtures of Lennard-Jones (LJ) particles,20hard charged hard sphere^,^^^^^ and carbon and o ~ y g e n , and ~ ~ .also ~ ~ to a polydisperse mixture of hard
sphere^.^^,^^ The freezing theory requires as input the correlation functions of the liquid and a choice of crystal symmetry. Freezing theory can thus be split into two independent steps, summarized by uij(r)
- gJr)
phase diagram
The potential energy between each pair of species i and j is used to calculate the pair correlation functions gij, which in turn are used in the density functional to predict the phase diagram. For all the calculations reviewed here, we have used an approximate integral equation, the mean spherical approximation (MSA), generalized when necessary for continuous potentials and binary mixtures, to obtain the liquid correlation functions. For the case of LJ mixtures, this integral equation input has been compared with the results of "exact" computer simulations of the same liquid mixture.2o The outline of this article is as follows. In section 2 we describe the density functional theory for the freezing of mixtures. We also discuss several crystal symmetries relevant to binary mixtures. The predicted phase diagrams for Lennard-Jones mixtures are presented in section 3. In section 4, we present the results for hard-sphere mixtures, including the extension of results beyond mixtures of fixed, 1 :1 composition to mixtures of arbitrary composition. Coulomb binary mixtures are discussed in section 5. (18) McCoy, J. D.; Rick, S. W.; Haymet, A. D. J . J . Chem. Phys. 1989, 90, 4622.
(19) McCoy, J. D.; Rick, S. W.; Haymet, A. D. J. J . Chem. Phys., in press. Rick, S. W.; McCoy, J. D.; Haymet, A. D. J . J . Chem. Phys., in press. (20) Rick. S. W.; Haymet, A. D. J. J . Chem. Phys. 1989. 90, 1188. (21) Smithline, S. J.; Haymet, A. D. J. J . Chem. Phys. 1987, 86, 6486; erratum 1988, 88, 4104. (22) Xu, H.; Baus, M. J . Phys. C 1987, 20, 2373. (23) Barrat, J. L.; Baus, M.; Hansen, J. P. J. Phys. C 1987, 20, 1413; Phys. Rea. Lett. 1986, 56, 1063. (24) Barrat, J . L. J . Phys. C 1987, 20, 1031. ( 2 5 ) Brami, B.; Joly, F.; Barrat, J. L.; Hansen, J . P. Phys. Lett. A 1988, 132. 187. . ~ .. , (26) Barrat, J. L.; Hansen, J. P.; Mochkovitch, R. Astron. Asirophys. 1988, 199, L15. I I , (27) Ichimaru, S.; Iyetomi, H.; Ogata. S.Astrophys. J . Left. 1988, 334, Li I .
(28) McRae, R.; Haymet, A. D. J. J . Chem. Phys. 1988, 88, 1 1 14. (29) Barrat, J. L.; Hansen, J. P. J . Phys. 1986, 47, 1547.
0 1990 American Chemical Society
The Journal of Physical Chemistry, Vol. 94,No. 13, 1990 5213
Feature Article Polydisperse mixtures of hard spheres (namely, mixtures with an "infinite" number of components) are reviewed in section 6. Some conclusions and speculations are collected in section 7 .
2. Freezing Theory for Mixtures The first density functional calculation for the freezing of a monatomic liquid, using accurate liquid input, was presented in 1983.4 Following related ideas from gas-liquid coexistence theory, developed among others by Yang, Fleming, and G i b b ~ the ,~~ functional used in this work is A m = Jdr1 Mr1) In [P(rl)/pLl - Ap(r1)I Y2Jdrl
Jdr,
dependent crystal density of component i. This functional, when minimized with respect to pi(rl), for all components i = 1, ..., Y , is the grand thermodynamic potential difference between the crystal and the liquid (within the stated approximations). The original freezing calculations proceeded by Fourier expansion of the liquid direct correlation functions and the density d i f f e r e n ~ e . ~A. ~subsequent, popular method of minimizing the functional (2.2), which is derived from the work of Jacobs,32is to parametrize the crystal density as a sum of Gaussians centered on lattice sites pi(r) =
c(lr1 - r2l)&(r1) Ap(r2) (2.1)
where pL is the density of the liquid, p(r,) is the spatially varying density of the crystal, Ap(rl) = [p(rl) - pL],and c(rI2)is the direct correlation function of the liquid. When minimized, the functional ApQ is the grand potential difference between the coexisting liquid at crystal phases at temperature T and chemical potential p . This approach was generalized to mixtures independently by two groups. Hansen and collaborators have studied freezing of binary mixtures using a canonical ensemble formalism, but at first23 restricted the calculation to substitutionally disordered crystals. In the grand canonical ensemble, this leads to a pseudo-one-component functional in which the direct correlation function c(r) in eq 2.1 is replaced by an "effective" function c*(r), which is a linear combination of the three binary mixture direct correlation functions cij(r). Our group developed a formalism for freezing into both ordered and disordered solids,20,21in the grand ensemble. Subsequent calculation^^^-^^ have used the grand canonical formalism, since it avoids unnecessary double-tangent constructions, and calculates directly the properties of interest. Here we summarize the essential ideas of density functional theory for mixtures in the grand canonical ensemble, complete details of which are given in refs 20 and 21. The central idea is to approximate the correlation functions of the crystal by the correlation functions of liquid by using thermodynamic perturbation theory. From knowledge of the correlation functions, it is straightforward to calculate the thermodynamic properties of the two phases and to find the thermodynamic conditions at coexistence (Le., at the freezing point), at which point the liquid and the crystal have the same pressure and temperature, and each component has equal chemical potential in all phases. The input required for the freezing theory is (i) the liquid correlation functions and (ii) the crystal symmetry. The lattice symmetry must be assumed and a variety of symmetries tested to find the most stable crystal. However, it is worth emphasizing that the lattice constant is predicted by the theory. Even for single-component systems, the density functional theory is, of course, not an exact description of the freezing transitions, and its strengths and weaknesses have been analyzed in detail recently.sv6 In addition, it is expected that the perturbation expansion, which lies at the heart of the density functional method, may be tested by using a very recent, nonperturbative freezing theory,31at least for certain "benchmark" systems. The generalization of the density functional (2.1) for a system of Y components is
-EPi
?r3I2ti R
exp[-(r - R
+ Ti)'/€?]
where {RJis the set of lattice vectors of the crystal, t i is the Gaussian width, and T~ is a basis vector of the crystal. The quantity Ai is a lattice-dependent normalization constant, determined by the requirement pi
1 Ai = - 1 d r pi(r) = V A
(2.4)
where A is the volume per lattice site and pi is the bulk crystal density. Here we have made no ad hoc attempt to account for lattice defects9 A formalism for incorporating vacancies and other point defects into DF theory has been proposed by McCoy, McRae, and Haymet." Stable crystal phases are found by minimizing eq 2.2 with respect to the widths t i and the bulk crystal densities pi. The functional (2.2) can be reduced to a more convenient form, writing ApQ as the sum of three terms Awl, AW2,and AW3. The first term Awl is
awl = i JI d r l
lpi(rl) in [pi(rl)/piL~- Api(rl)l (2.5)
I=
Since pi(r) is a periodic function, the integration in eq 2.5 can be reduced to an integral over a single unit cell
Without loss of generality, we take the unit cell to be centered on a molecule. If the Gaussian densities are sharply peaked (till?, 5 0.3), then the Gaussian centered on a neighboring lattice site will not overlap significantly with the Gaussian centered on the unit cell. Hence, only the R = 0 term contributes to eq 2.5 and we obtain
(2.7) where pL is the total liquid density, which for two components is just p I L + p2L. If the Gaussians do overlap then eq 2.7 is not valid; this case, which happens frequently for substitutionally ordered crystals, is discussed fully in ref 20. To evaluate the second part of the functional
we expand the difference between the crystal density and the liquid density in a Fourier series higher order terms (2.2) where
Api(r) =
$ p i exp(ik,-r)
(2.9)
is the liquid density of component i, and pi(rl) is the spatially
where the sum runs over all reciprocal lattice vectors of the crystal. The Fourier coefficients p' may be viewed as the order parameters of this freezing theory. %his expansion, combined with eq 2.3, determines the Fourier coefficients
(30) Yang, A. J. M.; Fleming, P. D.; Gibbs, J. H. J . Chem. Phys. 1977, 67, 74. See also: Saam, W. F.; Ebner, C. Phys. Rev. A 1977, I S , 3566. (31) McCoy, J. D.; Haymet. A . D. J . In?. J . Thermophys. 1989, 10, 87.
(32) Jacobs, R. L. J . Phys C 1983, 16, 273. (33) McCoy, J . D.; McRae, R. M.: Haymet, A. D. J . Chem. Phys. Lerr., submitted for publication.
Api(rl) = [pi(rl) - ~ piL
i
~
l
5214
The Journal of Physical Chemistry, Vol. 94, No. IS, 1990 pi
= pi exp(ik,.r,) exp(-kq2e?/4), = pi - piL, for q
z0
for q # 0 (2.10)
Further manipulations of this term, described in ref 20, lead to the final form of the functional
-a m - - 1 - -I X" p i [ 5 / 2 + In piL + 3 In (r1/2ci) - In A,] PLV
P L i=l
At the freezing point, the grand potential /3Q = - P V / k T has the same value in both liquid and crystal. The coexistence point is found when the minimized functional (2.1 1) is equal to zero, at which point both phases have the necessary requirements for phase coexistence, namely equal pressure, chemical potentials, and temperature. This method is different from, although thermodynamically equivalent to, the method of Barrat, Baus, and H a n ~ e nwho ~~ choose to work with the Gibbs free energy. The above method uses the Grand potential and is demonstrably much easier to use. The two main advantages of the above formulation are the following. First, in order to calculate the Gibbs free energy as a function of its natural variables, the density p must be eliminated in favor of P, T, and the mole fractions x, and x2. This requires an equation of state, which is known analytically only for hardsphere liquid mixtures in the Percus-Yevick (PY) approximation. In general it would have to be calculated numerically through the virial or compressibility equation. Subsequently, in the Gibbs ensemble, the coexistence point must determined by plotting the curve of the crystal Gibbs free energy and the liquid Gibbs free energy, at the same pressure and temperature, as a function of the mole fraction x . The concentration of the coexisting phases is given by the points on the two curves that have the same slope and therefore the same chemical potential. This requires a graphical analysis, or equivalent algebra, and is certainly more tedious and probably prone to numerical error than the grand ensemble formalism. The work of Barrat and c o - w o r k e r ~ ~also ~ , differs ~ ~ , ~ from ~ the approach presented above in that it relates the solid to a (hypothetical) reference liquid with the same composition as the solid, rather than the actual coexisting liquid. For systems which undergo a high degree of phase separation at the freezing point, this technique may prove to be especially useful, since the coexisting phases have significantly different compositions. However, to date applications of this scheme have been published for only two binary systems. hard spheres,23 yielding results similar to the above method, and charged hard spheres,24which do not undergo phase separation at all. For crystals of binary mixtures, five different symmetries have been examined to date: (1) Substitutionally disordered fcc crystal. A substitutionally disordered crystal has the diifferent atoms of the solid mixture distributed randomly on a common lattice, in this case a fcc lattice. The normalization constant, A,, determined by eq 2.22 is just the probability of finding a particle of type i on a particular lattice site. In the disordered fcc case the probability is the number of type i particles divided by the total number of particles, namely A, =
Pi PI
+ P2
(2.12)
Since p I and p2 are varied independently, the disordered crystal need not necessarily have the same mole fractions as the liquid. [In the grand ensemble, it is important to remember that the number of particles changes: only the chemical potentials remain constant.] In other words, this theory can and does describe phase separation. (2) Ordered cesium chloride symmetry. This structure consists of two interpenetrating simple cubic lattices, one for each type of atom. The basis vectors are r1 = '/,a(O,O,O) and 7 , = 1/2a(I , 1 , 1 ) , where a is the simple cubic lattice constant. Since CsCl is an ordered structure, the normalization constant, Ai, is
Rick and Haymet unity, and both bulk crystal densities are the same. (3) Ordered sodium chloride structure. This consists of two symmetrically placed interpenetrating fcc lattices. The basis vectors are T , = 1/2a(0,0,0)and r2 = '/2a(1,I , I ) , where a is the fcc lattice constant, and Ai = I . (4) Ordered zincblende structure. This crystal also consists of two interpenetrating fcc lattices, but with one displaced from the other along the body diagonal of the cubic cell by one quarter of the length of the diagonal. The basis vectors are 7, = '/,a(O,O,O) and r 2 = 1/2a(1/2,1/2,1/2), where a is the fcc lattice constant, and Ai = 1 . ( 5 ) "Fast sphere" phase. This is our terminology for a crystal with the larger atoms fixed on a lattice, in this case a fcc lattice, and the smaller atoms free to flow through the lattice of the larger atoms. This phase is analogous to a "fast ion" phase, although these atoms have no charge. The density of the small atoms is approximated as a constant but in reality the small sphere density should go to zero at the lattice points of the large atoms. Both bulk solid densities are varied independently and are not constrained to equal either each other or the liquid densities. The crystal normalization constant for the large atoms of this phase is unity. This completes the formalism used in this study. It is worth emphasizing that this theory starts from the partition function and makes two well-defined approximations: (1) the MSA for liquid structure, and ( 2 ) the (truncated) density functional theory for freezing. From these well-studied approximations, the phase diagram is calculated directly.
3. The Freezing of Lennard-Jones Binary Mixtures The interactions of many atoms and of simple, approximately spherical molecules can be characterized by Lennard-Jones potentials. The freezing of binary Lennard-Jones mixtures has been predicted by Rick and Haymet.20 These calculations focused on binary mixtures of three substances, argon, krypton, and methane, because these binary phase diagrams have been studied experimentally. The Lennard-Jones potential energy, ui,(r), between a particle of type i and a particle of type j is (3.1) where cij sets the energy scale and uv the length scale. For a binary mixture, there are three potentials: two potentials between "like" components, and an "unlike" or cross potential. In this work we have adopted standard "literature" values for the Lennard-Jones parameters of argon, krypton, and methane, and the simplest mixing rules for the cross potential, without further attempts to optimize the already good agreement with experimental data for the phase diagrams. The cross-interaction parameters are approximated by the combining rules (3.2) (3.3)
These are known as the Berthelot and Lorentz rules, respectively. The correlation functions for the Lennard-Jones mixtures are obtained by using the mean spherical approximation (MSA). The use of the MSA and the combining rules, which are not necessarily the most accurate rules, introduce additional approximations into the calculations. The results for the three binary mixtures of argon, krypton, and methane are shown in Figure I . The stable phase for all the systems is the disordered fcc phase. In Figure 1 the temperature-composition (T-x) phase diagram is displayed for krypton/methane, argon/krypton, and methane/methane mixtures at constant pressure (1 atm), compared to experimental results.w36 In order to calculate the constant-pressure phase diagram, the (34) Veith, H.;Schroder, E. 2.Phys. Chem. A 1937,179, 16. (35) Heastie, R.Norure 1955, 176, 747. (36) Van 'T Zelfde, P.; Omar, M. H.; Le Pair-Schroten, H . B. M.; Dokoupil, Z. Physica 1968, 38, 241.
The Journal of Physical Chemistry, Vol. 94, No. 13, 1990 5215
Feature Article 120
I
P)
110
Kr/CH4
I
0.60
I
c) C H 4 / A r
b) AI/KI
0.55 kT €1
DISORDERED
0.50
CRYSTAL
80
0.45 00
05 'Kr
00
05
XA,
00
05
10
XCH4
Figure 1. Phase diagram for Lennard-Jonesmixtures at p = 1 atm: (2) krypton/methane, theory (ref 20) (solid line) and the experiments of Veith and Schroder (ref 34) (dashed line): (b) krypton/argon, theory (ref 20) (solid line) and the experiments of Heastie (ref 35) (dashed line); and (c) methane/argon with present theory (ref 20) (solid line) and the experiments of Van 'T Zelfde (ref 36) (dashed line).
temperature and density of a liquid at fixed composition was varied until a freezing solution was found with a pressure of 1 atm (as calculated by the virial equation). As shown in Figure la, using even the simplest mixing rules for the potential energy, and literature values of the Lennard-Jones parameters, the theory predicts the general shape of the krypton/methane phase diagram in good agreement with the experiments of Veith and S ~ h r o d e r The . ~ ~ spindle-shaped phase diagram shows that the molecules are miscible in all proportions in the crystal, and that the fraction of krypton molecules increases upon freezing. When simple, Lennard-Jones potential functions are used, the predicted phase diagram for an argon/krypton mixture at constant pressure, shown in Figure 1 b, does not agree with the experimental data as well as the krypton/methane case. The experiments of Heastiejs find another spindle-type phase diagram. In contrast, the LJ potential functions predict the appearance of a eutectic point at xKr= 0.395, where the liquid is in coexistence with two solid phases, one rich in argon and one rich in krypton (Figure 1b). The eutectic temperature, the lowest temperature at which the mixture is at least partially liquid, is 74.7 K. This calculation shows that the predicted phase diagram is sensitive to the exact values of the correlation functions, which in turn depend on the precise functional form and parameters used in the potential energy functions. Additional sources of the discrepancies could be the MSA equation and the combining rules. These sources of error are expected to become more important as the size difference of the particles becomes larger and as the system moves away from the one-component limit. For the methane/argon mixture, the experiments of Van 'T Zelfde et al.36yield a phase diagram with an azeotrope, a minimum = 0.4. More experimental in the freezing temperature, at xCH4 data on these binary mixtures, including solid-solid transitions, has been obtained by The MSA freezing theory predicts a eutectic diagram, although stable solutions vanish between xCH, = 0.5 and 0.65 (Figure 1 c). The solidus for the argon-rich phase cannot be seen on the scale of the figure since the solid is almost a pure argon phase (xCH,).Within the approximations used here, the disappearance of stable solutions in this middle region is to be expected, since the theory treats the crystal as a "perturbation" about the liquid. As the differences in the compositions of the liquid and the solid increase, shown by the large miscibility gap between the liquidus and solidus curves, the solid becomes less well approximated by the coexisting liquid. Therefore, the (truncated) density functional is less valid. Within this formalism, (37) Greer, S. C. Ph.D. Thesis, University of Chicago. Greer, S. C.; Meyer, L.; Barrett, C. S . J. Chem. Phys. 1%9, 50, 4299. Greer, S. C.; Meyer, L. J . Chem. Phys. 1970, 52, 468. Greer, S. C. Phys. Left. 1973, 43A, 73.
0.0
0.5
x2
0.0
0.5
*e
0.0
0.5
1.o
x2
Figure 2. Phase diagram for Lennard-Jonesmixtures: (a) u1/u2= 0.93 and q / t 2 = 1.0; (b) u , / u 2 = 0.90 and q/t2 = 1.0: and (c) ul/u2 = 1.0 and q / t 2 = 0.9.
it is expected that higher order terms in the perturbation series become more important in these circumstances. The series may converge more rapidly if the expansion were made about a liquid with the same composition as the solid, rather than the coexisting liquid, and this idea is currently under investigation. The krypton/methane phase diagram (Figure la), where the = 1.198) and difference in energy parameters t is large (eKr/tCH, the u difference is small (uK~/uCH, = 0.965), undergoes very little phase separation and forms nearly ideal solid solutions. The methane/argon phase diagram (Figure IC), where the t difference is small (tCH4/tAr= 1.131) and the u difference is relatively large (uCH4/uAr = 1 . I 1 l), behaves very nonideally. From a comparison of the two phase diagrams, it is apparent that size differences play a larger role than energy differences. This is consistent with the notion that geometric factors play an important role in freezing. In order to separate the effect of one of the parameters from the other, systems with identical c parameters and different u parameters were examined. In Figure 2a the phase diagram of a mixture with u2/u1 = 0.93 is displayed. It is a completely miscible type phase diagram, with an azeotrope. The solidus curve is plotted as a dotted line and rises above the liquidus curve around x2 = 0.4. This is unphysical, of course, and is a result of a slightly inaccurate prediction of the solid composition. For example, at a liquid composition of x2 = 0.4, the predicted solid composition is 0.415. Since the slope of the liquidus and solidus curves is negative at this point, the solidus line rises above the liquidus. The important feature is that the liquid will freeze into a solid of about equal composition and small errors in the prediction of this composition give rise to unphysical looking phase diagrams. In the regions rich in small particles, there is some phase separation because, in general, large particles are not as soluble in a solid of small particles as small particles are soluble in a solid of big particles. Systems with u 2 / u 1= 0.90 have a very different type of phase diagram (Figure 2b). This phase diagram has a eutectic point at x2 = 0.39. The solid phase rich in small atoms contains very few large atoms: they are almost completely immiscible. In Figure 2c the phase diagram of a mixture with u 2 / u 1= 1 and t 2 / t , = 0.9 is shown. Very little phase separation takes place, even though the t parameters differ substantially. Energetically, there is nothing to gain from phase separation for this system, because from the Lorentz combining rule (eq 3.3), for all distances the potential for unlike particles is equal to the average of the two like interactions. It is only when the particles have different diameters that phase separation can become favorable. The density functional theory, combined with the MSA closure, is able to predict successfully the freezing of simple Lennard-Jones mixtures with similar sizes, such as Kr/CH4 (Figure la). For mixtures of less similar particles, such as CH4/Ar (Figure 2c), small errors in the representation of the forces between the molecules, and perhaps also errors arising from the MSA, lead to less accurate phase diagrams. As the difference in the molecular
5216
The Journal of Physical Chemistry, Vol. 94, No. 13, 1990
Rick and Haymet
1 .o
o.06P7
0.8
0.05
5
0.6
T/K
t
0.95
f
0.93
1
01
04
0.2
0 0.2 0.4 0.6 0.8 0 0.2 0.4 0.6 0.8 1.0 XI
0.0
0.4
0.6
1.o
0.8
1.2
1.4
f Figure 3. Phase diagram for the I:] liquid hard-sphere mixture with size ratio usrmsll/ubig. Plotted are the coexisting liquid and solid density for
the big spheres for the disordered fcc phase (solid line), as well as the metastable solutions for crystal of CsCl (dotted line), NaCl (short dashed line), and fast sphere (long dashed line) symmetry. size increases, phase separation in the solid phase becomes more favorable, due principally to geometric factors. This results in phase coexistence between a liquid and a solid which differ greatly, and therefore the present version of the perturbative density functional expansion becomes less accurate.
4. The Freezing of Binary Hard-Sphere Mixtures The freezing of binary mixtures of hard spheres with different diameters has been studied by several groups.2&z2*25 Hard-sphere mixtures are easier to study than Lennard-Jones mixtures because the temperature scales out of the problem, and because a relatively accurate, analytical expression for the liquid correlation function is available from the Lebowitz solution of the Percus-Yevick (PY) equation.38 The freezing results of Rick and Haymet20 for hard-sphere mixtures are shown in Figure 3, in which the density of the large spheres of the coexisting liquid and solid phases is plotted as a function of the size ratio, O,mall/ubi In the PY approximation, the stable solid for all size ratios is t8k disordered fcc phase. Also shown are the metastable solutions for other phases, including the fast sphere, NaCI, and CsCl symmetries. For the ordered NaCl and CsCl phases, the composition of the = 0.5. However, compositions of the solid is fixed at Xbig = xsmall stable fcc and metastable fast sphere phases are not fixed and phase separation can (and does) occur. As the size ratio of the spheres decreases, the coexisting solid becomes richer in the large spheres. Therefore, the density of the small spheres at phase coexistence drops considerably as the size ratio is decreased. At usmall/ubg= 0.75, for example, the composition of the solid is Xbig = 0.8. These results differ from the results of Brami et al.,25who studied charged hard-sphere mixtures, which will be discussed in the next section. In the limit of infinite temperature, charged hard spheres have the same structure as hard spheres, but with the added constraint, due to charge neutrality, that xbjg = xsmall = 0.5. Under this constraint, Brami et aLz5found the same ordered phase solutions as those shown in Figure 3, but the disordered fcc solutions vanish at usmall/ubig = 0.84. There is a well established, empirical rule, due to HumeR ~ t h e r ywhich , ~ ~ states that disordered metal alloys are unlikely to be stable if the size ratios differ by more than 15%. For PY hard spheres, we observe no sudden change in the solid-phase solubilities on the spheres near usmll/ubig = 0.85 and it seems that
Figure 4. Phase diagram for the hard sphere mixtures at the normal pressurep = lkB/ubiiK: (a) uyMll/ubis= 0.95, (b) u,,ll/ub = 0.93, (c) urmall/ubig = 0.92, and (d) uSmall/ubig = 0.90.
the Hume-Rothery rule is not a result of strictly geometric factors. Perhaps strong intermolecular interactions are needed to overcome the entropy gained in forming a disordered (rather than ordered) solid. For the case of the constrained hard spheres, the HumeRothery rule seems to be followed, implying that forces which discourage phase separation might be present in systems which obey the Hume-Rothery rule. The freezing of hard-sphere liquid mixtures with concentrations other than 1:1 has also been examined by Rick and Haymetm and Barrat, Baus, and Hanseat3 For a fixed size ratio, the freezing point was determined for all liquid compositions of that mixture. Shown in Figure 4 are four temperature-composition phase diagrams as calculated by Rick and Haymet.20 The temperature is calculated from the virial equation, fixing the pressure arbitrarily at p = l k e / U b i g 3 K . The phase diagram of a “small/abig = 0.95 mixture (Figure 4a) is a spindle-type phase diagram with very little phase separation, similar to (but with even less separation than) the Kr/CH4 mixture shown in Figure la. As the size ratio is decreased to 0.93 (Figure 4b), an azeotrope appears at xsmall = 0.17. When bs,,ll/bbig = 0.92, a eutectic point occurs at xsmall = 0.19 (Figure 412). The usmall/Ubi = 0.90 mixture also has an eutectic (at xmu = 0.20), but now t i e solubility of the big spheres in the solid of mostly small spheres is almost zero (Figure 4d). These results agree with the results of Barrat, Baus, and H a n ~ e n . ~ ~ The azeotrope for the 0.93 mixture and the eutectics for the 0.92 and 0.90 mixtures appear in a much different region of the phase diagram than the eutectic points of the Lennard-Jones mixtures. Both of the Lennard-Jones phase diagrams (Ar/Kr, Figure la; u2/uI = 0.9, t z / t l = 1 .O, Figure 2b) have eutectics closer to the middle of the phase diagram. Figure 2b has an eutectic where the composition of the larger particle equals approximately = 0.6. From this comparison, it 0.4, which corresponds to xsmall is clear that hard-sphere and Lennard-Jones mixtures, while displaying some similar characteristics, differ in terms of the solubilities in the solid phase. The hard-sphere mixtures show a greater miscibility than the Lennard-Jones mixtures. The greater miscibility of hard spheres is also demonstrated by the position of the eutectic points. The Lennard-Jones mixture of a 1:l liquid would freeze into an almost completely phase separated solid, rather than into a solid containing both spheres in near-equal amounts (as the hard-sphere mixture does). 5. The Freezing of Coulombic Mixtures At least two different Coulombic mixtures have been studied via the DF freezing theory: oppositely charged hard sphere^^^*^^^* and the cosmologically significant system of carbon/oxygen
mixture^.^^.^' (38) Lebowitz, J. L. Phys. Reo. 1964, 133, AB95 (39) Hume-Rothery, W.; Smallman, R. E.; Haworth, C. W.The Structure of Metals and Allovs: The Metals and Metallurgy Trust: London, 1969.
(40) Rick, S. W. Ph.D. Thesis, University of California, 1989
The Journal of Physical Chemistry, Vol. 94, No. 13, 1990 5217
Feature Article
50 fcc + liquid
1
-
CsCl + liquid
25
0.6
P’
0 2 4
0
0.4
25 0.2
0.0
0
0
10
20
30
40
0.5 1.0 1.5 2.0 2.5 0.5
50
P’
1.0 1.5 2.0 2.5 3.0
P q3
Figure 5. Phase diagram for the charged hard-sphere mixture showing the regions of solid-liquid coexistence for the different solid structures. Plotted are the size ratio of the spheres versus the reduced inverse tem-
perature, 8’.
Figure 6. Density-inverse temperature phase diagram for the charged hard-sphere mixture: (a) us,,l,/ubie = 1.0, (b) us,ll/abig = 0.85, (c) usmall/ubi8 = 0.75, and (d) usmadubig= 0.4.
where a, is the hard-sphere diameter of ion a, e is the elementary charge, and Z, is the charge of a in electron units. The mixtures discussed here possess unit charge (Z, = i l ) . In addition to the work of Brami et aLZSand Barrat,24Rick and Haymet@ have calculated independently the phase diagram for charged hard spheres of different diameters, using correlation functions from the analytic solution to the mean spherical approximation (MSA) of Hiroike.4’ The MSA solutions, which reduce to the PY hard-sphere solutions in the appropriate limit, have been shown to be reasonably accurate near the infinite temperature limit.42 The results of these calculations are shown in Figure 5, which is a plot of the range of solid-liquid coexistence with reduced inverse temperature p* = 2/(ubi,ks7‘)on the x axis and the size ratio on t h e y axis. The figure displays the regions where stable solutions exist for a particular symmetry. The regions where no stable solutions were found are marked by the parallel lines. In the /3* = 0 limit, as mentioned in the preceding section, the Coulomb interactions vanish, leaving only the hard-sphere = xbig = 1/2. With this constraint, parts, with the constraint xsmall the fcc results vanish below usmall/abig = 0.85. The fcc results are temperature independent, since the disordered nature of the fcc lattice cancels the Coulomb interactions. The ordered phases are not temperature independent, as can be seen in Figure 6, which displays density versus /3* for four sets of size ratios. At a,,,ll/abig = 1 (Figure 6a), the CsCI structure becomes the stable phase at @*= 21.6 in agreement with the results of Barrat,24who calculated the phase diagram for charged hard spheres with equal diameters. At ‘&mall/bbig = 0.85 (Figure 6b), the CsCl phase becomes stable at @*= 18.9. At this size ratio, which is the smallest size ratio with stable fcc solutions, the fractional density change into the fcc structure is very small. Below t&mll/abi = 0.85 and above a,,,ll/q,ig = 0.6, the CsCl structure is the only stable phase (Figure 6c). As the temperature is decreased (or o* increased) and the strength of the ionic interactions, which stabilize the CsCl lattice, are increased, the CsCl phase becomes more stable. It freezes at a lower density and stable solutions can be found for a wider range of size ratios. These
calculations are consistent with the results of Brami et The crystal of NaCl symmetry has the opposite temperature dependence (Figure 6d); as the temperature is decreased the crystal becomes less stable until stable solutions vanish a t Usmai\/Ubig = 0.39 and @*= 45. At the range of size ratios for which the NaCl lattice is the stable phase and by use of the MSA correlation functions, the repulsive same ion interactions outweigh the favorable different ion interactions. Therefore, the ionic interactions destabilize the NaCl lattice at usmll/u6g > 0.5 and increasing the temperature decreases the lattice stability. Brami et aLzs appear to have found the opposite dependence: the NaCl solid becomes more stable as temperature is decreased. The discrepancy between the two independent calculations is extremely difficult to explain, since (i) both use the same liquid input and (ii) both yield the same result at @*= 0 and also for CsCl and fcc at asmall/bbig > 0.6. It should be noted that the MSA solution at Usmall/abig < 0.5 and at finite temperatures is not proven to be accurate, and better liquid structural information may be required. A second Coulombic mixture which has been studied with density functional freezing theory is a model of the carbon/oxygen m i x t ~ r e . ~ ~Completely *~’ ionized mixtures of carbon and oxygen are important materials in white dwarf stars, and the C / O mixture in the dense interior of the white dwarf is believed to undergo a freezing transition. The composition of the resulting solid phase will effect the gravitational energy of the white dwarf, and hence affect the cooling rates, luminosity functions, and supernova mechanisms. From an analysis of the cooling rates and the observed luminosity function, one can even attempt to estimate an age of the universe! In the density functional calculations of Barrat, Hansen, and Mochkovitch26 and Ichimaru, Iyetomi, and Ogata,27 the C/O mixture is treated as a mixture of two positively charged ions (2, = 6, Zo = 8) in a neutralizing background of degenerate, incompressible electrons. The two studies differed slightly in the input used for the correlation functions and, not surprisingly, found different results. Barrat et a1.26found a spindle-type phase diagram (similar to Figure 4a) whereas Ichimaru et aLz7found an azeotropic phase diagram similar to Figure 4b, with the azeotrope at XO = 0.16. In the cosmologically significant region around equal composition of carbon and oxygen, both studies find the same qualitative behavior: the fluid freezes into a solid which is slightly richer in oxygen than the fluid phase. These results are used to support the assumptions made by Winget et aI.$j who calculated the white
(41) Hiroike, K. Mol. Phys. 1977, 33, 1195. (42) Abramo, M . C.; Caccamo, C.; Pizzimenti. G. J . Chem. Phys. 1983,
(43) Winget, D. E.; Hansen, C. J.; Liebert, J.; Van Horn, H.M.; Fontaine, G . ;Nather, R. E.; Kepler, S. 0.; Lamb, D.Q. Asrrophys. J . Letr. 1987, 315,
Charged hard spheres are a simple (albeit unrealistic) model for molten salts. They interact with the usual hard-sphere potential plus a Coulombic term for r C (u, + U b ) / 2 uab(r) = = Z,Zbe2/r
78, 357.
for r
> (a,
+ (Tb)/2
(5.1)
177.
5218
The Journal of Physical Chemistry, Vol. 94, No. 13, 1990
dwarf luminosity function. From these luminosity calculations, the age of the universe is estimated to be a relatively young (10-12) x lo9 years. 6. Polydisperse Hard Spheres The density functional theory described above has been generalized to study the freezing of a polydisperse liquid of hard spheres into both face centered cubic (fcc) and hexagonally close packed (hcp) crystal^.^**^^ Two, physically relevant, continuous distributions of particle size have been studied, the gamma (or Schulz) distribution and the Gaussian distribution. The structure of a liquid of polydisperse hard spheres can be calculated analytically, and quite accurately, from the approximate PercusYevick integral equation. For both distributions, it is found that, when the standard deviation of the particle size distribution exceeds approximately 5% of the mean size, the liquid no longer freezes into a crystalline array. Despite the approximations involved in the interactions between the particles in our model, this result is in relatively good agreement with experiments on real colloidal suspensions. The statistical mechanics of liquids with a continuous distribution of sizes has been studied carefully by O n ~ a g e rBlum, ,~~ Stell, and S a l a c u ~ e , 4and ~ ~ other^.^'-^^ ~~ Independently, Barrat and HansenB and McRae and HaymetZ8have studied the ordering of polydisperse hard spheres into crystalline arrays. Experimentally, colloidal suspensions of polystyrene spheres, and other materials,50are observed to form crystalline arrays as the thermodynamic conditions (temperature, electrolyte concentration) are changed. Although most colloids are modeled via the screened Coulomb potentialSo(which is not obviously corrects1),there are some colloids for which the hard-sphere interaction is truly app r ~ p r i a t e . ~Moreover, ~ , ~ ~ the hard-sphere model is also a useful approximation for developing theoretical techniques. Some earlier ~~ of models have used t r i a n g ~ l a r ~or” ~r e~ c t a n g ~ l a rdistributions particle size, even though the gamma and Gaussian distributions used by McRae and Haymet28are the experimentally relevant ones.s7 For these distributions, the results of Blum and Ste114s may be used to obtain analytic expressions for the direct correlation function of the liquid, which of course plays a crucial role in the DF freezing t h e ~ r y . ~ A system of polydisperse hard spheres has been studiedz8 with average diameter uo and standard deviation in particle size s, with the diameters distributed according to both the gamma distribution
fig) = ~ n a c z - le-Xu / U a )
(6.1)
where uo = a/A, s = ( c u / A ~ ) ’ / ~ r(a) , is the usual gamma function, and the Gaussian distribution Throughout this discussion, the average particle size is fixed, and the distance unit is such that the average diameter is unity. Strictly speaking, the Gaussian distribution is unphysical because it ascribes (a tiny) nonzero probability to particles with negative diameter. However, it is found that the upper limit of the standard deviation, beyond which the liquid will not crystallize, is small enough that this “negative tail” on the Gaussian distribution can be safely ignored. In fact, in the physically interesting region, (44) (45) 2212. (46) (47) (48) (49) 2197.
Onsager, L. Ann. N . Y . Acad. Sei. 1949, 51, 627. Blum, L ; Stell, G. J . Chem. Phys. 1979, 71, 42; erratum 1979, 72, Salacuse, J. L.; Stell, G.J . Chem. Phys. 1982, 77, 3714. Vrij, A. J . Chem. Phys. 1979, 71, 3267. van Beurten, P.; Vrij, A. J . Chem. Phys. 1981, 74, 2744. Griffith. W. L.; Triolo, R.; Compere, A . L. Phys. Rev. A 1986, 33,
(50) Pieranski. P. Contemp. Phys. 1983, 24, 25. ( 5 1 ) Vlachy, V.; Haymet. A. D . J . J . Chem. Phys. 1986, 8 4 , 5874. (52) Pusey, P. N . J . Phys. 1987, 48, 709. (53) Pusey, P. N.; van Megen, W. Nature 1986, 320, 340. (54) Dickinson, E.; Parker, R . J . Phys. Lett. 1985, 46, L229. ( 5 5 ) Dickinson, E.; Parker, R.; La], M . Chem. Phys. Lett. 1981, 79, 578. (56) Dickinson, E. Chem. Phys. Lett. 1978. 57, 148. (57) Hunter. R . J. Foundations of Colloid Science; Clarendon: Oxford, 1987
Rick and Haymet for fixed standard deviation, the two distributions are barely distinguishable. The polydisperse limit of the freezing functional (2.1 1 ) is just the analogue of the limit which takes a function of v variables into a function. All microscopic parameters in the liquid mixture which depend on a component index i are mapped into continuous analogues. For a liquid of N , polydisperse hard spheres, the number fraction N J N , of spheres with size ai is mapped into the probability density function (pdf)f(u) da. This quantity is just the probability of finding a sphere of diameter a in the liquid. Sums over discrete liquid densities in the functional (2.1 1 ) are mapped into integrals over the particle size pdf
The limit of sums over discrete partial crystal densities is more tricky to perform, due to the position dependence of the partial crystal density p j ( r ) . This gives rise to the wonderful complexity of polydisperse freezing. We define a position- and diameterdependent function p * ( r , a )
Plausible first approximations are introduced28for this crystal density. This set of approximations is appropriate for liquids which are not too polydisperse. First, only freezing into substitutionally disordered crystals is considered. Second, the shape of the crystal density about any lattice site is assumed to be independent of the size of the particle occupying it. The crystal density may then be written p*(r,o) = p ( r )r”(4 (6.5) where p ( r ) is the total number density at position r in the crystal and /s(u) is the hard sphere diameter pdf in the crystal phase. This is not necessarily equal to the liquid pdf the system may phase separate on freezing, as do binary mixtures of hard spheres.2w22 A simple version of the polydisperse freezing theory is obtained by making the Gaussian approximation for the crystal density and introducing the “constrained eutectic” (CE) approximation, which is physically meaningful for small polydispersities, in which the crystal pdf is the same as the liquid pdf
rs(d = P ( d (6.6) and hence ss = sL, uos = coL= I , and b = 0. The example is completed by choosing a probability distribution function for the diameters of the hard spheres in the liquid. Both gamma and Gaussian pdf (not to be confused with the Gaussian approximation for the average crystal density!) have been investigated. Choosing the Gaussian pdf, for both the polydisperse liquid and the crystal, leads to a simplified functional which has the form of the functional for an effective one-component system A/311/pLV = 1 - ( p s / p L ) [In p L j/27?2w
+ 3 In t d / * + 5/21 - 72(l + ?I2
c PL,ZW,) n
(6.7)
where P ( k ) = PLs(k) = PLL = Fss. In this example, the polydispersity enters the problem exclusively through the average structure of the liquid, defined by C(k). This liquid structure is the central input to the freezing theory and must now be determined. Note that this liquid structure is completely different from the average polydisperse liquid structure introduced by Griffith and c o - ~ o r k e r sfor ~ ~the analysis of experimental scattering data; the two types of averages are compared in ref 28. The partial direct correlation function between hard spheres of diameter uj and a,, in a polydisperse liquid with particle diameter distribution functionf(a), can be calculated analytically within the Percus-Yevick approximati~n,~~ which is known to be relatively accurate in the monodisperse limit.s9*60 By use of the methods (58) Percus, J. K.; Yevick, G. J. Phys. Reu. 1958, 110, 1 . (59) Wertheim, M . S.Phys. Reo. Lett. 1963, 10, 321. (60) Thiel, E. J . Chem. Phys. 1963, 39, 474.
Tht? Journal of Physical Chemistry, Vol. 94, No. 13, 1990 5219
Feature Article
I 0.030
I
I
0.035
0.040
I 0.045
S
Figure 7. Coexisting liquid and fcc crystal densities of the polydisperse hard-sphere system with both Gaussian (solid line) and gamma (dashed line) distributions of particle diameters, as a function of the standard deviation in particle diameter s.
of Blum and Ste11,29~4s the direct correlation function can be written in closed form. With this input, two crystal structures have been examined in detail, the face-centered-cubic (fcc) structure and the closely related hexagonal-close-packed (hcp) structure. Stable, substitutionally disordered crystals of both types are found to coexist with polydisperse liquid. Metastable crystals with the more open body-centered-cubic (bcc) structure were sought, but none were found. Since the bcc structure is unstable for monodisperse hard spheres, this is not a surprising result. A sample polydisperse hard-sphere phase diagram, in the CE approximation, is shown in Figure 7,which displays the coexisting liquid and crystal densities as a function of polydispersity s, where s is the standard deviation in the distribution of particle diameters. The data in Figure 7 is for the freezing of a polydisperse liquid with both Gaussian pdf (solid line) and gamma pdf (dashed line) into an fcc crystal. The density at which the liquid freezes increases markedly as the polydispersity s increases and the fractional density change on freezing generally decreases. At a critical polydispersity, here s = 0.045, the fcc crystal is no longer stable; a general feature of all such DF calculations is that, when the polydispersity exceeds about 5%, the liquid no longer freezes into a periodic array. Presumably this means that the high-density limit of even more polydisperse systems is some kind of “glass” or amorphous solid. This appears to be in fair agreement with experiments on real polydisperse systems.52 An important additional prediction from these calculations is that the liquid freezing density increases measurably as the polydispersity is increased, before the “freezing limit” is reached. In Figure 7, the fcc phase diagrams for a polydisperse liquid with Gaussian and gamma particle diameter distributions (pdf) are compared. In the monodisperse limit, these pdfs are of course identical. In the physically interesting region for freezing, these distributions remain almost the same. At the largest polydispersity shown there are discernible differences in the average direct correlation function transforms, but these differences have a minimal effect on the freezing predictions. At higher polydispersities, the gamma distribution becomes asymmetric about the average particle diameter, but for hard spheres this is beyond the predicted critical freezing limit. These calculations are an interesting first step in understanding the crystallization of polydisperse liquids. However, there remains the open question of characterizing the solid (amorphous) phase of more polydisperse liquids.
7. Conclusions The above calculations, by a number of different groups, indicate that the density functional theory of freezing has made substantial progress in the long-standing question of liquid/crystal (61) Madden, W. G.; Rice, S . A. J . Chem. Phys. 1980, 72, 4208. (62)Singh, Y.;Stoessel, J. P.;Wolynes, P. G. Phys. Reo. Lett. 1985, 54, 1059.
phase coexistence. It is perhaps now appropriate to summarize the shortcomings and prospects for further development. Since the density functional theory is by its very construction a mean-field theory, it suffers the same problems as any such theory when fluctuation effects become dominant. For example, Evans1 has summarized the problems in the analogous gas-liquid theory near the critical point. In two dimensions, the DF theory in its simple form does not predict the algebraic decay of correlation functions in two-dimensional crystals (although it still has some practical uses). In the isotropic/nematic transition in liquid crystals, the DF theory, along with any mean-field theory, cannot predict the small entropy change at this very “weak” first-order transition. It is important to remember that these interesting examples do not negate the central problem of describing the “usual”, “strong” first-order transitions found in most materials. At the heart of the DF theory is the truncation of a Taylor series expansion of the (grand) free energy functional. The vexing question of analyticity of the functional at a first-order transition remains. It is fair to say that most workers agree that the fact that there is a weak singularity does not prevent the application of the theory to practical problems. At least one nonperturbative method has been pr~posed,~’ although any such approach involves much more computational work. The optimization of the functional is an extremely active and promising area of research. In practical applications, one is forced (usually by ignorance of higher order correlation functions in the liquid) to truncate the above-mentioned Taylor expansion at second order. Given this constraint, at least two independent methods for improving the DF theory are available. Very recently, it has been realized that the ideal system, which underlies the derivation of the DF theory, can be chosen to minimize higher order contributions. In all the work presented here, the ideal system was chosen (almost by default) to be the usual, classical ideal gas of point particles. This is a sensible choice, since it implies that the liquid correlation function c ( r ) in the functional is the usual Orstein-Zernike direct correlation function. But recent work1*J9 shows that other choices of the ideal system are possible and even preferable. The initial calculations have been done for quantum freezing, for example for helium, but the applications to freezing of polymers63and other complex liquids are immediate. In simple terms, the procedure is simply to choose the ideal system to minimize higher order contributions to the functional. The second approach, which has become popular very recently, is the approximate inclusion of higher order terms in the functional. Most recently, at least three have proposed different “infinite order” functionals, which at the present time appear to possess very different properties. Originally, Tarazona8 used “smoothed density approximations” to build into the functional partial inclusion of higher order terms, involving higher order liquid correlation functions.67 This approach is related to an earlier fine grained “generalized van der Waals theory” of Nordholm and co-workers6*and has been generalized and improved in important calculations by several g r o ~ p s . ’ In ~ the , ~ present ~ ~ ~ article, ~ ~ ~ the ~ “smoothing” or “weighting” function never appears (equivalently, it is set equal to unity). Other, empirical or semiempirical choices can be made, for example, by choosing a polynomial weighting function such that for hard spheres either the Carnahan-Starling equation is recovered or thermodynamic consistency is enforced. The great value of this approach has been demonstrated by Tarazona? who has presented a functional which displays three-phase coexistence at the triple point. The activity in the above extensions of DF theory indicates that the renaissance of density functional theory in statistical mechanics, especially phase transitions, is far from complete. Unsolved (63) McMullen, W. E.;Freed, K. F. J . Chem. Phys. 1990, 92, 1413. (64) Rosenfeld, Y. Unpublished work. (65) Denton, A.; Ashcroft, W. Unpublished work. (66) Lutsko, J. F.; Baus, M. Unpublished work. (67) Curtin, W. A. Phys. Reo. B 1989, 39, 6775. (68)Nordholm, S.;Johnson, M.; Freasier, B. C. AUSI.J . Chem. 1980, 33, 2139. Nordholm, S . ; Haymet, A. D. J . Aust. J . Chem. 1980, 33, 2013. (69) Meister, T.F.; Kroll, D. M. Phys. Reo. A 1985, 31, 4055. (70) Igloi, F.; Hafner, J . J . Phys. C 1986,19, 5799.
J. Phys. Chem. 1990, 94. 5220-5221
5220
problems, and opportunities for extension beyond mean-field theory, remain. Acknowledgment. All of our research on freezing has been supported by grants from the donors of the Petroleum Research
Fund, administered by the American Chemical Society, and from NSF and Ford Research through the PYI program. We also acknowledge many helpful discussions over many years with Professor David Oxtoby, and Drs. Shep Smithline, John McCoy, Brian Laird, Robin McRae, and Chris Marshall.
ARTICLES Structures of Dhitrosobenrenes and 1,2-DinitrosoethyIenes Stephen Marriott, Ronald D. Topsom,* Department of Chemistry, La Trobe University, Bundoora, Australia 3083
and Brian G. Gowenlock Department of Chemistry, Heriot- Watt University, Edinburgh, EHI 4 4AS Great Britain (Received: August 14, 1989: In Final Form: January 19, 1990)
The structures of 0 - , m-, and p-dinitrosobenzenes have been studied at the ab initio STO-3G and 3-21G basis levels. There is no evidence for any additional effects compared to nitrosobenzene as a model. A study of 1,2-dinitrosoethylenes gave similar results, contrary to a report by earlier workers.
Introduction We were most interested in a recent work’ on the structures and charge distributions in 1,2-dinitrosoethylenes. Studies using ab initio molecular orbital calculations with various basis sets reported’ that the most stable structures for both the cis and trans isomers had carbon-carbon bonds of 1.43-1 $50A, compared to 1.32 A in the parent nitrosoethylene. Our interest was 2-fold. First, as far as we are aware,* this was the first report of significant lengthening of a double bond when conjugated to two apparently x-electron-withdrawing groups. We were interested to see if the effect was also found in 0-and p-dinitrosobenzenes. In earlier work,z3we had found that “through conjugation” by one r-electron-donating and one n-electron-accepting substituent had much less effect on geometry in parasubstituted benzenes than in 1,2-disubstituted ethylenes. Indeed, the effect was not great in either case and was insignificant in the benzenes, unless the two groups were powerful in n-donation and n-acceptance, respectively. Second, we have previously4” observed unusual behavior in para-substituted nitrosobenzenes (1). Thus, studies4 of the intensities of ring-breathing frequencies in the infrared region have indicated that the *-withdrawing behavior of the nitroso group could be modified by the para substituent and, at an extreme, the nitroso group could become n-donating. In other ~ t u d i e sit, ~was found that the carbon-I 3 substituent chemical shifts at carbon 1 were anomalous in that they shifted upfield, regardless of (1 1 Sedano, E.; Sarasola, C.; Ugaide, J. M; Irazabalbeitia, 1. X.;Guerro, A. G. J. f h y s . Chem. 1988, 92, 5094. Sedano, E.; Sarasola, C.; Ugalde, J. M. Tetrahedron 1989, 45, 6531. (2) Topsom, R. D. frog. f h y s . Org. Chem. 1987, 16, 85. (3) Marriott, S.; Topsom. R. D. J . Mol. Srruct. (THEOCHEM) 1984, 110, 3~. 37. (4) Katritzky, A. R.; Keogh, H. J.; Ohlenrott, S.; Topsom, R. D. J. A m . Chem. SOC.1970, 92, 6855. ( 5 ) AI-Tahou, 8 . M.; Gowenlock, B. G. Red. Trau. Chim. fays-Bas 1986, /OS. 253. . .., -. . (6) Butt, G.; Topsom. R. D.:Gowenlock, B. G.; Hunter, J . A. J. M o l . Srruct. (THEOCHEM) 1988, 164. 399.
NO
I / c ’
Q \
x (1)
whether X was a a-donor or n-acceptor. The upfield shift for the nitro substituent, for example, was comparable with that for the NMez substituents. This contrasts with the behavior of any other series of substituted benzenes so far studied. The general finding is that electron-withdrawing substituents at position 4 cause downfield shifts and electron-donating ones lead to upfield shifts. We examined6 the atomic electron populations at C1 in the para-substituted nitrosobenzenes at the STO-3G basis, and these showed the generally expected increase or decrease with substituent. However, standard geometries were employed. We wished to see if changes in structure were responsible for these effects. Accordingly, we examined the preferred structure of 0-,’m-, and p-dinitrosobenzenes at both the STO-3G and 3-21G basis levels. We also reexamined the previous work on cis- and frans-1,2-dinitrosoethylenes at the 4-3 IG and 4-3 lG* and 631G** basis levels. Calculations and Results All calculations were done at the ab initio STO-3G, 3-21G, 4-31G, 4-31G*, or 6-31G** bases, using the GAUSSIAN-86 series of programs.* Geometry optimization carried out was for all bond ( 7 ) I t is, of course, recognized that the o-dinitrosobenzene occurs only as a transient species (see, e.g.: Mallory, F. B.; Manatt. S.L.; Woods, C. S.J. A m . Chem. SOC.1965, 87, 5433. Mallory, F. B.; Cammarata, A. J . Am. Chem. SOC.1966, 88, 61. Hoffmann, R.;Gleiter, R.; Mallory, F. B. J. Am. Chew Soc 1970. 92. 1460) in the dissociation of benzofurazan oxide
0022-3654/90/2094-5220.$02.50/0 @ 1990 American Chemical Society
