- Email: [email protected]
Calculations of correlation functions for two-dimensional square ice
Chemical Physics 43 (1979) 329-339 0 North-HollandPublishingCompany
CALCULATIONS OF CORRELATION FUNCTIONS FOR TWO-DIMENSIONAL SQUARE ICE * A. YANAGAWA and J.F. NAGLE Deprtnrem
of Physics, Carnegie-Mellon University, Pittirtsburgh, Pennsylvania 15213, UsA
Received2 April 1979 Two approximation methods, the Monte Carlo method and the series expansion method, are used to calculate correlation functions, including the Kirkwood correlation functiong, in twodimensional square ice. The approximation methods are tested on the exact values of the polarization correlation function G and the self-correlation functiong(‘).Both approximate methods are reasonably precise and both methods give values of G andg(‘) in good agreement with the exact values. Using both methods we estimate that g = 1.5 e 0.3. Just as for three-dimensional ice this value ofg is significantly smaller than G, thus supporting the assumption of dimensional similarity in the behavior of the correlation functions. Our results strongly indicate thatg is a discontinuous function in the ideal ice limit as the Bjerrum fault concentration is varied. From these results, previous results for threedimensional ice and general considerations relatingg and G it is concluded that for real ice in three dimensions with a non-zero concentration of Bjerrum faults the values ofg and G are both close to three.
1. Introduction In the previous paper [1] the fundamental theory of the dielectric constant of ice was discussed. In contrast this paper is exclusively concerned with the computation of statistical quantities such as the Kirkwood correlation factor g which appears in the KirkwoodFrohlich theory [2,3] and the polarization factor G which appears in the Onsager-Slater theory [4-61. For isotropic systems g=
lim &J-+m
lim /J.-2 C (B1
‘Pj),
(1)
v+ca i
where the limits mean that the pair correlations of a central dipole labelled 1 are summed over a sphere w which becomes large, but only after the volume Y of the sample becomes infinitely large. Thus, there arealways an infinity of correlations in the system which are not summed in (1). In contrast, all correlations in the system are summed in (2) which gives G for an isotropic system, G=
lim p-2 ~(~l+jj. V+=
(2)
’ This work is based upon the Ph.D. thesis of A. Yanagawa, Carnegie-Mellon University, USA (1979).
As discussed in the preceding paper [l] the “dipole” moments pi may be merely indicators of the orientation of individual units and do not imply dipole-dipoIe interactions. In fact, in neither the KirkwoodFrohlich theory nor the Onsager-Slater theory are dipolar interactions considered in performing the statistical averages in (1) and (2). Such dipolar interactions are presumed to be accounted for in a different part of the Kirkwood-Frohlich theory [2,3,7] and are presumed to be absent in the Onsager-Slater theory ii 1. Therefore, only the short range ice rule interactions are included in the g and G calculations. Until the work of Rahman and Stillinger [8] and Stillinger and Cotter [9] it was believed that g = G because the decay of correlations at large distances would make the difference between g and G negligible, since the difference consists only of the sum of those correlations which are very far away. The conventionaj wisdom that supported this conclusion is that for most models, such as lattice gas or spin models which are not at their critical point, the correlations decay exponentially as (PI .P,) 0: emKr,
K > 0.
(3)
However, for the family of models which strictly obey the Bernal-Fowler-Pauling [lO,l l] ice rules the correlation functions decay more slowly [12] asf(8) rsd,
A. Yamgawa. J.F. iVagIe/Corre[ation functions for ice
330
where d is the dimensionality of the system. This family of models includes the F-model [I31 and the Slater KDP model [6] as well as the model of ice which we will henceforth call the “ideal ice model”. This slow decay of correlations in ideal ice permits the possibility, as Stillinger and Cotter [9] pointed out, that g may be smaller than G. (See the appendix for more discussion of this point and its relation to shape dependence of G.) This slower power law decay is related to the conservation law that the polarization in all layers of the crystal is the same and from this a convincing demonstration of the power law decay is provided by StilIinger and Cotter [9]. But in real ice with a finite concentration of Bjerrum defects the situation is quite different than in ideal ice as regards the decay of the correlation functions because the Bjerrum defects break the polarization conservation laws. Stillinger and Cotter’s demonstration of power law decay does not continue to hold. 1ndeed;if one follows that demonstration, then an atmosphere of one net L Bjerrum defect in front of a central dipole and one net D defect behind it nullifies the long-range power law behavior in the correlations, and this is precisely the atmosphere of defects requlred to maximize the entropy while providing overall charge neutrality. Thus, the relation (3) holds for real ice with a finite concentration of Bjerrum defectswith a very high degree of certainty, and further support for this comes from the series expansions of Stillinger and Cotter and of this paper. If the energy of a Bjerrum fault is 2w, then A = tanh@ introduced by Stillinger and Cotter [9] is a useful parameter to characterize ice models. The value X = 1 requires no Bjerrum faults, which is the ideal ice model. Real ice has X smaller than but very close to one. In terms of h. the conclusions of the preceding paragraph can be expressed as follows. Eq..(3), when inserted into (1) and (2), leads immediately (see appendix) to the conclusion that G(X) =g(A)
for h
(4)
which includes real ice. However, when X = 1 for ideal ice, the slower power law decay of the correaltion functions makes it possible that G(l) #g(l).
(5)
The result g(1) = 2.1 from Monte Carlo calculations [S] and G(1) = 3 from series expansions [14] suggests
that this possibility indeed occurs for threedimensional ice models. The conclusion drawn from eqs. (4) and (5) is that at least one of the correlation functions G or g is badly behaved, that is, discontinuous in the X = 1 Limit of ideal ice [14]. This discontinuity has an important implication for computational accuracy for real ice with small values of 1 -h because approximate methods, such as series,expansions or Monte Carlo calculations, cannot reproduce such a discontinuity and will give inaccurate results near it. However, if one of the functionsg or G is better behaved near X = 1, then the chance of computational success using approximate methods is enhanced for that function. Since by eq. (4),g = G for real ice, the choice of functions to compute is largely dependent upon this consideration of which function is better behaved at X = 1. It has been argued by Nagle [7,14] that the function C is better behaved at h = 1. The primary evidence came from the two-dimensional analogue of ice for which the series expansions for G adequately reproduce the exact result, G(1) =9/x=2.86 [14,15]. Assuming that twodimensional ice behaves qualitatively the same as threedimensional ice then suggests that g(X) is the badly behaved function at h = 1. However, it is clear that the preceding assumption of dimensional similarity of the correlation function behavior which leads to this conclusion should be tested as far as possible. For this purpose it is desirable to apply to two-dimensional ice all the approximate computational methods applied to three-dimensional ice. This computational program is also desirable because the approximate methods can be tested against exact results available for two-dimensional ice but not for threedimensional ice. The most conspicuous computational gap is that for ?wo-dimensional ice no Monte Carlo calculations have been performed, and there are no reliable calculations ofg. In this paper we fill this computational gap by performing Monte Carlo calculations ofg for twodimensional square ice and we also perform series expansions to check the values of close neighbor correlation functions. Basically, this closely parallels the work of Stillinger and co-workers [8,9] on tbreedimensional ice. However, we go beyond this work by relaxing an arbitrary constraint in previous Monte Carlo calculations [8] (introduced to save computer time) so that we can also calculate G to compare to the ex-
A. Yanagawa, J.F. Nagle/Correlation firnctions for ice
act value. As is discussed more fully in section 5, the agreement for G(1) is good by all methods andg(l), although much less accurate, appears to be substantially smaller than G(1). These conclusions strengthen the assumption of dimensional similarity and the advantageousness of computing G rather thang for real threedimensional ice. There is also an exact result for the non-trivial selfcorrelation functiong(“) of twodimensional ice [16]. Both approximate methods give results in very good agreement with this exact result.
2. Monte Carlo method of calculation The Bemal-Fowler-Pauling [lo,1 l] ideal ice model consists of intact water molecules fully bonded to four neighbors situated on a crystal lattice. Fig. 1 shows a typical proton state for the two-dimensional squareice lattice model that is the subject of this paper. All such proton states are assumed in this paper and in refs. [8,9] to have the same energy. The Monte Carlo method approximates statistical averages over all states in eqs. (1) and (2) by averages over a pseudorandom subset of such states. For the ice
problem the method for generating such a subset was pioneered by Rahman and Stillinger (RS) [8] _Starting from an arbitrary state, e.g. the one shown in fig. 1, a pseudorandom number generator (PRNG) selects a site, e.g. site a, and one of the two protons on site Q, e.g. the one forming the H-bond to site b. This proton
(al
331
(b)
Fig. 2. (a) A simple “tadpole” chain. Removal of the tailab produces a simple “cycle” bcb. (b) The chain brb’ is called a “loop”. b’ is the image of b in the next unit cell required by the periodic boundary conditions.
is then moved from a to b, forming a negative ion on a and a positive ion on 6. Next, one of the two protons previously on site b is chosen by the PRNG, e.g. the proton forming the bc bond. Moving this proton to site c advances the positive ion from b to c. Repetition of this process constructs a chain of displaced protons, with a negative ion at the beginning of the chain and a positive ion at the end. Eventually, the chain intersects itself, sometimes at the beginning site, but generally in a middle site to form a “tadpole” shown in fig. 2a,wirh a negative ion at a and a positive ion at b. When thisoccurs, the protons on the tail of the tadpole, the portion ab, are returned to their original positions, so the final change consists of a chain of displaced protons around the cycle bcb. An example of a permissible cycle in fig. 1 is abcdefa. Notice that abcda is not a permissible cycle. Incidentally, it has been proven by Griffiths that this state generation method is ergodic [17]. It is obviously desirable to employ periodic boundary conditions, which raises an important point. Sometimes the chain intersects itself in one of the image lattices as shown in fig. 2b. Such a chain, bcb’, is called a loop, in contrast to the cycle of fig. 2a. Permissible loops in fig. 1 are O&I and adehfa. The difference between loops and cycles is that cycles do not change the net polarization P of the finite lattice with N = IZ* sites, where P is defined as N
P=
cpj, i=l
whereas loops change the components
Fig. 1. A typical proton configuration in twodimensional square ice. The 3 X 3 lattice has periodic boundary conditions.
of the polaiiza-
tion P, and/or Py by multiples of 21j2 np. In fig. 1 the polarization is P, = 3 (p/21j2) and P,, = - 3 (p/21j2). Each of the loops odgu and ndehfa change PY to +3 (c(/2l/*)i (There are no Py = 0 states for an n X n lattice with n odd.)
332
A. Yanagawa.J.F. Nagle/Correlotio!l functionsfor ice
RS [8] deliberately eliminated loops in their program. Strictly speaking, that was not legitimate since there is no such restriction in eqs. (1) and (2). The reason for this restriction was to save computer time by reducing the number of states to be sampled. This is an unfortunate restriction because it requires G in eq. (2) to be identically zero. Also, in order to obtaing, RS used some plausible, but non-rigorous corrections. Fortunately, it is not necessary for us to make such a restriction of e!iminating loops for twodimensional ice. Thus, we calculate G as well as g. Furthermore, in our program we simultaneously calculate two sets of averages. The first set is for all states including those for which P # 0. Such averages will be designated with subscripts P, e.g. gp.The second set is identical to the averages of RS and includes only those states generated with P = 0. These averages will be designated with subscript 0, e.g.go. With both sets of averages we can test the restriction and the corrections used by RS [8]. Ideally, the lattice should be very large so as to minimize finite size errors. However, increasing size increases the number of states to be sampled and the computer time. This is especially serious because the percentage of loops decreases rapidly with larger lattice size and the system becomes frozen into a particular polarization for long computer times (many state changes). However, for valid comparison to the work of RS it is necessary to have a lattice of “comparable” size. Fortunately, from the point of view of finite size effects a lattice of comparable size in two dimensions is not one with the same number of lattice sites as the one in three-dimensions. Rather, it is one with the same linear dimensions, so an PIX n two-dimensional lattice in two dimensions is “comparable” to an n X n X n lattice in threedimensions. For such comparable lattices the shortest loop is the same length. RS [8] studied a cubic lattice with 4096 sites which has a shortest loop of 16 steps, so we began our work with a 16 X 16 square lattice. For our method which allowed loops this turned out to be the largest lattice for which we could obtain reliable results in a reasonable amount of computer time. For larger lattices the low frequency of loops is a severe problem which induces very long period fluctuations in the averages. However, for 12 X 12 and 8 X 8 lattices the statistical convergence is quite rapid. As was mentioned in the introduction eqs. (1) and (2) are valid for isotropic systems. However, for sys-
terns of lower symmetry, such as hexagonal ice in three dimensions or two-dimensional ice, the correlation functions xi (pl .P~)/P~ must be replaced by the more basic correlations [2,7]
where I’, is the unit vector along the electric field. By symmetry these correlations reduce to those in eqs. (1) and (2) for cubic syst,ems. Twodimensional square ice is an extreme example of asymmetry because two of the six water molecule configurations have no dipole moment. However, this does not qualitatively change the statistical problem so long as the proper orrelations are calculated. One feature induced by this asymmetry which is unlike cubic ice is that the six different water molecule configurations do not occur with equal probability. An exact calculation of Wu [16] shows that each of the water molecule configurations with zero dipole moment has occurrence probability of 19% compared to only 15.5% for the other four configurations. Consider the term in (7) when j = 1, which we will call the self-correlation function g(O). This is the first term in the sum for g. The exact result of Wu requiresg(O) = 0.930 for square ice, compared to g(O)= 1 for cubic ice. In our computer program the correlation functions in (7) are computed for E = E, along the x axis and also for E = I?,, along they axis for all states generated by the Monte Carlo method. By symmetry these correlations are equal for the sum of all states, but not for a subset of states. Thus, one monitor of the reliability of the sampling procedure is the difference between the correlations with the two different field directions.
3. Results of Monte Carlo calculations 3. I. Distribution of chain lengths Fig. 3 shows the probability p(l) of occurrence of cycles and loops of length 1 for a 16 X 16 lattice.Table 1 gives pertinent statistics on the average cycle length, the average loop length and the frequency of.occurrence of loops_ Most of the cycles are just the simplest square cycle of length 4. The fact that the average cycle length is less than 6 means that only a small fraction of the system is affected by a single configura-
A. Yanaga wa, J. F. Nagle/Conelation
0.8-
functions for ice
333
- 0.4
- 0.3
o8 >
-0.2
-Of
I
Fig. 3. Probability of chain Iengthsp(I) versus Iength 1. The solid curves show p (I) (open circles) and p(l) X 10 (solid circles) for cycles with the scale on the left. The dashed curve showsp(0 for loops with the scale on the right. The lattice sizeis16X 16.
tional change and so most of the correlations do not change. Since the computation of the correlations consumes a large fraction of the computer time, it was decided to compute correlations only after k = 25 successive configurational changes. In comparison, RS [S] used k = 200, but this was necessary because, even though the average cycle length is larger for three dimensional ice, there are many more molecules for lattices of comparable size. 3.2. Computation of G by Monte Carlo method G is technically an easy and inexpensive quantity to compute. Consider
(8)
Thus, for E in they direction one only needs to add they-components of the dipole moments in one row of the lattice, multiply by the number of rows, and square the result to pbtain P# for each state generated. Fig. 4 shows computed vaBues of G as a function of the number of configurations sampled for lattice sizes 8 X 8,12 X 12 and 16 X 16. The average over sums of 10000 configurations show typical fluctuations which become larger with larger n. This trend with larger n is to be expected from the occurrence of loops shown in the last column of table 1 because, as was discussed in section 2, the only way to change P ( and therefore from eq. (9) the value of G) is via loop changes. Fig.4 shows that for a 16 X 16 lattice one clearly does not sample as many different P’s in 10000 configurations as for an 8 X 8 lattice. Of course, the cumulative averages over all generated states, also shown in fig. 4, show sinaller fluctuations, akhough long fluctuation waves are clearly apparent for the 16 X 16 lattice. Averaging over the 15,20 and 24 runs, each run involving 10000 configuration changes, for the 8 X 8,12 X 12 and 16 X 16 lattices, respectively, gives our final estimates of G with ENS errors which are presented in the third column of table 2.
3.5 30 PC. Ee
Since the central molecule labelled 1 can be any ofthe N molecules, it is convenient to sum over all of them which gives
(8 x 8) ___a.-.
--
_L
.e-*-
3.0 25 .
3.0
11 116x161
20
I I
5
10
c/10000
% of loops
8X8 12x 12 16x 16
Average loop length
Average cycle length
14.1
15.33
5.42
3.0 0.5
23.83 32.62
5.53 5.53
a-4.
&---
--
.
Lattice size
-
3.5
4.0
Table 1 Distribution of cycles and loops
*
15
20
c.4
,
25
Fig. 4. Monte Carlo resulti for G for three sizes of lattice as a function of number of configurations C sampled. The solid line with sotid circles shows the value of G averaged eve: the preceding 10000 configurations sampled. The dashed linewith crasses shows the value of G averaged over all preceding confgurations sampled.
A. Yanagawa, J.F. Nagie/Correlation functiom for ice
334
Table 2 Monte Carlo results for the polarization function G, the self-correlation functiong (‘1 and the Kirkwood correlation functiong.The subscripts P and 0 refer to all states and states with zero polarization, respectively, with relative numbers n,/n~~ of the two kinds of states btticc
“&P
G
gp
gp
gP
0.33 0.35 0.33
2.88 f 0.04 2.78 c 0.06 2.86 r 0.11
0.942 f 0.002 0.937 I?0.001 0.933? 0.001
0.920 f 0.003 0.927 * 0.002 0.925 f 0.002
1.41 1.40 1.43
size
8X8
12x 12 16 x 16
3.3. Result for g(o) and error analysis
RO
1.52 1.53 1.41
given for all the entries in table 2 are t on2 as obtained by the second method.
The results for g(o) for two different samplingprocedures,&
for only P=G states and $1
for all states,
are given in table2. Although these estimates are in reasonably good agreement with the exact value, g(O) ~0.930, the statistics for a given lattice give standard
deviation errors that are generail smaller than the real error. However, it is ciear that gb) IS steadily decreasing toward the exact value and g&O)is steadiIy increasing as the size of the lattice increases. Doing a simple extrapolation versus l/r1 for these n X n lattices gives us a final Monte Carlo e;tim&e,g(@ = 0.930 -t-0.005. For completeness, let us men&on how the standard deviation errors were corn “ted. In our first method the individual values ofgi( B) for each state obtained (after 25 configuration changes as explained in section (2.1)) were used as the n data points in the usual formu-
la for normal distributions,
02 = m
2 [g,!0’-~q2/,(,-l).
i=l
However, even after 25 configuration
changes it might
be argued that the new state is not randomly independent of the previous one. Therefore, in our second method of obtaining statistics g(O)was first averaged over sequences of 10000 configuraticn changes(m = 400 states = 10000/25) and each of these averages $‘) should be much more independent than the individual gj”) used in the fikmethod. Then these&‘) replace the gi”) in (10) and n becomes of the order of 10 to 20 depending upon the time the computer was run which depended upon the lattice. Just as in the
first method jj”) is the grand mean. This second mkthod gave somewhat larger values of urnIbut by less than a factor of 2, than did the first method. The errors
3.4. Results for the correlation fitnctions and the Kirkwood correlation factor Following RS [8] it is useful to consider the sum of correlation functions over the tih shell of neighbors, defined as.g@. For an infinite system m
g=
c
p). v=o
(11)
The shell-v = 0 just gives g(O)which is the self-correlation and the shell Y = 1 gives g(l) which is the sum of the four nearest neighbor correlations, etc. The maximum coordination shell that we routinely consider is v = 16 for the 16 X 16 lattice. The number of neighborsin shells 0 to 16 is 10 1. In fig. Sa is shown gcV)versus v for the two different sampling procedures, gk) for P = 0 states and gg) for unrestricted P states. Both calculations give values of g(O) which are very close to the exact value 0.930 (see table 2) and both give a negative gc2). Further correlations are rather undistinguished, except that g;) tends to be greater than gt). This becomes apparent when the sum of alI correlations up to the vth shell, defined as ”
g””=
izdi),
02)
is plotted versus v as seen in fig. Sb. Clearly gp” continues to increase as a function of v and gr continues to decrease. Although this may seem alarming, it is expected. As v becomes large the vth coordination shell contains most of the molecules in the finite lattice. Therefore, one expects g? to become approximately
A. Yanagawa, LF. NaglefCorrelation
335
functions for ice
g can be made. The combined estimates are summarized in table 2 for the two casesgp andgo. No errors are given because of the intrinsic uncertainty in the procedure. The number of state changes made was 100000 for the 8 X 8 and for the 12 X 12 and 160000 for the 16 X 16 lattices with correlations calculated after every k = 25 states. The fact that estimates ofg made from bothgg and from our new gg are in substantial agreement as seen in fig. 5c lends support to the correction scheme of RS [8].
4. Series expansion calculations of the correlation functions The general method of weak graph series expanFig. 5. Monte Carlo results for correlation functions for a 16 X 16 lattice from a sample of 10000 configurationsversus shell number V. (a) The shell correlation functionsg$?) (solid line and circles) andgg) (dashed line and crosses). (b) The summed correlation functionsgy (solid line and circles) and gr (dashed line and crosses). (c) The corrected correlation functionsggi (solid line and circles) andggy (dashed line and crosses).
sions [M] has been applied to the correlation functions in cubic ice by Gobush and Hoeve [ 191 and Stillinger and Cotter [9] _In this section analogous computations are reported for square ice. The StillingerCotter [9] expansion variable, X = tanhlvfi, where 2w is the energy of a Bjerrum fault, will be used. The basic correlation functions are g’O”(h) = (3//J2) (PO- $T) (I$ - &)A
equal to G and gr to become approximately equal to 0 and we verified this for a few of our calculations. Since there is no clear value of n at which to cut off the summation, this poses the interesting problem of extracting from fig. Sb a reasonable estimate ofg. RS [8] proposed !m obviously non-rigorous correction to gr as follows. The zero P constraint was supposed to introduce a linear negative bias proportional to -l/N for each molecule within the uth shell so that
where p. and 4 are the dipole moments of the 0th and rth molecules. These correlation functions may be obtained most easily from the modified partition function,
Z($-J,~)=
C cnrv@
states
(13)
where NO, is the number of molecules in the system which are within the tih shell. It is hoped that gz will then give a better approximation tog for values of u such that 1 - (NW/N) is of the order of l/2. We now propose also to correct gp” in the only way consistent with the aforementioned linear RS correction, namely gp” =g;; (1 -A&/N)
+ G(N,,/N),
(16)
zo2,,
where the states have n Bjerrum faults and where ZO~~=~~P~~~@O.~E~IPOI)+~(J~,
gr =gg (1 -N,,JN),
(15)
~&/ly-I)]
(17)
by the formula a2 m Z(~o,o;) g(O” (1) = 3 -Qd-&oe=o=‘yr
(18)
Using the weak graph method eq. (16) becomes
(14)
where it is now supposed that each molecule within the uth shell adds a positive linear bias (G-gg)/N. Fig. SC shows gz and gE$ versus v. Both g$ and g$ now have reasonably flat portions from which estimates for
as is thoroughly explained in the literature [14,18]. The expression (19) is then expanded term by term where each term corresponds to a subgraph of the square lattice. After performing the operation in (18)
A. Yanagatvu. J.F. NaglejCorrelation functions for ice
336
i = vi
=
la
lb
-JF
fi
2a -1
Zb -I
3a
3b
~3E-3/z
4 9
Fig. 6. Vertex types i and vertex weights “i in the series expansions. one has a sum over the subgraphs of the lattice
g@‘)(X)= 5 w(~),
(20)
where w(((,) is the contribution or weight of subgraph ai. It is given by a simple product W(@)‘(-Z)’
n Vi(U),
vertices in U
(21)
where t =X/3, e is the number of edges in the graph aand ui((g) is the weight of the ith vertex of &These vertex weights vi(@) depend only upon the local coordination of the ith vertex in the subgraph @. Whenthe electric field is in the vertical direction, there are only seven different types of local coordination that have non-zero vertex weights and these are shown in fig. 6 with their weights. Types la, lb, 3a and 3b may only occur on vertices 0 and r.
The correlation function g(O)can be treated in a similar way except that r and 0 are the same vertex. In this case the vertex weights for the special vertex 0 are ~2, = 1, u’& = -2 and u4 =O. There is also a contribution w(0) = 1 from the null subgraph. Having found the graph weights w(($) in (20) the calculation reduces to counting the subgraphs Qj of the different types and details will not be given. Results are expressed for gc”) which is the sum over alI g@) where r is in the tih shell around 0. We have obtained g(O) = 1 -4z4 -68~~ - 16~1~ - 1384z12, g(11=4z-8z3-32z5-8Clz7-648zg,
x
g(*) = -16~~ -32~~ - 176~8 -52&O, gt3) = 4z* - 8z4 - 64z6 t 8z8 - 744z’O, g(4) = 8,s -40~5 - 72~’ - 40829, gts) = 8z4 -3226 - 80~8 - 145zfi,
(22)
where z =X/3, to the highest order shown. Unlike the threedimensional cubic ice case [9] it seems that one of the correlations, gt3) has a higher order term, +8z8,. which is not negative, though its magnitude is very small. In general, hcwever, the higher order terms are negative as in the case for cubic ice. An independent check on these series is obtained by adding them algebraically to obtain the series G(A). This latter series has been obtained using a renormalized weak graph method [14] which requires less graph counting and is therefore less subject to errors. From
1141 C(X) = 1 t 4z t 4z2 t 4z3 - 12z5 -20~~ - 12z7 t 16z8 -(32/3)zg
+ (704/3)z10 f 9762” •r-. . . .
(23)
In addition to theg@‘) shown in (22), all other terms up to and including z7 were found in the remaining correlation functions and the sums for each order agree with (23) up to and including 7th order. Numerical estimates of the correlation functions may be obtained in a number of ways. A particularly illuminating procedure is shown in fig. 7 for ideal ice (X= 1). Defimeg$) to be the partial sum of the terms up to and including the zn term in the series for gt”). In fig.? &‘I for zj=O and 1 are plotted versus l/n. For v = 0 fig. 7 shows that gA”) appears to be converging rather well in the limit l/,t + 0 to the exact value g(O)= 0.930 although our best extrapolation slightly overestimates the exact value by 0.007. The Monte Carlo estimate &) is also shown in fig. 7. For v = 1 ii 7 shows that g I) has much more variation than did ’ value as l/n + 0 appears in this gn8) , but the limiting case to be consistent with the Monte Carlo estimate. For v = 2 to 5 similar plots (not shown) suggest limiting values of g(y) which are very close to, although slightly larger than, the Monte Carlo results. Other ways to estimateg(V) involve examining the ratios of successive terms or the use of Pade approximants PO]. The ratios of successive terms alternate between large and small, so one has to attempt to extrapolate two very short subsequences. Also, the series are a bit short for the full power of the Pade approximant to be high!y beneficial. Therefore, these techniques do not improve much upon the simple graphical technique illustrated in fig. 1. Our numerical estimates for the individual correlation functions g@)are collected in table 3. The values
A. Yanagarua. J.F. NaglefCorrelation
functions
337
for ice
3.0
2.0
1 .o
0
0.2
0.4
0.6
0.8
1.0
Fig. 8. G(k) andggE(X) versus h from series expansions. The number v of correlations summed is indicated for eachgz&Q
MC cases. The last figure in each column g$ is an approximation for the Kirkwood correlation function
Fig. 7. Analysis of the convergence of correlation function series using a l/n pl0t.g:) IS . the partial sum of the first n terms in the series for ideal ice (h = 1). The arrow shows the exact value ofg(O) and the error bars show the range of values given by the hlonte Carlo estimates.
& are extrapolated from the series expansions. The vaiues g(v) ,lma?: are obtained by adding the available terms in the series; this gives an upper bound because the higher terms in the series are redominantly negative. For comparison the valuesg&P obtained from the Monte Carlo method are also given. In table 3 the sum of g'$from n = 0 to II = u are also given in the columns labelled g,“, where x refers to.the nmax, ex or
Table 3 Results of series and Monte Carlo methods. LJis the shell number,g@) give the sum of correlation functions in each shell and guV give the sum of correlation functions over shells up to and including the vth shell. The subscripts nmax, ex and hfC refer to truncated series approximations, sxtrapolated series approximations and Monte Carlo approximations, respectively. All results are for the limit X = 1
g(l). It is most instructive to consider the series approximation as a function of X. Fig. 8 shows G(A) and gg(X) for v = 1 to 5 versus X. The curve for G(X) proceeds monotonically as a function of X to the e.xact value G = 9/n = 2.865 at h = 1. Earlier estimates using three different extrapolation methods [14] gave G(1) = 2.848,2.855 and 2.898 which are closely clustered around the exact value. In contrast, the estimated errors for g::(X) shown in fig. 8 only for&(X) are considerably greater than for G(X). Most important is the behavior of g’$(‘R) as a function of h. As can be seen in fig. 8 each gg(h) has a maximum for h < 1 which becomes larger with increasing Y and whose posiiion shifts to larger h. We interpret this as a precursor of a discontinuous decrease ing(A) at X= 1. For example, since g::(X) is a truncated approximation to g(X) and since the individual correlation functions g::(X) are likely to be continuous functions of X,gi:(X) should be continuous, A continuous approximation to a discontinuous g(X) would have precisely the behavior shown by g’$(h) in fig. 8 because in order to simulate the discontinuous decrease at h= 1, the continuous approximation would reach a maximum for X < 1 and then decrease.
av
gex 0
0.937
0.937
0.930
0.937
0.937
1 0.836 2 -0.277 3 0.247 4 0.078 5 0.018
0.681 -0.287 0.195 0.058 -0.010
0.768 -0.318 0.128 -0.064 -0.026
1.77 1.50 1.74 1.82 1.84
1.62 1.33 1.53 i.58
1.57
0.930 1.68 1.36 1.49 1.43 1.40
5. Conclusions First, we have tested two approximate methods, the Monte Carlo method and the series method, against some exact resuits for correlation functions in twodimensional square ice. The series method gives very
A. Yanagawa, J.F. Nag/e/Correlation
338
good estimates [14] of the polarization factor G which equals 9/n exactly and of the self-correlation function g(O)(see table 3) which equals O.Y30 exactly. The Monte Carlo method is less accurate for G (see table 2) than the series method, but the standard deviation errors include the exact result. For g(n) the Monte Carlo method is even more accurate than the series method, although not much effort was expended to extend the length of the series.
Second, we have used both the Monte Carlo and the series methods to calculate the Kirkwood correlation factor g for which there is no exact result. These estimates of g are considerably poorer than our estimates of G and g(o). Part of the reason for the relative inaccuracy in the estimates ofg is intrinsic to the calculation methods. In the Monte Carlo method the finite Iattice size imposes a finite, but arbitrary, cut-off in summing the correlation functions. It also requires an ad hoc correction such as the one introduced by RS [8], but our inclusibn of states with non-zero polarization in the calculation lends support to this correction. In the series method the calculation of g requires not only series extrapolations, similar to those in the G or g(O)calculations, but also an additional summation over correlation functicns. Although these considerations make it difficult to estimate g and especially difficult to assign errors, both methods clearly indicate that g(l) is much smaller than G(1) with even more certainty than for threedimensional ice. From tables 2 and 3 the Monte Carlo method appears to give g(1) about 1.4 and the series method appears to give g(l) about 1.6, we estimate that g(l)=
1.5kO.3
(24)
for square ice. We can now conclude that the correlation func. tions in twodimensional square ice behave similarly to those in the threedimensional ices, supporting the assumption of dimensional similarity discussed in the introduction. In particular, the previous conclusion that eqs. (4) and (5) held in three-dimensional ice 17, 141 now hold at least as firmly in two dimensions. Therefore, either g(A) and/or G(X) is discontirmous at X = 1. The fact that both the Monte Carlo method and the series method accurately reproduce the exact value of G(1) and that the series give a smooth monotonic curve for G(1) suggests that G(X) is not strongly discontinuous at h = 1. In contrast, the relatively in-
fitnctions for ice
accurate estimates of g(1) may not be due entirely to intrinsic computational problems. Instead, they are likely to be due in part to g(h) being discontinuous and the difficulties in approximating such a function with continuous approximations. This conclusion is strongly indicated by the series expansion method which would be expected to approximate a discontinuousg(h) by the sequence of approximationsgi?(X) shown in fig. 8. Now that the approximate methods for correlation functions have been tested against exact results in twodimensional square ice and now that the resuits of the calculation have been shown to be similar to the re-
sults for three-dimensional ice and to indicate that has a discontinuous decrease at h = 1, the previous conclusion [14], that for real three-dimensional ice with a small concentration of Bjerrum defects
g(h)
g(X) = G(A) = 3,
(25)
can be drawn.even more strongly.
Acknowledgement
We wish to thank Professors J.S. Langer and R.B. Griffiths for considerable discussionand R.B. Griffiths for his ergodicity proof which was instrumental in excisinga bias from an earlier version of the Monte Carlo program. This research was supported by the United States National Science Foundation grant DMR 78-10226.
Appendix: Long range correlations and shape dependence
If the correlation functions behave as (JII .I+> a f(0) r,F*-“, then for a system with periodic boundary conditions, G-g=
lim
‘K *-
lim
c
Vdrn
rj>rK
f(e) r,+n
2K rw z
15, ‘K
>i_
$
s
f(0)r-2-~rdrdB,
(26)
’ ‘K
is clearly zero for 7 > 0 or for correlations given by
eq. (3); For ideal ice n = 0, so in order for the integral
A. Yarugawa, J.F. Nagle/Correlation fiolctions for ice in (26) to remain finite, there must be strong cancellations in the 0 integral, which in fact occur in Sutherland’s formula [12]. However, this raises the .jnteresting question whether G depends upon the shape r,(e) of the sample. There are four points to be made: (1) The exact calculation of G involves the use of the transfer matrix for rectangular n X m shapes [16]. Taking the largest eigenvalue corresponds to letting m + m before n -+ ~0, so this method might be expected to give G for long thin rectangles, although the final formula is isotropic. The series method never refers to shape of the sample. The fact that the Monte Carlo result for G for square n X n shapes agrees so
well with the other two methods is inconsistent with a shape dependence unless it can be argued that the exact method and the series method implicitly, and mysteriously, require an n X n lattice. (2) We performed some Monte Carlo calculations for n Xm lattices, jz >m. For the field along the short m side, for which one obtains the best statistics because loop changes are more frequent, the value of G was 9/n within statistical error. However, for the field along the long IZside, for which the statistics were’ poorer,. the value of G decreased as n/m increased, reaching G = 1.85 +0.26 for n/m = 16/6, although we
do not have great confidence in this last result. Experimentally, fields are usually applied along the short dimension for which our calculations indicate G to be independent of n/m. (3) The asymptotic expression for the correlations used in (26) is only expected to hold when thejth molecule is far away from any surface, and surface corrections might conspire to negate any shape dependence. (4) One does not expect any shape dependence in the theoretical dielectric constant in either the Kirkwood-Frohlich theory or in the Onsager-Slater theory. In particular, the usual macroscopic depolarization factors have been taken into account by using
339
the internalE field in these theories. After consideration of these points we are inclined to discount shape dependence in G.
References [l] J.F. Nagle,Chem. Phys. 43 (1979) 317, preceding paper. [2] J.G. Kirkwood, J. Chem. Phys. 7 (1939) 9 11. [3] H. Frohlich, Theory of dielectrics. Dielectric constant and dielectric loss (Clarcndon Press, Oxford, 1949). [4] L. Onsager and M. Cupuis, in: Electrolytes, The electri cal properties of ice, ed. B. Pesce (Pergamon Press, New York, 1962) pp. 27-46. [5] L. Onsager and L.K. Runnels, J. Chem. Phys. 50 (1969) 1089. 16) J.C. Slater, J. Chem. Phys. 9 (1941) 16. [7] J.F. Nagle, I. Glaciology 21 (1978) 73. [8] A. Rahman and F.H. Stillinger, J. Chem. Phys. 57 (1972) 4009. [9] F.H. Stillinger and M.A. Cotter, J. Chem. Phys. 58 (1973) 2532. [IO] L. Pauling, J. Am. Chem. Sot. 57 (1935) 2680. [ll] J.D. Bernal and R.H. Fowler, J. Chem. Phys. 1 (1933) 515. [12] B. Sutherland, Phys. Letters 26A (1968) 532. [I31 F. Rys, Helv. Phys. Acta 36 (1963) 537. [14] J.F. Nagle, J. Chem. Phys. 61 (1974) 883. [15] B. Sutherland, C.N. Yang and C.P. Yang, Phys. Rev. Letters IV (1967) 588. [16] E.H. Lieb and F.Y. Wu, in: Phase transitions and critical phenomena, Vol. 1, eds. C. Domb and M.S. Green (Academic Press, London, 1972) ch. 8. See especiaUy p. 450 forg(O)(l). [I71 A. Yanagawa, Ph.D. Dissertation, Carnegie-Mellon University, USA (1979). [IS] J.F. Nagle, in: Phase transitions and critical phenomena, Vol. 3, eds. C. Domb and M.S. Green (Aczxiemic Press, London, 1974) p. 653. [19] W. Gobush and C.A.J. Hoeve, J. Chem. Phys. 57 (1972) 3416. 1201 D.S. Gaunt and A.J. Guttmann, in: Phase transitions and critical phenomena, Vol. 3, cds. C. Domb and h1.S. Green (Academic Press, London, 1974) ch. 4, p. 181.
CALCULATIONS OF CORRELATION FUNCTIONS FOR TWO-DIMENSIONAL SQUARE ICE * A. YANAGAWA and J.F. NAGLE Deprtnrem
of Physics, Carnegie-Mellon University, Pittirtsburgh, Pennsylvania 15213, UsA
Received2 April 1979 Two approximation methods, the Monte Carlo method and the series expansion method, are used to calculate correlation functions, including the Kirkwood correlation functiong, in twodimensional square ice. The approximation methods are tested on the exact values of the polarization correlation function G and the self-correlation functiong(‘).Both approximate methods are reasonably precise and both methods give values of G andg(‘) in good agreement with the exact values. Using both methods we estimate that g = 1.5 e 0.3. Just as for three-dimensional ice this value ofg is significantly smaller than G, thus supporting the assumption of dimensional similarity in the behavior of the correlation functions. Our results strongly indicate thatg is a discontinuous function in the ideal ice limit as the Bjerrum fault concentration is varied. From these results, previous results for threedimensional ice and general considerations relatingg and G it is concluded that for real ice in three dimensions with a non-zero concentration of Bjerrum faults the values ofg and G are both close to three.
1. Introduction In the previous paper [1] the fundamental theory of the dielectric constant of ice was discussed. In contrast this paper is exclusively concerned with the computation of statistical quantities such as the Kirkwood correlation factor g which appears in the KirkwoodFrohlich theory [2,3] and the polarization factor G which appears in the Onsager-Slater theory [4-61. For isotropic systems g=
lim &J-+m
lim /J.-2 C (B1
‘Pj),
(1)
v+ca i
where the limits mean that the pair correlations of a central dipole labelled 1 are summed over a sphere w which becomes large, but only after the volume Y of the sample becomes infinitely large. Thus, there arealways an infinity of correlations in the system which are not summed in (1). In contrast, all correlations in the system are summed in (2) which gives G for an isotropic system, G=
lim p-2 ~(~l+jj. V+=
(2)
’ This work is based upon the Ph.D. thesis of A. Yanagawa, Carnegie-Mellon University, USA (1979).
As discussed in the preceding paper [l] the “dipole” moments pi may be merely indicators of the orientation of individual units and do not imply dipole-dipoIe interactions. In fact, in neither the KirkwoodFrohlich theory nor the Onsager-Slater theory are dipolar interactions considered in performing the statistical averages in (1) and (2). Such dipolar interactions are presumed to be accounted for in a different part of the Kirkwood-Frohlich theory [2,3,7] and are presumed to be absent in the Onsager-Slater theory ii 1. Therefore, only the short range ice rule interactions are included in the g and G calculations. Until the work of Rahman and Stillinger [8] and Stillinger and Cotter [9] it was believed that g = G because the decay of correlations at large distances would make the difference between g and G negligible, since the difference consists only of the sum of those correlations which are very far away. The conventionaj wisdom that supported this conclusion is that for most models, such as lattice gas or spin models which are not at their critical point, the correlations decay exponentially as (PI .P,) 0: emKr,
K > 0.
(3)
However, for the family of models which strictly obey the Bernal-Fowler-Pauling [lO,l l] ice rules the correlation functions decay more slowly [12] asf(8) rsd,
A. Yamgawa. J.F. iVagIe/Corre[ation functions for ice
330
where d is the dimensionality of the system. This family of models includes the F-model [I31 and the Slater KDP model [6] as well as the model of ice which we will henceforth call the “ideal ice model”. This slow decay of correlations in ideal ice permits the possibility, as Stillinger and Cotter [9] pointed out, that g may be smaller than G. (See the appendix for more discussion of this point and its relation to shape dependence of G.) This slower power law decay is related to the conservation law that the polarization in all layers of the crystal is the same and from this a convincing demonstration of the power law decay is provided by StilIinger and Cotter [9]. But in real ice with a finite concentration of Bjerrum defects the situation is quite different than in ideal ice as regards the decay of the correlation functions because the Bjerrum defects break the polarization conservation laws. Stillinger and Cotter’s demonstration of power law decay does not continue to hold. 1ndeed;if one follows that demonstration, then an atmosphere of one net L Bjerrum defect in front of a central dipole and one net D defect behind it nullifies the long-range power law behavior in the correlations, and this is precisely the atmosphere of defects requlred to maximize the entropy while providing overall charge neutrality. Thus, the relation (3) holds for real ice with a finite concentration of Bjerrum defectswith a very high degree of certainty, and further support for this comes from the series expansions of Stillinger and Cotter and of this paper. If the energy of a Bjerrum fault is 2w, then A = tanh@ introduced by Stillinger and Cotter [9] is a useful parameter to characterize ice models. The value X = 1 requires no Bjerrum faults, which is the ideal ice model. Real ice has X smaller than but very close to one. In terms of h. the conclusions of the preceding paragraph can be expressed as follows. Eq..(3), when inserted into (1) and (2), leads immediately (see appendix) to the conclusion that G(X) =g(A)
for h
(4)
which includes real ice. However, when X = 1 for ideal ice, the slower power law decay of the correaltion functions makes it possible that G(l) #g(l).
(5)
The result g(1) = 2.1 from Monte Carlo calculations [S] and G(1) = 3 from series expansions [14] suggests
that this possibility indeed occurs for threedimensional ice models. The conclusion drawn from eqs. (4) and (5) is that at least one of the correlation functions G or g is badly behaved, that is, discontinuous in the X = 1 Limit of ideal ice [14]. This discontinuity has an important implication for computational accuracy for real ice with small values of 1 -h because approximate methods, such as series,expansions or Monte Carlo calculations, cannot reproduce such a discontinuity and will give inaccurate results near it. However, if one of the functionsg or G is better behaved near X = 1, then the chance of computational success using approximate methods is enhanced for that function. Since by eq. (4),g = G for real ice, the choice of functions to compute is largely dependent upon this consideration of which function is better behaved at X = 1. It has been argued by Nagle [7,14] that the function C is better behaved at h = 1. The primary evidence came from the two-dimensional analogue of ice for which the series expansions for G adequately reproduce the exact result, G(1) =9/x=2.86 [14,15]. Assuming that twodimensional ice behaves qualitatively the same as threedimensional ice then suggests that g(X) is the badly behaved function at h = 1. However, it is clear that the preceding assumption of dimensional similarity of the correlation function behavior which leads to this conclusion should be tested as far as possible. For this purpose it is desirable to apply to two-dimensional ice all the approximate computational methods applied to three-dimensional ice. This computational program is also desirable because the approximate methods can be tested against exact results available for two-dimensional ice but not for threedimensional ice. The most conspicuous computational gap is that for ?wo-dimensional ice no Monte Carlo calculations have been performed, and there are no reliable calculations ofg. In this paper we fill this computational gap by performing Monte Carlo calculations ofg for twodimensional square ice and we also perform series expansions to check the values of close neighbor correlation functions. Basically, this closely parallels the work of Stillinger and co-workers [8,9] on tbreedimensional ice. However, we go beyond this work by relaxing an arbitrary constraint in previous Monte Carlo calculations [8] (introduced to save computer time) so that we can also calculate G to compare to the ex-
A. Yanagawa, J.F. Nagle/Correlation firnctions for ice
act value. As is discussed more fully in section 5, the agreement for G(1) is good by all methods andg(l), although much less accurate, appears to be substantially smaller than G(1). These conclusions strengthen the assumption of dimensional similarity and the advantageousness of computing G rather thang for real threedimensional ice. There is also an exact result for the non-trivial selfcorrelation functiong(“) of twodimensional ice [16]. Both approximate methods give results in very good agreement with this exact result.
2. Monte Carlo method of calculation The Bemal-Fowler-Pauling [lo,1 l] ideal ice model consists of intact water molecules fully bonded to four neighbors situated on a crystal lattice. Fig. 1 shows a typical proton state for the two-dimensional squareice lattice model that is the subject of this paper. All such proton states are assumed in this paper and in refs. [8,9] to have the same energy. The Monte Carlo method approximates statistical averages over all states in eqs. (1) and (2) by averages over a pseudorandom subset of such states. For the ice
problem the method for generating such a subset was pioneered by Rahman and Stillinger (RS) [8] _Starting from an arbitrary state, e.g. the one shown in fig. 1, a pseudorandom number generator (PRNG) selects a site, e.g. site a, and one of the two protons on site Q, e.g. the one forming the H-bond to site b. This proton
(al
331
(b)
Fig. 2. (a) A simple “tadpole” chain. Removal of the tailab produces a simple “cycle” bcb. (b) The chain brb’ is called a “loop”. b’ is the image of b in the next unit cell required by the periodic boundary conditions.
is then moved from a to b, forming a negative ion on a and a positive ion on 6. Next, one of the two protons previously on site b is chosen by the PRNG, e.g. the proton forming the bc bond. Moving this proton to site c advances the positive ion from b to c. Repetition of this process constructs a chain of displaced protons, with a negative ion at the beginning of the chain and a positive ion at the end. Eventually, the chain intersects itself, sometimes at the beginning site, but generally in a middle site to form a “tadpole” shown in fig. 2a,wirh a negative ion at a and a positive ion at b. When thisoccurs, the protons on the tail of the tadpole, the portion ab, are returned to their original positions, so the final change consists of a chain of displaced protons around the cycle bcb. An example of a permissible cycle in fig. 1 is abcdefa. Notice that abcda is not a permissible cycle. Incidentally, it has been proven by Griffiths that this state generation method is ergodic [17]. It is obviously desirable to employ periodic boundary conditions, which raises an important point. Sometimes the chain intersects itself in one of the image lattices as shown in fig. 2b. Such a chain, bcb’, is called a loop, in contrast to the cycle of fig. 2a. Permissible loops in fig. 1 are O&I and adehfa. The difference between loops and cycles is that cycles do not change the net polarization P of the finite lattice with N = IZ* sites, where P is defined as N
P=
cpj, i=l
whereas loops change the components
Fig. 1. A typical proton configuration in twodimensional square ice. The 3 X 3 lattice has periodic boundary conditions.
of the polaiiza-
tion P, and/or Py by multiples of 21j2 np. In fig. 1 the polarization is P, = 3 (p/21j2) and P,, = - 3 (p/21j2). Each of the loops odgu and ndehfa change PY to +3 (c(/2l/*)i (There are no Py = 0 states for an n X n lattice with n odd.)
332
A. Yanagawa.J.F. Nagle/Correlotio!l functionsfor ice
RS [8] deliberately eliminated loops in their program. Strictly speaking, that was not legitimate since there is no such restriction in eqs. (1) and (2). The reason for this restriction was to save computer time by reducing the number of states to be sampled. This is an unfortunate restriction because it requires G in eq. (2) to be identically zero. Also, in order to obtaing, RS used some plausible, but non-rigorous corrections. Fortunately, it is not necessary for us to make such a restriction of e!iminating loops for twodimensional ice. Thus, we calculate G as well as g. Furthermore, in our program we simultaneously calculate two sets of averages. The first set is for all states including those for which P # 0. Such averages will be designated with subscripts P, e.g. gp.The second set is identical to the averages of RS and includes only those states generated with P = 0. These averages will be designated with subscript 0, e.g.go. With both sets of averages we can test the restriction and the corrections used by RS [8]. Ideally, the lattice should be very large so as to minimize finite size errors. However, increasing size increases the number of states to be sampled and the computer time. This is especially serious because the percentage of loops decreases rapidly with larger lattice size and the system becomes frozen into a particular polarization for long computer times (many state changes). However, for valid comparison to the work of RS it is necessary to have a lattice of “comparable” size. Fortunately, from the point of view of finite size effects a lattice of comparable size in two dimensions is not one with the same number of lattice sites as the one in three-dimensions. Rather, it is one with the same linear dimensions, so an PIX n two-dimensional lattice in two dimensions is “comparable” to an n X n X n lattice in threedimensions. For such comparable lattices the shortest loop is the same length. RS [8] studied a cubic lattice with 4096 sites which has a shortest loop of 16 steps, so we began our work with a 16 X 16 square lattice. For our method which allowed loops this turned out to be the largest lattice for which we could obtain reliable results in a reasonable amount of computer time. For larger lattices the low frequency of loops is a severe problem which induces very long period fluctuations in the averages. However, for 12 X 12 and 8 X 8 lattices the statistical convergence is quite rapid. As was mentioned in the introduction eqs. (1) and (2) are valid for isotropic systems. However, for sys-
terns of lower symmetry, such as hexagonal ice in three dimensions or two-dimensional ice, the correlation functions xi (pl .P~)/P~ must be replaced by the more basic correlations [2,7]
where I’, is the unit vector along the electric field. By symmetry these correlations reduce to those in eqs. (1) and (2) for cubic syst,ems. Twodimensional square ice is an extreme example of asymmetry because two of the six water molecule configurations have no dipole moment. However, this does not qualitatively change the statistical problem so long as the proper orrelations are calculated. One feature induced by this asymmetry which is unlike cubic ice is that the six different water molecule configurations do not occur with equal probability. An exact calculation of Wu [16] shows that each of the water molecule configurations with zero dipole moment has occurrence probability of 19% compared to only 15.5% for the other four configurations. Consider the term in (7) when j = 1, which we will call the self-correlation function g(O). This is the first term in the sum for g. The exact result of Wu requiresg(O) = 0.930 for square ice, compared to g(O)= 1 for cubic ice. In our computer program the correlation functions in (7) are computed for E = E, along the x axis and also for E = I?,, along they axis for all states generated by the Monte Carlo method. By symmetry these correlations are equal for the sum of all states, but not for a subset of states. Thus, one monitor of the reliability of the sampling procedure is the difference between the correlations with the two different field directions.
3. Results of Monte Carlo calculations 3. I. Distribution of chain lengths Fig. 3 shows the probability p(l) of occurrence of cycles and loops of length 1 for a 16 X 16 lattice.Table 1 gives pertinent statistics on the average cycle length, the average loop length and the frequency of.occurrence of loops_ Most of the cycles are just the simplest square cycle of length 4. The fact that the average cycle length is less than 6 means that only a small fraction of the system is affected by a single configura-
A. Yanaga wa, J. F. Nagle/Conelation
0.8-
functions for ice
333
- 0.4
- 0.3
o8 >
-0.2
-Of
I
Fig. 3. Probability of chain Iengthsp(I) versus Iength 1. The solid curves show p (I) (open circles) and p(l) X 10 (solid circles) for cycles with the scale on the left. The dashed curve showsp(0 for loops with the scale on the right. The lattice sizeis16X 16.
tional change and so most of the correlations do not change. Since the computation of the correlations consumes a large fraction of the computer time, it was decided to compute correlations only after k = 25 successive configurational changes. In comparison, RS [S] used k = 200, but this was necessary because, even though the average cycle length is larger for three dimensional ice, there are many more molecules for lattices of comparable size. 3.2. Computation of G by Monte Carlo method G is technically an easy and inexpensive quantity to compute. Consider
(8)
Thus, for E in they direction one only needs to add they-components of the dipole moments in one row of the lattice, multiply by the number of rows, and square the result to pbtain P# for each state generated. Fig. 4 shows computed vaBues of G as a function of the number of configurations sampled for lattice sizes 8 X 8,12 X 12 and 16 X 16. The average over sums of 10000 configurations show typical fluctuations which become larger with larger n. This trend with larger n is to be expected from the occurrence of loops shown in the last column of table 1 because, as was discussed in section 2, the only way to change P ( and therefore from eq. (9) the value of G) is via loop changes. Fig.4 shows that for a 16 X 16 lattice one clearly does not sample as many different P’s in 10000 configurations as for an 8 X 8 lattice. Of course, the cumulative averages over all generated states, also shown in fig. 4, show sinaller fluctuations, akhough long fluctuation waves are clearly apparent for the 16 X 16 lattice. Averaging over the 15,20 and 24 runs, each run involving 10000 configuration changes, for the 8 X 8,12 X 12 and 16 X 16 lattices, respectively, gives our final estimates of G with ENS errors which are presented in the third column of table 2.
3.5 30 PC. Ee
Since the central molecule labelled 1 can be any ofthe N molecules, it is convenient to sum over all of them which gives
(8 x 8) ___a.-.
--
_L
.e-*-
3.0 25 .
3.0
11 116x161
20
I I
5
10
c/10000
% of loops
8X8 12x 12 16x 16
Average loop length
Average cycle length
14.1
15.33
5.42
3.0 0.5
23.83 32.62
5.53 5.53
a-4.
&---
--
.
Lattice size
-
3.5
4.0
Table 1 Distribution of cycles and loops
*
15
20
c.4
,
25
Fig. 4. Monte Carlo resulti for G for three sizes of lattice as a function of number of configurations C sampled. The solid line with sotid circles shows the value of G averaged eve: the preceding 10000 configurations sampled. The dashed linewith crasses shows the value of G averaged over all preceding confgurations sampled.
A. Yanagawa, J.F. Nagie/Correlation functiom for ice
334
Table 2 Monte Carlo results for the polarization function G, the self-correlation functiong (‘1 and the Kirkwood correlation functiong.The subscripts P and 0 refer to all states and states with zero polarization, respectively, with relative numbers n,/n~~ of the two kinds of states btticc
“&P
G
gp
gp
gP
0.33 0.35 0.33
2.88 f 0.04 2.78 c 0.06 2.86 r 0.11
0.942 f 0.002 0.937 I?0.001 0.933? 0.001
0.920 f 0.003 0.927 * 0.002 0.925 f 0.002
1.41 1.40 1.43
size
8X8
12x 12 16 x 16
3.3. Result for g(o) and error analysis
RO
1.52 1.53 1.41
given for all the entries in table 2 are t on2 as obtained by the second method.
The results for g(o) for two different samplingprocedures,&
for only P=G states and $1
for all states,
are given in table2. Although these estimates are in reasonably good agreement with the exact value, g(O) ~0.930, the statistics for a given lattice give standard
deviation errors that are generail smaller than the real error. However, it is ciear that gb) IS steadily decreasing toward the exact value and g&O)is steadiIy increasing as the size of the lattice increases. Doing a simple extrapolation versus l/r1 for these n X n lattices gives us a final Monte Carlo e;tim&e,g(@ = 0.930 -t-0.005. For completeness, let us men&on how the standard deviation errors were corn “ted. In our first method the individual values ofgi( B) for each state obtained (after 25 configuration changes as explained in section (2.1)) were used as the n data points in the usual formu-
la for normal distributions,
02 = m
2 [g,!0’-~q2/,(,-l).
i=l
However, even after 25 configuration
changes it might
be argued that the new state is not randomly independent of the previous one. Therefore, in our second method of obtaining statistics g(O)was first averaged over sequences of 10000 configuraticn changes(m = 400 states = 10000/25) and each of these averages $‘) should be much more independent than the individual gj”) used in the fikmethod. Then these&‘) replace the gi”) in (10) and n becomes of the order of 10 to 20 depending upon the time the computer was run which depended upon the lattice. Just as in the
first method jj”) is the grand mean. This second mkthod gave somewhat larger values of urnIbut by less than a factor of 2, than did the first method. The errors
3.4. Results for the correlation fitnctions and the Kirkwood correlation factor Following RS [8] it is useful to consider the sum of correlation functions over the tih shell of neighbors, defined as.g@. For an infinite system m
g=
c
p). v=o
(11)
The shell-v = 0 just gives g(O)which is the self-correlation and the shell Y = 1 gives g(l) which is the sum of the four nearest neighbor correlations, etc. The maximum coordination shell that we routinely consider is v = 16 for the 16 X 16 lattice. The number of neighborsin shells 0 to 16 is 10 1. In fig. Sa is shown gcV)versus v for the two different sampling procedures, gk) for P = 0 states and gg) for unrestricted P states. Both calculations give values of g(O) which are very close to the exact value 0.930 (see table 2) and both give a negative gc2). Further correlations are rather undistinguished, except that g;) tends to be greater than gt). This becomes apparent when the sum of alI correlations up to the vth shell, defined as ”
g””=
izdi),
02)
is plotted versus v as seen in fig. Sb. Clearly gp” continues to increase as a function of v and gr continues to decrease. Although this may seem alarming, it is expected. As v becomes large the vth coordination shell contains most of the molecules in the finite lattice. Therefore, one expects g? to become approximately
A. Yanagawa, LF. NaglefCorrelation
335
functions for ice
g can be made. The combined estimates are summarized in table 2 for the two casesgp andgo. No errors are given because of the intrinsic uncertainty in the procedure. The number of state changes made was 100000 for the 8 X 8 and for the 12 X 12 and 160000 for the 16 X 16 lattices with correlations calculated after every k = 25 states. The fact that estimates ofg made from bothgg and from our new gg are in substantial agreement as seen in fig. 5c lends support to the correction scheme of RS [8].
4. Series expansion calculations of the correlation functions The general method of weak graph series expanFig. 5. Monte Carlo results for correlation functions for a 16 X 16 lattice from a sample of 10000 configurationsversus shell number V. (a) The shell correlation functionsg$?) (solid line and circles) andgg) (dashed line and crosses). (b) The summed correlation functionsgy (solid line and circles) and gr (dashed line and crosses). (c) The corrected correlation functionsggi (solid line and circles) andggy (dashed line and crosses).
sions [M] has been applied to the correlation functions in cubic ice by Gobush and Hoeve [ 191 and Stillinger and Cotter [9] _In this section analogous computations are reported for square ice. The StillingerCotter [9] expansion variable, X = tanhlvfi, where 2w is the energy of a Bjerrum fault, will be used. The basic correlation functions are g’O”(h) = (3//J2) (PO- $T) (I$ - &)A
equal to G and gr to become approximately equal to 0 and we verified this for a few of our calculations. Since there is no clear value of n at which to cut off the summation, this poses the interesting problem of extracting from fig. Sb a reasonable estimate ofg. RS [8] proposed !m obviously non-rigorous correction to gr as follows. The zero P constraint was supposed to introduce a linear negative bias proportional to -l/N for each molecule within the uth shell so that
where p. and 4 are the dipole moments of the 0th and rth molecules. These correlation functions may be obtained most easily from the modified partition function,
Z($-J,~)=
C cnrv@
states
(13)
where NO, is the number of molecules in the system which are within the tih shell. It is hoped that gz will then give a better approximation tog for values of u such that 1 - (NW/N) is of the order of l/2. We now propose also to correct gp” in the only way consistent with the aforementioned linear RS correction, namely gp” =g;; (1 -A&/N)
+ G(N,,/N),
(16)
zo2,,
where the states have n Bjerrum faults and where ZO~~=~~P~~~@O.~E~IPOI)+~(J~,
gr =gg (1 -N,,JN),
(15)
~&/ly-I)]
(17)
by the formula a2 m Z(~o,o;) g(O” (1) = 3 -Qd-&oe=o=‘yr
(18)
Using the weak graph method eq. (16) becomes
(14)
where it is now supposed that each molecule within the uth shell adds a positive linear bias (G-gg)/N. Fig. SC shows gz and gE$ versus v. Both g$ and g$ now have reasonably flat portions from which estimates for
as is thoroughly explained in the literature [14,18]. The expression (19) is then expanded term by term where each term corresponds to a subgraph of the square lattice. After performing the operation in (18)
A. Yanagatvu. J.F. NaglejCorrelation functions for ice
336
i = vi
=
la
lb
-JF
fi
2a -1
Zb -I
3a
3b
~3E-3/z
4 9
Fig. 6. Vertex types i and vertex weights “i in the series expansions. one has a sum over the subgraphs of the lattice
g@‘)(X)= 5 w(~),
(20)
where w(((,) is the contribution or weight of subgraph ai. It is given by a simple product W(@)‘(-Z)’
n Vi(U),
vertices in U
(21)
where t =X/3, e is the number of edges in the graph aand ui((g) is the weight of the ith vertex of &These vertex weights vi(@) depend only upon the local coordination of the ith vertex in the subgraph @. Whenthe electric field is in the vertical direction, there are only seven different types of local coordination that have non-zero vertex weights and these are shown in fig. 6 with their weights. Types la, lb, 3a and 3b may only occur on vertices 0 and r.
The correlation function g(O)can be treated in a similar way except that r and 0 are the same vertex. In this case the vertex weights for the special vertex 0 are ~2, = 1, u’& = -2 and u4 =O. There is also a contribution w(0) = 1 from the null subgraph. Having found the graph weights w(($) in (20) the calculation reduces to counting the subgraphs Qj of the different types and details will not be given. Results are expressed for gc”) which is the sum over alI g@) where r is in the tih shell around 0. We have obtained g(O) = 1 -4z4 -68~~ - 16~1~ - 1384z12, g(11=4z-8z3-32z5-8Clz7-648zg,
x
g(*) = -16~~ -32~~ - 176~8 -52&O, gt3) = 4z* - 8z4 - 64z6 t 8z8 - 744z’O, g(4) = 8,s -40~5 - 72~’ - 40829, gts) = 8z4 -3226 - 80~8 - 145zfi,
(22)
where z =X/3, to the highest order shown. Unlike the threedimensional cubic ice case [9] it seems that one of the correlations, gt3) has a higher order term, +8z8,. which is not negative, though its magnitude is very small. In general, hcwever, the higher order terms are negative as in the case for cubic ice. An independent check on these series is obtained by adding them algebraically to obtain the series G(A). This latter series has been obtained using a renormalized weak graph method [14] which requires less graph counting and is therefore less subject to errors. From
1141 C(X) = 1 t 4z t 4z2 t 4z3 - 12z5 -20~~ - 12z7 t 16z8 -(32/3)zg
+ (704/3)z10 f 9762” •r-. . . .
(23)
In addition to theg@‘) shown in (22), all other terms up to and including z7 were found in the remaining correlation functions and the sums for each order agree with (23) up to and including 7th order. Numerical estimates of the correlation functions may be obtained in a number of ways. A particularly illuminating procedure is shown in fig. 7 for ideal ice (X= 1). Defimeg$) to be the partial sum of the terms up to and including the zn term in the series for gt”). In fig.? &‘I for zj=O and 1 are plotted versus l/n. For v = 0 fig. 7 shows that gA”) appears to be converging rather well in the limit l/,t + 0 to the exact value g(O)= 0.930 although our best extrapolation slightly overestimates the exact value by 0.007. The Monte Carlo estimate &) is also shown in fig. 7. For v = 1 ii 7 shows that g I) has much more variation than did ’ value as l/n + 0 appears in this gn8) , but the limiting case to be consistent with the Monte Carlo estimate. For v = 2 to 5 similar plots (not shown) suggest limiting values of g(y) which are very close to, although slightly larger than, the Monte Carlo results. Other ways to estimateg(V) involve examining the ratios of successive terms or the use of Pade approximants PO]. The ratios of successive terms alternate between large and small, so one has to attempt to extrapolate two very short subsequences. Also, the series are a bit short for the full power of the Pade approximant to be high!y beneficial. Therefore, these techniques do not improve much upon the simple graphical technique illustrated in fig. 1. Our numerical estimates for the individual correlation functions g@)are collected in table 3. The values
A. Yanagarua. J.F. NaglefCorrelation
functions
337
for ice
3.0
2.0
1 .o
0
0.2
0.4
0.6
0.8
1.0
Fig. 8. G(k) andggE(X) versus h from series expansions. The number v of correlations summed is indicated for eachgz&Q
MC cases. The last figure in each column g$ is an approximation for the Kirkwood correlation function
Fig. 7. Analysis of the convergence of correlation function series using a l/n pl0t.g:) IS . the partial sum of the first n terms in the series for ideal ice (h = 1). The arrow shows the exact value ofg(O) and the error bars show the range of values given by the hlonte Carlo estimates.
& are extrapolated from the series expansions. The vaiues g(v) ,lma?: are obtained by adding the available terms in the series; this gives an upper bound because the higher terms in the series are redominantly negative. For comparison the valuesg&P obtained from the Monte Carlo method are also given. In table 3 the sum of g'$from n = 0 to II = u are also given in the columns labelled g,“, where x refers to.the nmax, ex or
Table 3 Results of series and Monte Carlo methods. LJis the shell number,g@) give the sum of correlation functions in each shell and guV give the sum of correlation functions over shells up to and including the vth shell. The subscripts nmax, ex and hfC refer to truncated series approximations, sxtrapolated series approximations and Monte Carlo approximations, respectively. All results are for the limit X = 1
g(l). It is most instructive to consider the series approximation as a function of X. Fig. 8 shows G(A) and gg(X) for v = 1 to 5 versus X. The curve for G(X) proceeds monotonically as a function of X to the e.xact value G = 9/n = 2.865 at h = 1. Earlier estimates using three different extrapolation methods [14] gave G(1) = 2.848,2.855 and 2.898 which are closely clustered around the exact value. In contrast, the estimated errors for g::(X) shown in fig. 8 only for&(X) are considerably greater than for G(X). Most important is the behavior of g’$(‘R) as a function of h. As can be seen in fig. 8 each gg(h) has a maximum for h < 1 which becomes larger with increasing Y and whose posiiion shifts to larger h. We interpret this as a precursor of a discontinuous decrease ing(A) at X= 1. For example, since g::(X) is a truncated approximation to g(X) and since the individual correlation functions g::(X) are likely to be continuous functions of X,gi:(X) should be continuous, A continuous approximation to a discontinuous g(X) would have precisely the behavior shown by g’$(h) in fig. 8 because in order to simulate the discontinuous decrease at h= 1, the continuous approximation would reach a maximum for X < 1 and then decrease.
av
gex 0
0.937
0.937
0.930
0.937
0.937
1 0.836 2 -0.277 3 0.247 4 0.078 5 0.018
0.681 -0.287 0.195 0.058 -0.010
0.768 -0.318 0.128 -0.064 -0.026
1.77 1.50 1.74 1.82 1.84
1.62 1.33 1.53 i.58
1.57
0.930 1.68 1.36 1.49 1.43 1.40
5. Conclusions First, we have tested two approximate methods, the Monte Carlo method and the series method, against some exact resuits for correlation functions in twodimensional square ice. The series method gives very
A. Yanagawa, J.F. Nag/e/Correlation
338
good estimates [14] of the polarization factor G which equals 9/n exactly and of the self-correlation function g(O)(see table 3) which equals O.Y30 exactly. The Monte Carlo method is less accurate for G (see table 2) than the series method, but the standard deviation errors include the exact result. For g(n) the Monte Carlo method is even more accurate than the series method, although not much effort was expended to extend the length of the series.
Second, we have used both the Monte Carlo and the series methods to calculate the Kirkwood correlation factor g for which there is no exact result. These estimates of g are considerably poorer than our estimates of G and g(o). Part of the reason for the relative inaccuracy in the estimates ofg is intrinsic to the calculation methods. In the Monte Carlo method the finite Iattice size imposes a finite, but arbitrary, cut-off in summing the correlation functions. It also requires an ad hoc correction such as the one introduced by RS [8], but our inclusibn of states with non-zero polarization in the calculation lends support to this correction. In the series method the calculation of g requires not only series extrapolations, similar to those in the G or g(O)calculations, but also an additional summation over correlation functicns. Although these considerations make it difficult to estimate g and especially difficult to assign errors, both methods clearly indicate that g(l) is much smaller than G(1) with even more certainty than for threedimensional ice. From tables 2 and 3 the Monte Carlo method appears to give g(1) about 1.4 and the series method appears to give g(l) about 1.6, we estimate that g(l)=
1.5kO.3
(24)
for square ice. We can now conclude that the correlation func. tions in twodimensional square ice behave similarly to those in the threedimensional ices, supporting the assumption of dimensional similarity discussed in the introduction. In particular, the previous conclusion that eqs. (4) and (5) held in three-dimensional ice 17, 141 now hold at least as firmly in two dimensions. Therefore, either g(A) and/or G(X) is discontirmous at X = 1. The fact that both the Monte Carlo method and the series method accurately reproduce the exact value of G(1) and that the series give a smooth monotonic curve for G(1) suggests that G(X) is not strongly discontinuous at h = 1. In contrast, the relatively in-
fitnctions for ice
accurate estimates of g(1) may not be due entirely to intrinsic computational problems. Instead, they are likely to be due in part to g(h) being discontinuous and the difficulties in approximating such a function with continuous approximations. This conclusion is strongly indicated by the series expansion method which would be expected to approximate a discontinuousg(h) by the sequence of approximationsgi?(X) shown in fig. 8. Now that the approximate methods for correlation functions have been tested against exact results in twodimensional square ice and now that the resuits of the calculation have been shown to be similar to the re-
sults for three-dimensional ice and to indicate that has a discontinuous decrease at h = 1, the previous conclusion [14], that for real three-dimensional ice with a small concentration of Bjerrum defects
g(h)
g(X) = G(A) = 3,
(25)
can be drawn.even more strongly.
Acknowledgement
We wish to thank Professors J.S. Langer and R.B. Griffiths for considerable discussionand R.B. Griffiths for his ergodicity proof which was instrumental in excisinga bias from an earlier version of the Monte Carlo program. This research was supported by the United States National Science Foundation grant DMR 78-10226.
Appendix: Long range correlations and shape dependence
If the correlation functions behave as (JII .I+> a f(0) r,F*-“, then for a system with periodic boundary conditions, G-g=
lim
‘K *-
lim
c
Vdrn
rj>rK
f(e) r,+n
2K rw z
15, ‘K
>i_
$
s
f(0)r-2-~rdrdB,
(26)
’ ‘K
is clearly zero for 7 > 0 or for correlations given by
eq. (3); For ideal ice n = 0, so in order for the integral
A. Yarugawa, J.F. Nagle/Correlation fiolctions for ice in (26) to remain finite, there must be strong cancellations in the 0 integral, which in fact occur in Sutherland’s formula [12]. However, this raises the .jnteresting question whether G depends upon the shape r,(e) of the sample. There are four points to be made: (1) The exact calculation of G involves the use of the transfer matrix for rectangular n X m shapes [16]. Taking the largest eigenvalue corresponds to letting m + m before n -+ ~0, so this method might be expected to give G for long thin rectangles, although the final formula is isotropic. The series method never refers to shape of the sample. The fact that the Monte Carlo result for G for square n X n shapes agrees so
well with the other two methods is inconsistent with a shape dependence unless it can be argued that the exact method and the series method implicitly, and mysteriously, require an n X n lattice. (2) We performed some Monte Carlo calculations for n Xm lattices, jz >m. For the field along the short m side, for which one obtains the best statistics because loop changes are more frequent, the value of G was 9/n within statistical error. However, for the field along the long IZside, for which the statistics were’ poorer,. the value of G decreased as n/m increased, reaching G = 1.85 +0.26 for n/m = 16/6, although we
do not have great confidence in this last result. Experimentally, fields are usually applied along the short dimension for which our calculations indicate G to be independent of n/m. (3) The asymptotic expression for the correlations used in (26) is only expected to hold when thejth molecule is far away from any surface, and surface corrections might conspire to negate any shape dependence. (4) One does not expect any shape dependence in the theoretical dielectric constant in either the Kirkwood-Frohlich theory or in the Onsager-Slater theory. In particular, the usual macroscopic depolarization factors have been taken into account by using
339
the internalE field in these theories. After consideration of these points we are inclined to discount shape dependence in G.
References [l] J.F. Nagle,Chem. Phys. 43 (1979) 317, preceding paper. [2] J.G. Kirkwood, J. Chem. Phys. 7 (1939) 9 11. [3] H. Frohlich, Theory of dielectrics. Dielectric constant and dielectric loss (Clarcndon Press, Oxford, 1949). [4] L. Onsager and M. Cupuis, in: Electrolytes, The electri cal properties of ice, ed. B. Pesce (Pergamon Press, New York, 1962) pp. 27-46. [5] L. Onsager and L.K. Runnels, J. Chem. Phys. 50 (1969) 1089. 16) J.C. Slater, J. Chem. Phys. 9 (1941) 16. [7] J.F. Nagle, I. Glaciology 21 (1978) 73. [8] A. Rahman and F.H. Stillinger, J. Chem. Phys. 57 (1972) 4009. [9] F.H. Stillinger and M.A. Cotter, J. Chem. Phys. 58 (1973) 2532. [IO] L. Pauling, J. Am. Chem. Sot. 57 (1935) 2680. [ll] J.D. Bernal and R.H. Fowler, J. Chem. Phys. 1 (1933) 515. [12] B. Sutherland, Phys. Letters 26A (1968) 532. [I31 F. Rys, Helv. Phys. Acta 36 (1963) 537. [14] J.F. Nagle, J. Chem. Phys. 61 (1974) 883. [15] B. Sutherland, C.N. Yang and C.P. Yang, Phys. Rev. Letters IV (1967) 588. [16] E.H. Lieb and F.Y. Wu, in: Phase transitions and critical phenomena, Vol. 1, eds. C. Domb and M.S. Green (Academic Press, London, 1972) ch. 8. See especiaUy p. 450 forg(O)(l). [I71 A. Yanagawa, Ph.D. Dissertation, Carnegie-Mellon University, USA (1979). [IS] J.F. Nagle, in: Phase transitions and critical phenomena, Vol. 3, eds. C. Domb and M.S. Green (Aczxiemic Press, London, 1974) p. 653. [19] W. Gobush and C.A.J. Hoeve, J. Chem. Phys. 57 (1972) 3416. 1201 D.S. Gaunt and A.J. Guttmann, in: Phase transitions and critical phenomena, Vol. 3, cds. C. Domb and h1.S. Green (Academic Press, London, 1974) ch. 4, p. 181.