6292
J. Phys. Chem. 1993,97, 62926303
Equilibrium, Stability, and Density Anomalies in a Lattice Model with Core-Softening and Directional Bonding Steven S. Borick and Pablo G. Debenedetti’ Department of Chemical Engineering, Princeton University, Princeton, New Jersey 08544-5263 Received: February 3, 1993; In Final Form: March 23, I993
The cluster-variation method has been used to investigate the phase behavior, thermal expansion characteristics, and stability limits of the Bell-Meijer lattice model; the model includes directional bonding and second neighbor repulsion. Both of these mechanisms are in principle capable of stabilizing open ground states and can therefore give rise to negative thermal expansion. A canonical implementation of the natural iteration solution technique facilitated the calculation of metastable states and stability limits. The four types of phase behavior previously found by Van Royen and Meijer34 have been mapped as a function of the relative strength of the second neighbor repulsion and the bonding interaction over a wide range of parameter space. One of the types of phase behavior is water-like: it exhibits expansion upon freezing. As the strength of the bonding interaction increases, so does the range of second neighbor repulsion within which expansion upon freezing occurs. Density anomalies appear as a low-temperature locus of density minima that merges with a higher-temperature locus of density maxima. Of the two, density minima are the more prevalent feature, always spanning a greater range of pressure. Density maxima disappear at low enough values of the second neighbor repulsion. The peculiar features of thermal expansion in this model bear little relation to the known behavior in water, heavy water, and silica, and they persist throughout the wide range of parameter space explored. An undesirable feature of the Bell-Meijer model is that bonding only contributes to long-range orientational order and not to ground-state stability.
1. Introduction In this paper we investigate the phase behavior and stability limits of a model fluid that expands when cooled at constant pressure. The model incorporates two features (directional bonding and spherically symmetric softening of the repulsivecore) that can cause negative thermal expansion. A major goal of this study is to investigatehow the relative importance of coresoftening and directional bonding influences the shape and location of the line of density extrema and how density anomalies affect the stability limits of the model’s fluid phases. Supercooled liquids will either crystallize or vitrify, depending on the rate at which they are cooled. In most cases, the thermophysical properties of supercooled liquids show no anomalies or peculiarities. There are however, important exceptions. Anomalous increases in the isobaric heat capacity’and isothermal compressibility2 have been reported for supercooled water. The isobaric heat capacities of cesium, potassium, and mixtures of potassium/cesium and sodium/potassium all increase anomalously near free~ing.~Microscopically, this increase of the response functions corresponds to growth of density and enthalpy fluctuations. Understanding the behavior of fluctuation-dominatedsupercooled liquids is a major challenge. Speedy4ss has interpreted the behavior of water by postulating the existence of a continuous spinodal curve (locus of limits of stability) bounding the superheated, stretched (negative pressure), and supercooled states. Both the isothermal compressibility and the isobaric heat capacity must diverge at the spinodal, should one exist. Power-law extrapolations of these quantities for supercooled water are not inconsistent with Speedy’s conjecture. In the P-T projection of Speedy’s picture, water’s negatively-sloped locus of density maxima intersects the spinodal at two points; the two points correspond to the liquid’s condition of maximum tensile strength and to the low-temperature/high-pressureend of the supercooled liquid spinodal. The former intersection causes the P-Tprojection of the spinodal to change slope. Recently, Poole et a1.6 have proposed a different picture. In an extensive study of ST2 water7 by molecular dynamics, these authors found no evidence of a retracing spinodal. Instead, they proposed that a novel, metastable
critical point is responsible for water’s anomalous behavior. This critical point is associated with the first-order transition between two metastable amorphous forms of ice.819 Poole et al.6 observed that the P-T projection of the density maxima locus became infinitely sloped, then positively sloped, and did not intersect the spinodal. The absence of a retracing spinodal is a major difference between Speedy’s and Poole et al.’s interpretation of the properties of supercooled water. The general question of the stability of supercooled liquids was investigated by Debenedetti and DAntoniolOJl using thermodynamic consistency arguments. These authors showed that Speedy’s picture can apply to any liquid that expands when cooled isobarically. They did not consider the possibility that the locus of density maxima could become infinitely sloped, and thus their analysisdid not address the behavior observed by Poole et al. The close connection between density anomalies and stability is a key conclusionof the Debenedetti-DAntonio analysis. Accordingly, Debenedetti and co-workers then examined the possible molecular mechanisms that can lead to density anomalies. For additive and spherically symmetric forces, Debenedetti et al.12 showed that a necessary condition for negative thermal expansion is a decrease in the repulsive force between two molecules as they approach each other. This translates into an inflection point in the repulsive core of the pair potential. Such behavior is known as core-softening. The effective pair potential of many liquid metals is core-softened.’) Stillinger and Weber14 reported both negative thermal expansion and expansion upon freezing in their molecular dynamics study of a model core-softened potential: the Gaussian core fluid. In models with hard cores, softening has been introduced as a sudden collapse of the core. Such systems have been studied extensively by Stell and co-workers,ls20 using perturbation theory, and by Young and Alder,21 using molecular dynamics. Density anomalies were reported in both cases.16921 Debenedetti et al.12 modeled core-softening in a lattice model through competition between first neighbor attraction and second neighbor repulsion. Using a first-order (quasichemical) approximation in their calculations, they found both a negatively-sloped locus of density maxima and a retracing spinodal. Negative thermal expansion in this model was caused by the appearance
0022-3654/93/2097-6292%04.00~0 0 1993 American Chemical Society
Equilibrium, Stability, and Density Anomalies of open (non-close-packed) ground states stabilized by second neighbor repulsion. The thermal or mechanical disruption of these open ground states causes the model fluid to collapse into denser and more disordered configurations. Other studies of lattice models with competing first and second neighbor interactions have focused mainly on ground-state degeneracy and the resulting phase behavior (e.g., see refs 22-25). For the present purposes, competitive interactions are important because they stabilize low-density ground states and cause density anomalies: they are a convenient representation of core-softening,the simplest spherically symmetric mechanism by which a fluid can attain negative thermal expansion. In non-spherically-symmetricsystems, density anomalies can be caused by a different mechanism: tetrahedral coordination. In both water and silicon dioxide, open structures are formed by tetrahedrally-coordinated networks stabilized by directional bonding. Various lattice-based models have been proposed to explain the propertiesof water. Realistic ones include molecular orientation in some way. proposed a model with molecules centered on the sites of a body-centered cubic lattice; each molecule could be oriented in 12 ways. In addition to a favorable bonding energy, the model included an attractive first neighbor interaction and a three-body repulsion, which discouraged close packing. Bell solved the problem in the quasichemicalapproximation,and he observed density maxima. Bell and Salt2' used the zerothorder (Bragg-Williams) approximation to solve theoriginal model of Be11.26 Density anomaliesdisappeared in this case, but explicit inclusion of long-range order allowed the formation of two icelike solid phases. Of particular interest was the appearance of an open solid which melted into a denser liquid. Fleming and G i b b ~ studied 2 ~ ~ ~a~model similar to Bell's,26 except that their energy interactions depended only on the relative orientation between two first neighbor molecules. Using a perturbation scheme,they observed density maxima as well as water-likevaporliquid coexistence. Sasai30solved a field theory model of water's hydrogen bond network in the mean-field lattice gas approximation. Though he obtained spinodals in both the supercooled liquid and the stretched liquid region, the merging of the two spinodals was thermodynamically inconsistent with the form of the density maxima locus in real water. Sastry et aL31 have recently reported water-likephase behavior, densitymaxima, and a retracing spinodal in a mean-field, two sublattice model. Net repulsion between first neighbor molecules on different sublattices accounts for the hydrogen bonding energy benefit. The effective bonding energy parameter is temperaturedependent, which results from the approximate summation over orientational degrees of freedom in the evaluation of the partition function. A slight variation of Bell's modelz6was proposed by Meijer et al.,32 who introduced a second neighbor repulsion instead of the three-body force used by Bell. Results for this interesting model included the appearance of density anomaliesand of a negativelysloped melting curve between the liquid and a less dense ice-like phase. The possible types of phase behavior were identified,34 but limits of stability were not studied. Furthermore, density anomalies were not explored systematically. With hopes of reproducing the phase diagram of water, a new study was begun using a more complex m0del,~3-~~ which included sublattices and order parameters. The goal of this model was to obtain a more realistic representation of solid phases. Though stability limits were calculated for the sublattice model?5 no attempt was made to investigate the relationship between density anomalies and spinodals. In this paper we reexamine in detail the model proposed by Be11,26 with Meijer et al.'s modifi~ation.3~An interesting and useful feature of the model is the independent parametrization of orientational energy and core-softening (second neighbor repulsion) effects. Our aim is not to describe water specifically but rather to understand in general how molecular orientation and core-softening both affect equilibrium, thermal expansion, and stability upon supercooling. Using a highly effective (but
C
b
Figure 1. Molecule with its four 'bonding arms" within a bcc cell. Molecules can also be located on each of the ell's eight vertices.
not yet widely known) solution to the cluster-variation method in the canonical ensemble26 we have thoroughly investigated the equilibrium and stability behavior of the model. Our results include quantitative parameter mapping of phase behavior and characterization of density anomalies and stability limits. In principle, one would expect a model incorporating both mechanisms capable of stabilizing high-volume ground states to yield valuable insights into what are currently poorly understood aspects of thermal expansion and stability in supercooled liquids. The rather limitedextent,as will beshown, to which thisis possible with this model suggests that more than the mere superposition of these two mechanisms is required to capture the subtleties of negative thermal expansion in real liquids. 2. The Model and Its Solution
The simple directed-bonding lattice model was originally developed by We review the details here. Molecules are centered on the sites of a body-centered cubic (bcc) lattice. Each molecule has four bonding directions, or "bonding arms",% of which two are positive and two are negative. As shown in Figure 1, the arms are directed to the vertices of a tetrahedron contained in a bcc cell. Sites u and b comprise a first neighbor pair, while sites b and c are second neighbors. A bond is formed only when a positive arm of one molecule and a negative arm of a first neighbor are directed toward each other. A molecule can be oriented in 12 different ways (because a cube contains 12 face diagonals), so a pair of first neighbor molecules has 144 mutual orientation possibilities. Only 18 of these lead to bonding, which can be understood as follows. For a given bonding arm pointing to a given vertex, the remaining arms can be oriented in three distinct ways. Thus, two molecules in bonding orientation can be oriented in 2(32) ways, where the 2 appears because a bond is still formed if the polarity of the bonded arms is reversed. The energy interactions we consider include a first neighbor energy Ele, a second neighbor energy E ~ cand , a bonding energy E b e that occurs when two first neighbor molecules are correctly oriented. Note that El, E2, and E b are pure numbers while c has units of energy. In previous work the second neighbor interaction has been incorporated in two different ways. Bell's treatment26 includes an effectivethreebody force in which the second neighbor energy is only active if the two second neighbor molecules also share a common first neighbor molecule. In Meijer et a1.k presentation32 the second neighbor energy is active for all second pairs. Meijer et al. used second neighbor repulsions to account for so-called pseudo-three-body forces found from dipole considerations: two second neighbor molecules sharing a common first neighbor experience a net repulsive dipole force when the three molecules are oriented such that the overall first neighbor dipole energy is a minimum. This work uses the approach of Meijer et al., but not for reasons based on three-body dipole forces. As shown by Debenedettiet al.,12 first neighbor attraction and second neighbor repulsion is the simplest way to approximate core-softening in a lattice model: the first neighbor attraction represents the softened part of the repulsive core. Hard core repulsion is included through the single occupancy of sites.
6294
The Journal of Physical Chemistry, Vol. 97, No. 23, 1993
Borick and Debenedetti
Figure 3. Construction of a basic tetrahedron within the bcc lattice.
(c> F l p e 2. (a) Open ground-stateconfigurationof molecular centers. (b) Close-packed ground-state configuration of molecular centers. (c) Directional bonding consistent with open ground-state configuration. Long-ranged orientational order of ice I(c) type. All occupied sites are equivalent;open circles highlight molecules at the body center of a cube.
An important feature of the model is that an open ground state can be stabilized by second neighbor repulsion. Figure 2a shows the configuration of molecular centers for the open ground state on the bcc lattice. All occupied sites are equivalent; open circles merely highlight molecules at the body center of a cube. No second neighbor pairs are present, and only I/.I of the lattice's first neighbor contacts are active. Denoting the number of lattice sites by Nand the number of molecules by M we have, for the case in which every bond is active,
2M=N
(1)
U = N(E,e + Ebe) = 2Mt(E, + Eb)
(2) (3)
where U and H are the configurational energy and enthalpy, P is the pressure, and vo is the volume per site. The first neighbor energy is always attractive; for reference we fix El = -1. For the close-packed structure with maximum possible bonding, shown in Figure 2b,
M=N
(4)
U = Ne(4 + 3E2 + 2EJ
(5)
where E2 > 0 is the magnitude of the second neighbor repulsive barrier relative to lEll = 1. When all possible bonds are active, there are two bonds per molecule on both the open and close packed ground states; for a given lattice size, then, there are twice as many bonds per site in the close-packed ground state. At absolute zero and fixed pressure, the ground state has the lowest enthalpy. Comparing eqs 3 and 6 (with El = -1) we see that the close-packed structure will be stable if (7) For E2 I2/3, the close-packed ground state is favored for all positive pressures. Note that the bonding energy benefit has no effect on the ground-state stability. Rather, favorable bonding energy results in long-range orientational order of the ice I(c) type, shown in Figure 2c. To obtain the equilibrium properties of the model, we use the cluster-variation method (CVM).37.38 The essence of the CVM
is as follows. A so-called basic cluster is chosen. It is representative of the entire lattice. A site on the basic cluster can be occupied or empty; each complete characterization of every site of the basic cluster (including assigning orientations to the occupied sites) comprises a basic configuration. Any configuration of the entire lattice can be generated by tessellating with an appropriate distribution of basic configurations. For a fixed number of molecules M on N lattice sites ( p = M/N), a first improvement over the Bragg-Williams approximation to the correct enumeration of lattice configurations is given by an independent distribution of basic configurations. This is the counting approximation used in the cell-based quasichemical method.12-26 Now, such an independent distribution overcounts the number of lattice configurations,because it allows overlapping faces, edges, or single sites of the basic configurations to be unequally occupied. For example, an independent distribution of basic configurations can lead to conflicting assignments for the state of a site shared by two or more basic clusters. The CVM attempts to account for overlapdiscrepanciesby considering not only the distribution of basic cluster configurations, but also the distributions of the so-called subclusters that result when basic clusters overlap. Hijmans and DeB0er3~present a thorough discussion of systematicsubcluster generation, but wecan illustrate the process for the present model. We choose a tetrahedron (denoted by w ) as the basic cluster. Figure 3 shows how a basic tetrahedron fits into the lattice. Two adjacent tetrahedra can touch along a triangular face, a second neighbor bond, a first neighbor bond, or a site. Since no new geometric figures result from the adjacent contact of any of these four figures, the subclusters for the tetrahedron basic cluster are just the triangle ( u ) , second neighbor bond ( z ) ,first neighbor bond (y), and single site ( x ) (see note in ref 40). When the tetrahedron basic cluster has been chosen and the appropriate subclusters found, the expression for the number of lattice configurations available to M molecules on N lattice sites ( p = M/N)becomes
where, in general, fr is the occurrence probability of the ith configuration of clusterf and d i s the weight factor (degeneracy) of the same configuration. Each of the probability distributions f~ (f1, fi, ...) is normalized according to
(9) Equation 8 is density dependent through normalization of the wf which, as will be shown below, depends explicitly upon density. Upper bounds on the product symbols are the number of distinguishable configurations for each cluster; these configu-
The Journal of Physical Chemistry, Vol. 97, No. 23, 1993 6295
Equilibrium, Stability, and Density Anomalies
TABLE k Configurations and Degeneracies for Basic Cluster and Subclusters. Energies for Basic Configurations basic configurations subfigure configurations no. config energy degeneracy no. config degeneracy 0;QO
1
1(1)
2 2
0; 0; 0
3
1 0; -; 0 4
4~2)
2(144)
3
4 5
-.1 0 -1
4(18)
6
4( 126)
7
6’ ’6 1
;; 0;0
8 6
-1.1-. -1 3’ 4’ 6
4(432)
1
2
-.1 -.1 0 3’ 4’
4( 1296)
3 1
-2.1-. -1 3’ 2’ 3
2(324)
2 3
1 1 -.2 _. 3’ 2’ 6
4(2268)
4 1
2.1.0 3’ 2’
1 ( 1 1016)
rations will be enumerated shortly. The exponents for each term on the right hand side are the resuls of the CVM. A detailed discussionof the CVM is beyond the scope of this paper. Hijmans and DeBoer39 have shown how the exponents can be derived systematically by considering successive cluster assemblies (comprising the basic cluster, basic cluster plus largest subcluster, basic cluster plus largest subcluster plus next largest subcluster, etc.) in which the exponents arise to ensure correct counting of all clusters. Barker” formalized the CVM as a generalization of the quasichemical method. Configurations and weight factors for the basic tetrahedron and subclusters are shown in Table I. Energies for the basic tetrahedra are also shown. Black and white circles denote filled and empty sites. Solid lines are first neighbors, dashed lines are second neighbors, and bold lines are bonded pairs. In determining the weight factors, we need to consider not only the number of ways in which a configuration can be oriented in space but also the number of orientations of the molecules consistent with the number of bonds in the configuration. The weight factors for configurationswithout bonds are easy todetermine. For example,
2
tetrahedron configuration w3 can be realized in 2 ways in space, and each molecule can be oriented in 12 different directions; the weight factor is thereforeg? = 2( 122) = 288. For configurations in which bonding is possible, the situation is more complicated. Consider configuration w4 and its nonbonded analogue w5. Each can be spatially oriented in 4 ways. As shown earlier, only 18 of 144 mutual orientation possibilities for the two molecules lead to bonding. Thus, the weight factor for w4 is gr = 4(18) = 72 and for w5, gy = 4( 144 - 18) = 504. The weight factors for the remaining basic and subclusters are similarly derived. Also shown in Table I are energies for the tetrahedron configurations. The energy terms are derived as follows. Because there are six tetrahedra per lattice site, we write the lattice energy as
Now ifthelattice werecompletelyfilledandbonded(configuration
6296 The Journal of Physical Chemistry, Vol. 97, No. 23, 1993
Borick and Debenedetti
8 only; grw8 = l), its energy would be
U = 4NE,e + 3NE2e + 2NEbe
(11)
Thus
6Ne8 = 4NE1e
+ 3NE2e + 2NEbe
(12)
so
= (2/3)E1€+ (1/2)E2€+ (1/3)Ebe
(13) In configuration 8, the basic cell has 4 first neighbor contacts, 2 second neighbor contacts, and 2 active bonds. To get the correct value for ea, we need todividethe number of first neighbor contacts by 6, the number of second neighbor contactsby 4, and thenumber of active bonds by 6. If we perform the same exercise with any other basic configuration, we find that the same divisions by 6 , 4, and 6 are necessary. The subcluster probability variables are related to the basic cluster variablesby the consistency relations (subclustermaterial balances) e8
where odw is the total number of subclusters of type f contained is the number of times the jth in a tetrahedron and configuration of subcluster f occurs in the ith configuration of the basic cluster w. Equation 14 is a mathematical statement of the fact that subfigureoccurrence probabilities depend explicitly on the probabilities of the basic configurationsfrom which they come. Because of the consistency relations, eq 8 is a function only of the probability variables = w l r w2, ,..,w10 for the basic cluster. Applying Boltzmann's entropy formula (S = Nk In Q) to eq 8 and using Stirling's approximation, the entropy for the tetrahedral model is
4,
in
9
8
2
4
lnpt -pt. Combining eqs 10 and 15, we write the with L,(pt) dimensionless Helmholtz free energy per site as
E
A
where a* z A/(Ne), et* = etle, and T+ = kT/e. Note that a* is a function of the fractional occupancy ( p ) through the term for subcluster x (single site) because, from the normalization condition, x1 = 1 - p and x2 = p/12. Furthermore, density dependenceis implicit through the dependenceof the y o n density. For fixed temperature and density, the equilibrium state is found as the distribution of probability variables wtthat minimizes the Helmholtz energy. To find the equilibrium wt, we use the natural iteration (NI) method of Kiku~hi.~2 NI is distinguished by its remarkable effectiveness; the free energy decreases with each iteration, and convergence to the minimum is guaranteed. To date, most NI calculations have been performed in the grand canonicalrepresentation (first Legendre transform of Helmholtz energy), in which temperature, volume, and chemical potential are independent variables and the lattice free energy (grand potential) is given by -PV( T,V,p).32,43344 Nevertheless, we have applied NI in the canonical ensemble, using an approach similar
P Figure 4. Schematic isotherm for the Bell-Meijer model in p vs P projection. Branch 1 is the vapor phase. Branch 2 is an intermediatedensity phase (shown here as metastable) which in this model exhibits both solid-like ordered behavior at low temperatures and disordered behavior at higher temperatures. Branch 3 is the liquid phase. Point d is a limit of stability.
to that of Lawrence and Rossiter.36 Though the solution in the canonical ensemble requires slightly more programming effort, it has several advantages over the grand canonical formulation. Using temperature and fractional occupancy (density) as independent variablesensureS a single free energy extremum for given state conditions, even in the presence of metastable or unstable states. In the grand canonical formulation, it becomes difficult to locate metastable states for more complex phase behavior, because for a given value of chemical potential the metastable state may not be the state of lowest grand potential. This is an important limitation for this work, since we are specifically interested in metastable states. In Figure 4 a schematic p-P isotherm projection shows stable equilibrium at point 6 and a metastablestateat point c. Themetastablestateccanbe reached by starting at point a on branch 1 (easily found by fixing p = p' in the grand canonical ensemble) and stepping in p to p". For each new value of p the initial probability distribution of basic configurations fed to the iteration scheme is the equilibrium distribution from the previousvalue. Aslong as theinitial guesses are chosen in this way, the calculation stays on the metastable segment bc until a limit of stability is at which point the calculation jumps to branch 3. Near a spinodal point, the number of iterations required for convergence to equilibrium increasesdramatically. Using the same type of reasoning, we see that there is no convenient way within the grand canonical implementation to calculate the metastable branch 2. Furthermore, Figure 4 illustrates another problem that we have encounteredin the grand canonical calculations. If the chemical potential step interval is much larger than the chemical potential change over the metastable second bd, the calculation can jump from branch 1 to branch 3 without any indicationof the spinodal point d or the equilibrium intersection b. This problem is particularly troublesome for isotherms near the critical temperature. In the canonical implementation of NI it is easy to generate entire p-P isotherm projections, with no missing sections or branches. Another advantage of the canonical implementation is that isochores and isobars are easily computed. To generate isochores,we simply fix p and calculate pressure (as shown below) for a range of P.The calculation of isobars is outlined in the Appendix. Finally, we have found that the iterations generally converge much faster in the canonical ensemble, especially near limits of stability. Given T1 and p, we minimize eq 16 with respect to the wt, subject to two constraints. The first constraint is normalization of the wt:
The grand canonical formulation uses this constraint as well. Normalization is a natural and necessary constraint, since (w1, w2, ..., WIO)is a probability-of-occurrence distribution for the
Equilibrium, Stability, and Density Anomalies basic cell configurations. The second constraint reads
which is simply the consistencyrelation for the filled site subcluster; this constraint allows us to specify density. Introducing the Lagrange multipliers p and A*, we now define in
The Journal of Physical Chemistry, Vol. 97,No. 23, 1993 6291 minimum for thegiven density and temperature. The convergence criterion is the change in the y rather than the change in a*, because the wt are the direct results of the iterations. We have observed cases where a* was relatively constant while one or more of the wt was still changing significantly. In the Appendix we show that, at equilibrium, a* =
+ 4pX* - T+ [p In(;)
+ (1 - p ) In( 1 - p ) ]
(26)
IO
and In differentiating eq 19 with respect to wIrwe note two things. First, using the consistency relations for subcluster f we have
Second, the terms for subcluster x drop out upon differentiation, since x2 = p/12 and p is constant. Taking (dii*/aWt),,p = 0 and solving for WI, we obtain the iteration equation
with
Given current values of wt, p , and A* we use eq 21 to determine a new set @,ofbasic cluster probabilities. When first beginning a calculation, it is sufficient to use infinite temperature values for the wt and zero for both p and A* as initial guesses. For subsequent state points, it is efficient and most convenient to use the equilibrium values from the previous state point as initial guesses. The iteration eq 21 is similar to that obtained in the grand canonical formulation. A most important difference, however, is that the grand canonicalequation does not contain the variable A*, which appears when we include the fixed composition constraint. Prior to each major iteration step in which we calculate the new set GI,a minor iteration or Newton-Raphson step is required to determine the values of p and A* that satisfy the normalization and density constraints. To this end, we proceed as follows. Substituting the expression for f i t into the normalization and density constraints gives
and
where x exp[p/(6T+)]. It is nearly always possible to solve the system of equations 23 and 24 for x and A* using the NewtonRaphson method. Sometimesa minor iteration step is necessary to obtain good guesses for 'f and A*, after which the Newton routine converges quickly. For p values less than about 0.999, we have always been able to successfully solve the simultaneous equations using iterations, Newton-Raphson, or a combination of both. The major iterations continue until the convergence criterion
is satisfied. At this point the system's Helmholtz energy is a
Pin[ w 1 - P )]
P*=P/c=~X*+
(27)
Finally, a first Legendretransform of the Helmholtz energygives the pressure
where uo is the volume per site.
3. Results and Discussion To find phase equilibrium, we examined p* versus P* projections of isotherms,in which binodal points are easily found as intersections of different branches. As explained above, such isotherms were generated with density and temperature as independent variables. With the canonical implementation of NI, we were often able to generate completeisotherms,including unstable portions of van der Waals loops. At much lower temperatures,the cluster-variation free energy sometimesexhibits pathological behavior, characterized by discontinuous jumps in free energy and probability variables. This occurs at limit points where isotherms have a singular slope. Pathological behavior in the CVM is attributable to defects in the approximation,not the solution technique; Kiku~hi3~ and Barkefll have documented related behavior for other lattice systems solved by different techniques. In most cases we found that loss of solutions occurred within the unstable region of an isotherm. Thus, the usefulness of the CVM was not compromised. Though we used the canonical implementation of NI for most work, it was necessary, for numerical reasons, to use the grand canonical implementation when vapor and liquid coexisted at extremes of density. In many such cases, an intermediate-density phase also existed in either stable or metastable equilibrium with the vapor and liquid, as exemplified in Figure 4. As explained earlier, without a good initial guess of probability variables it is impossible to calculate branch 2 of Figure 4 in the grand canonical implementation. Nevertheless, since branch 2 is of moderate density it can easily be calculated in the canonical ensemble. We use the canonical ensemble to obtain a biased initial guess that ensuresconvergence to branch 2 in the grand canonical ensemble. Figure 5 illustrates the model's four different types of phase behavior and the region of parameter space where each typeoccurs. The four types are as follows: vapor-liquid equilibrium (type 1 behavior), one triple point and one critical point with the appearance of an ice-like ordered phase of density intermediate between that of the two fluid phases (type 2 behavior), two triple points and two critical points with the appearance of a dense fluid phase (type 3 behavior), and two separate fluid-fluid transitions, each with a critical point (type 4 behavior). When E2 < 2/), simple vapor-liquid equilibrium (type 1 behavior, shown in Figure 6 ) occurs for all Eb less than zero. Meijer et al.32 reported this result as well. As E2 is increased, a third phase appears at lower temperatures (type 2 behavior), as shown in Figure 7. This third phase has a density that is intermediate between that of the liquid and thevapor; hence, we refer toit as theintermediate phase. We shall see that, depending on the parameter values and the temperature, the intermediate phase can behave like a liquid or a solid. The solid-like behavior is interesting, considering that we have not attempted to impose long-range order in the model.
Borick and Debenedetti
6298 The Journal of Physical Chemistry, Vol. 97, No. 23, 1993
,901
I
E,
= 0.72,
E,
-4
=
i1
I
V
t
1
,,II Y A k -
,65
-4
-5
-2
-3
-1
0
E b Figure 5. Parameter mapping of the model's phase behavior. E2 and Eb
are the second neighbor repulsion and bonding energies, respectively. (El = first neighbor attraction = -1). 1.8
i
j 5t
E, = 0 . 6 5 , E, = -5
I S
01
1,2i ,
1 .oo
,
.2
.4 . 6
.2
0
. 8 1.0
P Figure 7. Binodals for type 2 phase behavior. (a) P-T projection (b) T-p projection. s, I, and v denote the intermediate, liquid, and vapor
phases, respectively. ,
,
,
,
.6
.4
,
,
.8
,\,
1
P Figure 6. T-p binodal for type 1 phase behavior.
++ In general order parameters, which indicate preferential occupation of sublattices, must be included to generate long-range order in lattice models describing solid phase^.^'.^^ The present model includes neither order parameters nor sublattices. For some values of the parameters and low enough temperatures, however, all of the basic cluster probabilities for the intermediate phase become negligibly small, with the exception of w4, which is the only basic configuration found in the ice I(c) structure (see Table I and Figure 2). The cluster variation method, then, closely approximates the appearance of long-range order without the explicit introduction of order parameters. Solid-likebehavior of the intermediate phase is possible because, for E2 > 2/3, an open ground state can be stabilized at positive pressures, as discussed in section 2. In Figure 7 the value of E2 is 0.72,so from eq 7 the open ground state is stable below P* = 0.16. This is not inconsistent with the extrapolation to T* = 0 of the solid-liquid binodal in Figure 7a. Figure 7b shows that, below the triple point, the density of the liquid in equilibrium with the intermediate phase is very close to 1. For this part of the phase diagram we use the grand canonical formulation. The density difference between the intermediate phase (here solidlike) coexisting with the vapor and with the liquid is too small to be seen on the figure; in T-p, the intermediate phase appears as single line, labeled s. Initial probability distributions for the intermediate phase isotherm are generated in the canonical ensemble, as described previously. Further increase of E2 results in the phase behavior shown in Figure 8 (type 3 behavior). At higher temperatures there is a triple point involving only fluid phases. Stell and co-workers have reported such triple points in their studies of repulsive core
a
I -.I'
.1
V '
"
.2
'
"
"
.4 . 5
.3
'
.6
T* .6,
E,
i
'
'0
'
.'2
'
.'4
E,
= 0.83,
'
.'6
'
= -1
~
I
!
.'8 1 0 '
P Figure 8. Binodals for type 3 phase behavior. (a) P-T projection (b) T-p projection. dl denotes the dense liquid phase.
collapse.1~19Debenedettiet a1.12 have reported the same behavior in cluster quasichemical calculations on a lattice model with attractive first neighbor and repulsive second neighbor interactions. Phaseequilibrium above the high-temperature triple point
Equilibrium, Stability, and Density Anomalies
The Journal of Physical Chemistry, Vol. 97,No. 23, 1993 6299 58
1
E,
= 0.83,
E, = -1
A
,j
~
~
*
b .
.20
.2
0
.6
.4
.8
1.0
42i5'
io'
2 5 ' 3 0 ' 3 5 ' 40' A5
T*
P Figure 9. T-p binodal for type 4 phase behavior.
Figure 11. Behavior of isobars showing density extrema. The label on each curve is the dimensionless pressure P* Pvo/z.
comprehensive calculations aimed at mapping the dependence of density anomalies (and of their relationship to phase boundaries and stability limits) upon model parameters. Figure 11 shows isobars for a potential within the type 3 phase behavior region. For pressures above P* = 4 . 0 3 , the isobars exhibit a minimum at lower temperatures and a maximum at higher temperatures. A smooth curve can be drawn through the extrema, generating E, = -5 I( a high-temperature locus of density maxima which merges with a low-temperature locus of density minima at T+ = 0.3 on the P* = 4 . 0 3 isobar. The high-temperature end of the density maxima locus terminates at a limit of stability (not shown) on the P* = 0.01 isobar. Debenedetti and D'Antonio" have used thermodynamic consistency arguments to show that a density maxima locus can terminate in only two ways: by intersecting a spinodal line or by merging with a locus of density minima. At the intersection of a locus of density maxima with a spinodal curve, the slope of the spinodal changes sign in the P-Tprojection. reported simultaneous density maxima and minima in his quasichemicalstudy;we find, in fact, that this is the distinguishing feature of density anomalies in the directed-bonding model. We now examine the behavior of the density extrema locus and its relationship to stability limits. Numerical details on the E, = 0 . 6 5 , E, = 0 calculation of stability limits in the cluster variation method are ,5/r------, , given in the Appendix. Density anomalies and stability limits can be understood most clearly by following their evolution from one type of phase behavior to another. Within the parameter range of interest, phase behavior is most sensitive to changes in 0 0 . 2 . 4 .6 . 8 1.0 E2, as seen in Figure 5. Our discussion, therefore, will begin at higher values of E2 and progress to lower values; Eb, however, P need not remain constant as we vary E2. Figure 12a shows the Figure 10. Effect of bonding interaction on type 1 phase and stability P-T projection for a potential within the type 4 phase behavior behavior. (a) P-Tprojection (b) T-p projection. Full lines are binodals; region. The two separate fluid-fluid binodals are shown as solid dotted lines are spinodals. lines labeled b. Spinodals are labeled s. The most interesting feature of Figure 12 is the locus of density minima mn which is difficult to detect in grand canonical calculations, for reasons discussedearlier regarding metastable continuationsof isotherms. intersects a locus of density maxima mx at point c. At higher pressures the density maxima locus attains a negative slope, and By stepping in p in the canonical ensemble, however, it is easy to generate van der Waals loops and to obtain binodal points. At the locus ends upon intersection with the spinodal curve along which the fluid phase that is stable at intermediate density and lower temperatures we encounter another triple point, again between the low- and high-density fluids (vapor, dense liquid) pressure reaches its supercoolinglimit with respect to the denser and the intermediate phase. The overall phase diagram is shown fluid phase. As demanded by thermodynamic consistency, the in Figure 8. As E2 is increased within region 3, the upper and slope of the P-T projection of the spinodal changes sign upon intersecting the locus of density extrema. As shown in Figure lower triple points approach each other. When they finally merge, the phase behavior shown in Figure 9 results (type 4 behavior). 12b, the locus of density minima (and most of the locus of density Type 4 behavior persisted even up to E2 = 10. maxima) lies in the stable region of the intermediate-densityfluid; The effect of the bonding attraction ( E b ) is shown in Figure extension into the metastable region where the intermediatedensity fluid becomes supercooled with respect to the high-density 10. When E b < 0, the T-p projections become asymmetric with respect to p = 0.5, and the vapor side of the phase diagram is fluid precedes termination of the density anomaly locus upon intersectionwith a spinodal curve. Low-temperaturecontinuation enhanced. Meijer et ales2attribute the asymmetry to increased bonding character in the lower density phases. When E b = 0 the of the density anomaly locus was precluded by numerical difficulties. asymmetry disappears (Figure lob), even though the bonding feature is included in the entropy expression. As E2 is decreased further, behavior like that shown in Figure Density anomalies have been reported in previous ~ t u d i e s ~ ~ . 13 ~ ~ occurs. Here we have moved into the type 3 region, characterized by two triple points. The most important features of the directed-bonding lattice model; here, we describe more
1 E,
=
I
0.65
,
~
Borick and Debenedetti
6300 The Journal of Physical Chemistry, Vol. 97, No. 23, 1993 1.0
PI
E, -.
4
0
I,.
.1
= 1,
.2
E,
I ,
0
=
.4
.3
.5 sid
T*
2
'
4
'
'6
'
'
'8 1 ' 0
P Figure 12. Binodals,spinodals, and density anomaly loci for type 4 phase behavior. (a) P-T projection (b) 7'-p projection. Symbols are as follows: b, binodal; s, spinodal; mn,density minima locus; mx, density maxima locus; c, merging of density extrema loci.
1
.1
\
'
I
'Isid
' , I
.2
'
'
.3
'
'
.4
'
'
.5
'
1 .6
T* Figure 13. Binodals, spinodals,and density anomaly loci in P-Tprojection for type 3 phase behavior. Symbols are as follows: siv, spinodal for the limit of superheating of the intermediate phase with respect to thevapor; sid, spinodal for the limit of supercoolingof the intermediate phase with respect to the dense liquid; sdi, spinodal for the limit of superheating of the dense liquid with respect to the intermediate phase. Other symbols as in
T*N 0.23
T*m 0.2
P
P
Figure 14. Schematicpvs Pisotherms correspondingto Figure 13.Points r and m are density extrema. Other symbols as in Figure 13.
0'
-.2'
sid
Figure 12.
are summarized as follows. A negatively-sloped locus of density minima mn intersects (at point c) a locus of density maxima mx, which terminates upon intersecting the intermediate liquiddense liquid spinodal sid (limit of stability of the supercooled intermediate liquid with respect to the denser liquid phase). The density anomaly locus does not touch the intermediate liquidvapor spinodal siv and thus does not cause it to change slope toward positive pressures. The relationship between density anomalies, binodals, and spinodalscan be understood more clearly by considering Figures 13 and 14. At T1 = 0.43, both the schematic p-Pisotherm in Figure 14a and the P-Tphase diagram in Figure 13 show the two equilibrium transitions. Density anomalies have not yet appeared at T1 = 0.43. Figures 13 and 14b show that, at T* = 0.35, the intermediate phase is now
metastable; a density maximum appears in this phase at point m. Below the lower temperature triple point, at T1 = 0.23, Figure 14c shows that the intermediate phase has again become stable. A density minimum now lies in this phase at point m. As discussed earlier, it is possible for the intermediate phase to exhibit solidlike character at low temperatures. Here, however, the presence of a density minimum shows that the phase is still fluid-like. At T1 = 0.2 the density minimum has moved to the now metastable intermediate phase, as shown at point r in Figure 14d. Figure 15a shows density anomalies within type 2 phase behavior. The most obvious feature in the P-T projection is a negatively-sloped density minima locus mn which intersects a positively-sloped density maxima locus mx at point c. Note the gradual disappearance of the locus of density maxima in going first from type 4 to type 3 phase behavior and then from type 3 to type 2 (Figures 12, 13, and 15). To understand the behavior of spinodals, we refer to the T-p projection in Figure 15b. A distinguishing feature is the metastable binodal p for the phase transition between the metastable intermediate phase and the more-dense liquid. Because this transition occurs only over a small temperature range, the binodal cannot be seen in P-T on the scale of Figure 15a. The metastable phase transition is essentially the same as the intermediate liquiddense liquid transition in Figure 13, except that in Figure 13 the transition is stable. Below the metastable critical point a, the liquid loses stability with respect to the intermediate phase (spinodal sli). Spinodals sli and sil, which join at the metastable critical point a in T-p, meet at a cusp in P-T. At temperatures higher than that corresponding to point a, the intermediate phase no longer exists, and the liquid now loses stability with respect to the vapor. Thus, a true liquid-vapor spinodal exists only over a small temperature interval. In P-T, the transformation of the liquidvapor spinodal slv to spinodal siv at the temperature of point a is smooth. As in Figure 13, density anomalies occur only in the intermediate phase; unlike Figure 13, they remain metastable, even when this phase becomes stable below the triple point. A further difference is that the density maxima locus in Figure 15a never attains negative slope. In Figure 16 EZ has been decreased further, but the phase behavior is still type 2. Density maxima have disappeared completely; all that remains of density anomalies in P-T is a
The Journal of Physical Chemistry, Vol. 97, No. 23. 1993 6301
Equilibrium, Stability, and Density Anomalies
Figures 14c and 14d, in which the density anomaly point approaches the spinodal point sid as the temperature is lowered. 4. Conclusions
,
,
;
,,,;', E,
0 8 2 , Eb
i
- 4
.2
.6
.4
=,
-3
.8
1
1.0
T*
4t u .2 . 4 . 6 . 8 1.0
'0
P Figure 15. Binodals,spinodals, and density anomaly loci for type 2 phase behavior. (2) P-T projection (b) T-p projection; line p is a metastable binodal with critical point a. Symbols are as follows: sil, spinodal for thelimit of supercoolingof themetastableintermediatephasewith respect
to the dense liquid;sll, spinodal for the limit of superheatingof the dense liquid with respect to the metastable intermediate phase; slv, spinodal for the limit of superheating of the dense liquid with respect to the vapor. Other symbols as in Figures 12 and 13.
\ E, = f$>
0.72, E, = -4
I'
-If!
'
.'2 ".4
'
.6
'
.'8 ~ Y O ~ l : 2 ~ l ! 4
T' Figure 16. Further evolution of binodals, spinodals, and density anomaly loci in P-T projection for type 2 phase behavior. Symbols as in Figures 12, 13, and 15. negatively-sloped density minima locus mn. Here again, density minima occur in the intermediate phase, which only becomes stable below the triple point. For this potential, the intermediate phase exhibits solid-like character. This can be seen in the T-p binodal in Figure 7b, where the density of the intermediate phase in equilibrium with both liquid and solid is very close to 0.5. In Figure 16, the liquid-vapor spinodal slu decreases monotonically with temperature until T* = 0.77; below this, the liquid no longer loses stability with respect to the vapor, but rather with respect to the metastable intermediate phase. Curve si1 is the spinodal of the intermediate phase with respect to the liquid, while siu is the spinodal of the intermediate phase with respect to the vapor. Note that the density minima locus mn and the spinodalsi1 appear to approach an intersection at lower temperature. Such an intersection is suggested by the progression of points m and r in
The possibility that liquids capable of expanding when cooled can lose mechanical stability upon supercooling is one of the most interesting open questions in liquid-state physics. There are two known mechanisms by which negative thermal expansion can occur. One involves the formation of tetrahedrally-coordinated open structures whose thermal disruption leads to denser and more disordered configurations. Water, heavy water, and silica, all of which can expand when cooled, form such tetrahedrally-coordinated networks. Alternatively, it has been shown theoretically12 that softening of the repulsive part of the pair potential also stabilizes open ground states. This mechanism may be important in liquid metals, many of whose effective pair potentials are core-~oftened.'3 Although three liquid metals (Ga, Sb, and Bi) expand upon freezing at atmospheric pressure, negative thermal expansion in the liquid state has not been reported to date. In light of these considerations, Bell's directed-bonding lattice in the slightly modified form used by Meijer et al.,32 is especially interesting because it incorporates both directional bonding and core-softening (the latter approximated by second neighbor repulsion). In this work, we have investigated how core-softening and directional bonding influence equilibrium, thermal expansion, and stability in the Bell-Meijer model. Each of the four types of phase equilibrium found had been identified previously by Van Royen and Meijer;34we have now mapped this phase behavior as a function of the relative strengths of the second neighbor repulsion and the bonding interaction energy. The model, like water, Ga, Sb, and Bi, exhibits expansion upon "freezing". This occurs because of the formation of an intermediate-density phase that, at low temperatures, forms a tetrahedrally-coordinatedcubic ice-like network and behaves like a solid. The intermediatedensity phase appears when the secondneighbor repulsion exceeds 2/3, regardless of the strength of the bonding interaction. As the strength of the bonding interaction increases, so does the range of second neighbor repulsion in which water-like phase behavior occurs. The existence of an open ground state causes negative thermal expansion (density anomalies). These appear as a low-temperature locus of density minima with a negative slope in the P-T projection. At high values of the second neighbor repulsion, the locus of density minima merges into a high-temperature locus of density maxima with initially positive slope in the P-Tprojection. The locus of density maxima tends to disappear at low values of the second neighbor repulsion. Four aspects of density anomalies in this model are unrelated to the behavior observed in water, heavy water, and silica. These aspects are as follows: the locus of density maxima, when it exists, always merges with a locus of density minima; the locus of density minima is a prevalent feature, always spanning a much greater range of pressures than the locus of density maxima; the major portion of the locus of density maxima has positive slopein the P-Tprojection; density anomalies always occur in the intermediate-density phase. This contrast between water-like phase behavior (in certain regionsof parameter space) and peculiar thermal expansion behavior is surprising, because the latter persists throughout the wide ranges of second neighbor repulsion and directional bonding strengths explored. It is apparent that more than the mere superposition of two mechanisms capable in principle of stabilizing open ground states is required to model a negatively-sloped locus of density maxima, without density minima. In this respect, an undesirable feature of the present model is the fact that bonding plays no role in stabilizing the open ground state but merely favors orientational order. In the present model the open ground state is stabilized only by second neighbor repulsion. Work is now in progress on a model that allows the bonding benefit to contribute to the
6302 The Journal of Physical Chemistry, Vol. 97, No. 23, 1993
stabilization of an open ground state. A mean-field model including this feature31 has recently been shown to give waterlike density anomaliesand stability limits that include a retracing spinodal, in agreement with Speedy's stability limit conjecture4 and thermodynamicconsistency arguments.' The more rigorous cluster-variationcalculations on a model whose open ground state is stabilized both by core-softeningand orientational interactions should therefore clarify the extent to which water-like features appear as a result of mathematical approximations or whether they are indeed inherent to the model.
Acknowledgment. P.G.D. gratefully acknowledgesthe support of the John Simon Guggenheim Memorial Foundation (1991 Fellowship), the Camille and Henry Dreyfus Foundation (1989 Teacher-Scholar Award), and the US.Department of Energy, Office of Basic Energy Sciences, Division of Chemical Sciences (Grant DE-FG02-87ER13714). S.S.B.is grateful to Professor Paul Meijer for helpful discussions.
Borick and Debenedetti Differentiating eq 26 gives
The derivatives ( a F / a p ) , , and (aA*/ap)p,w are obtained from implicit differentiation of 6 s 23 and 24. It iseasily shown, then, that (A.11) and thus eq 27 is verified. To calculate isobars, we start by rearranging eq 28 0 = pp* - a * - p*
(A.12)
Substituting for p* and a* we get
Appendix To derive eq 26, we start by writing
[
+ (1 - p) In( 1 - p)]
T* p In($)
- P* (A. 13)
Final simplification gives which holds at equilibrium because (81*/awi)p,p= 0. Using eq 19 and evaluating eq A.1, we get 10
2
2
Since the normalization condition holds for each set of probabilities, the term in square brackets is equal to 0. Also 2
~ & x ,In x, = 1 4 m=l
h)
ln(
s>
+
O = l - p - e x p [ 6 1 n ~ + j =P* ] =f(p)
(A.14)
where we have noted that F = 6P In x . For fixed P* and P, eq A.14 is an implicit function of p through the terms p and P. To solve it, we use Newton's method with a finite difference derivative. Given the current value of p, x (see eq 24) is determined by the natural iteration method. To find spinodal points, we need to evaluate (aza*/apz)p = (ap*/ap),. At a spinodal (ap*/ap)p = 0. We start with dp*
=
(
$)p,F
dT*
+
(5)
dp
+
T,W
(1 - P)
- P) (A.3)
so finally
a* =
P + 4ph* - T* [p ln(fi) + (1 - p) In( 1 - p)]
(A.4)
so
which is eq 26. To derive p*, we note first that from eq 28
+
p* = (1 /p)(a* P*) (A.5) The following identity holds for M molecules on N lattice sites
Usingeq27ineqA.16,notingthat thep termcan bedifferentiated directly, gives
T* P
P(1-P)
($)
Comparing the previous two equations gives P*
= ($$)P
We now write
(A.7)
4g(E) awi p,r
P (A.17)
The derivatives (aX*/ap)p,, and (aX*/aw,),,p are obtained implicitly from eqs 23 and 24: To determine (awllap),. we note that, at equilibrium,
At equilibrium (aa*/awl)p,p = 0 so -
(A.18)
Equilibrium, Stability, and Density Anomalies Differentiating with respect to p at fixed T* gives
(A. 19) In this equation all quantities are known except the vector (&vi/ which can now be determined using matrix algebra. Having evaluated ( d W , / d p ) p , ( a X * / d p ) p , w , and (aX*/aw,),,p, eq A.17 gives ( a p * / a p ) p .
ap),,
References and Notes (1) Angell, C. A.; Oguni, M.; Sichina, W. J. J. Phys. Chem. 1982, 86, 998. (2) Speedy, R. J.; Angell, C. A. J. Chem. Phys. 1976, 65 (3), 851. (3) Voronel, A.; Steinburg, V.; Sverbilova, T. Phys. Lerr. 1980, 79A (2,3), 183. (4) Speedy, R. J. J. Phys. Chem. 1982.86, 986. ( 5 ) Speedy, R. J. J. Phys. Chem. 1982,86, 3002. (6) Poole, P. H.; Sciortino, F.;Essmann, U.; Stanley, H. E. Nature 1992, 360, 324. (7) Stillinger, F. H.; Rahman, A. J. Chem. Phys. 1974, 60 (4), 1545. (8) Mishima, 0.;Calvert, L. D.; Whalley, E. Nature 1984, 310, 393. (9) M,ishima, 0.;Calvert, L. D.; Whalley, E. Nature 1985, 314, 76. (IO) D Antonio, M. C.; Debenedetti, P. G. J . Chem. Phys. 1987.86 (4), 2229. (11) Debenedetti, P. G.; D'Antonio, M. C. AIChE J. 1988,34 (3), 447. (12) Debenedetti, P. G.; Raghavan, V. S.;Borick, S . S. J . Phys. Chem. 1991, 95, 4540. (13) March, N. H. Liquid Metals: Concepts and Theory; Cambridge University Press: Cambridge, 1990. (14) Stillinger, F. H.; Weber, T . A. J. Chem. Phys. 1978,68 (8), 3837. (15) Hemmer, P. C.; Stell, G. Phys. Reu. Lerr. 1970, 24, 1284. (16) Stell, G.; Hemmer, P. C. J . Chem. Phys. 1972, 56 (9), 4274. (17) Kincaid, J. M.; Stell, G.; Hall, C. K. J. Chem. Phys. 1976, 65 (6). 2161.
The Journal of Physical Chemistry, Vol. 97, No. 23, 1993 6303 (18) Kincaid, J. M.; Stell, G.; Goldmark, E. J . Chem. Phys. 1976,65 (a), 2172. (19) Kincaid, J. M.; Stell, G. J . Chem. Phys. 1977, 67 (2), 420. (20) Kincaid, J. M.; Stell, G. Phys. Lett. A 1978, 65, 131. (21) Young, D. A.; Alder, B. J. J. Chem. Phys. 1979, 70 (l), 473. (22) Landau, D. P. Phys. Reu. B 1979, 21 (3). 1285. (23) Binder, K.; Landau, D. P. Phys. Reo. B 1980, 21 ( 5 ) . 1941. (24) Selke, W.; Binder, K.; Kinzel, W. Surf. Sci. 1983, 125, 74. (25) Schweika, W.; Binder, K.;Landau, D. P. Phys. Rev. Lcrt. 1990,65 (26), 3321. (26) Bell, G. M. J. Phys. C: Solid Srare Phys. 1972, 5, 889. (27) Bell, G. M.; Salt, D. W. J. Chem. Soc., Faraday Trans. 2 1976,72, 76. (28) Fleming, P. D., 111; Gibbs, J. H. J . Srat. Phys. 1974, I O (2), 157. (29) Fleming, P. D., 111; Gibbs, J. H. J . Srat. Phys. 1974, 10 ( S ) , 351. (30) Sasai, M. J. Chem. Phys. 1990, 93 (IO), 7329. (31) Sastry, S.;Sciortino, F.;Stanley, H. E. Submitted for publication in J. Chem. Phys. (32) Meijer, P. H. E.; Kikuchi, R.; Papon, R. Physica 1981,109A, 365. (33) Meijer, P. H. E.; Kikuchi, R.; Van Royen, E. Physica 1982, IISA, 124. (34) Van Royen, E.; Meijer, P. H. E. Physica 1984, 127A, 87. (35) Van Royen, E.; Meijer, P. H. E. Physica 1986, 139A, 412. (36) Lawrence, P. J.; Rossiter, P. L. J . Phys. F: Met. Phys. 1986,16,543. (37) Kikuchi, R. Phys. Rev. 1951, 81 (a), 988. (38) Kurata, M.; Kikuchi, R.; Watari, T. 1.Chem. Phys. 1953,21 (3), 434. (39) Hijmans, J.; De Boer, J. Physica 1955, 21, 471. (40) In some lattices it is possible for overlapping subclusters to generate
new subclusters. For example, in a triangular lattice with a hexagon (six triangles with a " m o n vertex) basic cluster, two hexagons may overlap to generate thediamondsubcluster (two triangles sharing anedge). Twodiamond subclusters can overlap to generate the triangle subcluster, which cannot be generated from the hexagon basic cluster. (41) Barker, J. A. Proc. R. Soc. London 1953, A216, 45. (42) Kikuchi, R. J. Chem. Phys. 1974,60 (3), 1071. (43) Kikuchi, R. Acra Merall. 1977, 25, 195. (44) Ackermann, H.; Inden, G.; Kikuchi, R. Acra Merall. 1989,37 (I), 1.
