J. Phys. Chem. 1996, 100, 1197-1205
1197
Solution of the Ornstein-Zernike Equation for Spheres with Octupolar Surface Adhesion: Toward a Simple Model of Water† L. Blum* Department of Physics, P.O. Box 23343, UniVersity of Puerto Rico, Rio Piedras, PR 00931-3343
F. Vericat Instituto de Fı´sica de Lı´quidos y Sistemas Biolo´ gicos (IFLYSIB), c.c 565 (1900), and Departamento de Fı´sicomatematicas, Facultad de Ingenierı´a, UniVersidad Nacional de La Plata, La Plata, Argentina ReceiVed: March 22, 1995X
An analytic solution of the molecular Ornstein-Zernike equation for spheres with octupolar sticky potentials is given explicitly. In its most general form the closure of the direct correlation function is of the form of the mean spherical approximation for arbitrary multipolar interactions, and the total correlation function contains terms that arise in the Percus-Yevick approximation for spheres with anisotropic surface adhesion. In addition to generalizing several earlier analyses of special cases of this closure, the solution presented here contains new simplifying insights that reduce the complexity of the resulting algebraic equations. The tetrahedral octupole case can be explicitly solved in terms of the inverses of two 3 by 3 matrices. We give explicit solution for a model that has the nearest-neighbor geometry of water and show that the atom-atom pair correlation functions are in rather fair agreement with the neutron scattering experiments, considering the shortcomings of the sticky potential.
1. Introduction Analytic models with explicit, closed form expressions for thermodynamic and structural properties can be very useful in understanding the nature of some important processes in molecular fluids, such as solvation and solvation dynamics. In previous work we have discussed formal solutions of the anisotropic Ornstein-Zernike (OZ) equation2-7 with sticky and electrostatic closures using the BT invariant expansion formalism.2 Understanding water and its unique properties from a molecular point of view is a challenge.8,9 The water-water intermolecular potential cannot be measured directly, and it has been shown that there are important three and multibody effects that need to be accounted for. The potentials used to calculate the properties of water are thus effective pair potentials. Effective potentials for liquid water consisting of a spherical shell plus 3-7 fixed point charges were extensively studied. Early attempts were made by Bernal and Fowler and Weissmann and Blum (ref 10). The excellent work of Stillinger, Ben-Naim, and Rahman provided much insight and understanding about this problem.11,12 There is today a rather extensive literature of excellent work in this area,13-16 and potentials such as the TIPS and especially SPC/E reproduce very well with the known properties of water. The class of potentials that can be solved analytically consists of a spherical hard repulsive exclusion shell plus a potential expanded in spherical harmonics.2 The simplest model is the hard spheres with point dipoles, first solved analytically by Wertheim17 in the mean spherical approximation (MSA). A model which includes a permanent point dipole, polarizability, † It is a great pleasure to contribute to this issue in honor of Harold Friedman. Harold is one of the founding fathers of modern solution theory. The work that we are presenting here has in its precursors a seminal paper of Harold and Don Jepsen,1 in which the ion-dipole interactions were discussed for the first time. X Abstract published in AdVance ACS Abstracts, November 15, 1995.
0022-3654/96/20100-1197$12.00/0
and a structure breaking or forming parameter explains very well the dielectric behavior of water.18 More advanced models which include quadrupoles were investigated by Carnie, Chan, and Walker19 and by Carnie and Patey,20 by numerical solution of the hypernetted chain equation and related closures. An analytically solvable model of water consisting of hard spheres with a point dipole and a narrow attractive potential well with the symmetry of a tetrahedral octupole was proposed sometime ago by Bratko, Blum, and Luzar21 and studied by Monte Carlo. In a recent communication we have shown that this very simple model agrees surprisingly well with the experimental pair correlation functions of Soper and Phillips.6,22 The softened version (with a parabolic tetrahedral attractive well and a spherical repulsive core) of the potential discussed here was studied by Monte Carlo simulation (Degreve and Blum, ref 23). The agreement for all three correlation functions (OO, OH, and OH) is just short of spectacular. This model consists of hard spheres with and embedded dipole and a tetrahedral octupolar attractive well. The angular part is of the form
ν(θ, φ) ) A(r) cos θ sin2 θ cos 2φ Blum, Cummings, and Bratko5 solved formally this model in the sticky limit of an infinitely narrow potential well. We present here a simplified solution of this model, which requires only three structural parameters. As we discussed already in our first papers,2,3 the pair correlation functions derived from diffraction experiments depend only on the axisymmetric irreducible representation that with χ ) 0. Since in the sticky potential PY approximation the irreducible representations are not coupled, for the neutron diffraction atom-atom pair correlation functions we need not consider the other representations. This is an important simplification, since it cuts the number of parameters from 6 to 3. In the next section the basic formalism is summarized. In section 3 we discuss the Baxter-Wertheim factorization and © 1996 American Chemical Society
1198 J. Phys. Chem., Vol. 100, No. 4, 1996
Blum and Vericat
derive the full set of closure conditions, for the first time. In section 4 we study the tetrahedral octupole and find also the atom-atom pair correlation functions.
core exclusion condition and the sticky interactions lead to the following conditions for the invariant coefficients of the correlation function mnl mnl hhµν (r) ) hˆ µν (r), r > σ
(7)
We study the Ornstein-Zernike (OZ) equation for a molecular fluid
mnl mnl hhµν (r) ) hˆ µν (r) ) 0, r < σ, m, n, l > 0
(8)
h(12) ) c(12) + F ∫ d3h(13) c(32)
hh000 ˆ 000 00 (r) ) h 00 (r) ) -1, r < σ
(9)
2. Basic Formalism
(1)
where h(12) is the molecular total correlation function and c(12) is the molecular direct correlation function, F is the number density of the molecules, and i ) 1, 2 symbolizes the position B i ) (ψi, θi,φi) of molecule i ) 1, 2. The b ri and orientation Ω approximations that will be discussed here are based on Baxter’s sticky potential,24 in which the pair correlation function is of the form -
h(12) ) hh(12) + Λ h (12)ξ(r12 - σ )
(2)
where hh(12) and Λ h (12) are functions of the orientation of b1 - b r2|, and σ is the distance of molecules 1 and 2, r12 ) |r closest approach of the two molecules. ξ(x) is a singular function related to Dirac’s delta δ(x). The function hh(12) is nonsingular: that is, it has no delta functions and at worst has step discontinuities. While the discussion of the relation of the function Λ h (12) to acutal intermolecular sticky potentials has been dealt with in a brilliant series of papers by Wertheim,25 in this paper we will be concerned with this point in detail. We will regard the quantities in Λ h (12) as parameters only whose physical significance will be discussed elsewhere. The inclusion of the solution of the arbitrary multipole mean spherical approximation (MSA) does not require addition work. In the MSA the direct correlation function is
c(12) ) -βu(12), r > σ
mnl mnl Λ h (12) ) Λ h (21) ) ∑λµν (r) Φ ˆ µν (Ω1Ω2Ωr)
(10)
mnl µν
For the direct correlation function the explicit form of the boundary condition is
[
]
1/2 (2l + 1)! 1 mnl (r) ) -β l+1QµmQnν ) cˆ µν (2m + 1)! (2n + 1)! r 1 mnl - l+1µµν , r > σ (11) r
To solve the OZ equation, we need the Fourier transform of the pair correlation functions. If we let f(12) stand for either c(12) or h(12), then the Fourier transform of f(12), ˜f(12), is defined by
r 12 exp(ik B‚b r 12) f(12) ˜f (12) ) ∫db
(12)
After straightforward operations we obtain mnl mnl ˜f (12) ) ∑˜f µν (k) Φ ˆ µν (Ω1Ω2Ωk)
(13)
mnl µν
where k is the absolute value of B k and Ωk its orientation. Here, mnl mnl (r) jl(kr) (k) ) 4πil ∫0 dr r2 ˆf µν ˜f µν ∞
(3)
(14)
where β ) 1/kBT is the Boltzmann thermal factor and u(12) is the electrostatic interaction between the molecules. The key ingredient of the analytical solution process is the BT invariant expansion of the correlation functions.2 The invariant expansion of the functions h(12) and c(12) in terms mnl (Ω1Ω2Ωr) is of the rotational invariants Φ ˆ µν
The invariant expansions (4) and (5) of the correlation functions are independent of the choice of the reference frame. Taking the z-axis along the intermolecular axis, the vector b r12 yields the irreducible representation in which the OZ equation is decomposed in a set of matrix equations.2 We get from eqs 4-6
mnl mnl (r12) Φ ˆ µν (Ω1Ω2Ωr) h(12) ) ∑hˆ µν
(4)
mn mn (k) Φ h µνχ (Ω1Ω2) ˜f (12) ) (-)χ∑˜f µνχ
mnl mnl c(12) ) ∑cˆ µν (r12) Φ ˆ µν (Ω1Ω2Ωr)
(5)
(15)
mn µνχ
mnl µν
mnl µν
where mn m n (Ω1Ω2) ) [(2m + 1)(2n + 1)]1/2Dµχ (Ω1) Dν-χ (Ω2) Φ h µνχ
(16)
where mnl (Ω1Ω2Ωr) ) [(2m + 1)(2n + 1)]1/2 × Φ ˆ µν
∑ (µ′ µ′ν′λ′
)
m n l Dm (Ω ) Dn (Ω ) Dl (Ω ) (6) ν′ λ′ µµ′ 1 νν′ 2 0λ′ r
m (Ω) is the Wigner rotation matrix and In this equation, Dµµ′
(
m n l µ′ ν′ λ′
)
is the usual notation for the 3 - j symbol, Ω1 and Ω2 are the Euler angles that give the orientation of molecules 1 and 2, r12 respectively, and Ωr gives the orientation of the vector b joining the centers of mass of molecules 1 and 2. The hard-
and mn (k) ) (-)χ∑ ˜f µνχ l
(
)
m n l mnl ˜f (k) χ -χ 0 µν
(17)
Substitution of these definitions and expansions into the OZ equation (1) yields the desired irreducible matrix form of the OZ equation2
˜χ ) H ˜ χC ˜χ H ˜χ-C
(18)
where the matrix elements of F ˜χ ) H ˜ χ, C ˜ χ are given by mn F ˜ χ ) Ff˜µνχ (k)
(19)
A Simple Model of Water
J. Phys. Chem., Vol. 100, No. 4, 1996 1199
We convert the Fourier Bessel transform of eq 14 into a onedimensional Fourier transform. Using the integral representation of the spherical Bessel function
jl(kr) ) (1/2)i-l ∫0 dt Pl(t)[eikrt - (-)le-ikrt] 1
(20)
where Pl(t) is the Legendre polynomial of order l. Substituting into eq 14 and also using the transforms (17) and (19), we obtain for the pair correlation function h(12) nm nm nm H ˜ νµχ (k) ) ∫0 [eikrJνµχ (r) + e-ikrJνµχ (r)] dr ∞
nm Sνµχ (r) ) 0
l
(
inf(m,n) mnl (r) ) -(2l + 1) 2πFfˆµν
{
(21)
nm (r1) ∫0∞F νµχ
)
∞ m n l mnl 2πF ∫r dr1 r1Pl(r/r1)hˆ µν (r1) χ -χ 0
(22) Explicitly, we get for r < σ nm Jνµχ (r)
) (-)
χ
∑l (
m n l χ -χ 0
)∑
-
nm Λνµχ θ(r
-
-σ )
i)0
where ali is the coefficient of the Legendre polynomial which is nonzero and given by
(l + i)! 2 i! [(l - i)/2]! [(l + i)/2]! l
l
(
)
m n l mnl λ χ -χ 0 µν
(25)
mnl are defined by The parameters bµν,i
mnl mnl ) 2πF ∫σ dr1 r1-i ˆ µν (r1) bµν,i 1 h ∞
nm mn;i i nm Jνµχ (r) ) ∑Bµν,χ r - Λνµχ θ(r - σ-)
)
m n l l mnl ab χ -χ 0 i µν,i
nm nm nm (k) ) ∫0 [eikr Sνµχ (r) + e-ikr Sνµχ (r)] dr C ˜ νµχ
l
(
(33)
(34)
where Q ˜ Tχ (-k) is the complex conjugate and transpose of ˜ χ(k) is even Q ˜ χ(k). In essence, it is because the matrix I - C and symmetric as a function of k that it can be factorized into ˜ Tχ (-k). The first matrix the two matrices I - Q ˜ χ(k) and I - Q is nonsingular in the upper half complex k-plane, while the second is nonsingular in the lower half complex k-plane. From eq 31 it can be shown3-5 that the factored correlation functions must be of the form
˜ χ(r) Q ˜ χ(k) ) I - ∫0 dr eikrQ σ
(35)
Fourier inversion of eq 34 leads to
(27)
Sχ(r) ) Qχ(r) - ∫r dr1 Qχ(r1) QTχ (r1-r) σ
(28)
˜ χ(k)] ) [I - Q ˜ Tχ (k)]-1 [I - H ˜ χ(k)][I - Q (29)
σ
)
(30)
(38)
This equation must be used to determine the functions Qχ(r). Taking the lm successive derivatives of eq 38 and using eqs 23 and 40, we find that Qχ(r) must be of the form
Qχ(r) ) Q0χ(r) + Λχθ(σ- - r), 0 e r e σ The MSA condition on the direct correlation function, eq 3, yields
(37)
where by Cauchy’s theorem the right-hand side makes no contribution, so that
Jχ(r) ) Qχ(r) + ∫0 dr1 Jχ(r-r1) Qχ(r1)
∞ m n l mnl 2πF ∫r dr1 r1Pl(r/r1)cˆ µν (r1) χ -χ 0
(36)
where QTχ (r) is the transpose of Qχ(r). A key step in the BW factorization is the Fourier inversion of
where the matrix S has elements nm (r) ) (-)χ∑ Sνµχ
}
[I - C ˜ χ(k)] ) [I - Q ˜ χ(k)][I - Q ˜ Tχ (k)]
Similarly, the direct correlation function yields ∞
()
where I is the identity matrix. The second term in this equation can be factorized
(26)
with
l
1 1 δ′(r - r1) - P′l(1)δ(r - r1) + r r2 r1 1 P′′l θ(r - r1) dr1 (32) r r3
˜ χ(k)] ) I [I + H ˜ χ(k)][I - C
i)0
(
)
m n l (-)χ × χ -χ 0
In this section, we describe the Baxter-Wertheim26,27 (BW) factorization of the molecular OZ equation. Following earlier work2,5 the OZ equation (18) can be rewritten as
Equation 23 can be rearranged to
mn;i Bµν,χ ) (-)χ∑
χ)-inf(m,n)
(
mnl as We find explicit expressions for the coefficients bµν,i mnl mnl functions of the parameters µµν and λµν in the next section.
(24)
for l - i g 0 and even. For other values of l and i, ali ) 0. The elements of the matrix Λχ are given by nm Λνµχ ) 2πF(-)χ∑
∑
3. Baxter-Wertheim Factorization of the OZ Equation mnl i ailbµν,i r
(23)
ail ) (-)(l-i)/2
(31)
From the orthogonality properties of the Legendre polynomials, we invert the integral equations (21) and (29). Using the mnl mnl mnl generic notation ˆf µν (r) to represent either hˆ µν (r) or cˆ µν (r) and nm nm nm F νµχ(r) to represent Jνµχ(r) or Sνµχ(r), we have
where the matrix J has elements nm Jνµχ (r) ) (-)χ∑
for r > σ
with
(39)
1200 J. Phys. Chem., Vol. 100, No. 4, 1996
Blum and Vericat
l
m 1 Q0χ(r) ) ∑ qi,χ(ri - σi), 0 e r e σ i)1 i!
The equations for the coefficient matrices qi,χ are obtained by, for example, by setting r ) 0 in all the derivatives. We have (n) (n) J(n) χ (r) ) Qχ (r) + ∫0 dr1 Jχ (-r1) Qχ(r1) σ
Bni,χ ) (-)i-n
(40)
Bni,χ ) 0, i < n For lm ) 3 (our model of water) we get
[
(41)
B1,χ -B2,χ B3,χ Bχ ) 2B2,χ -3B3,χ 0 3B3,χ 0 0
which gives
nBn,χ )
1
qn,χ + (n - 1)! i!
(-)i-nBn,χ × ∑ (i - n + 1)! (n - 1)! i)n
1
[q0] ) B0[I0 - H0B0]-1{[I0] - [Λ0]}
-
σj+1
i
∑
1
2lm
σj+1
2lm
qi,χ
, j ) 0, ..., lm (43) i+j
1
(44)
This can be written in matrix notation as
n[Kχ] ) -Hχn[qχ]
[ ]
1 0 n) . .
0 . . 2 0 . 0 3 . . . . K0,χ K1,χ [Kχ] ) . . K(lm-1),χ
[ ] []
q1,χ q2,χ/2! qχ ) . . qlm,χ/lm!
(45)
{
)nm
-
2l + 1
inf(m,n)
2πF
χ)-inf(m,n)
1 ∂Sνµχ(r) r
∂r
∑
(
)
m n l (-)χ × χ -χ 0
() }
r1 1 r nm 1 nm dr1 (r) + ∫0 Sνµχ (r1) P′′l - P′l(1) Sνµχ r r2 r3
(54) (47)
Since from eq 31 we know that Sχ(r) ) 0, r > σ, and from the definitions of Sχ and Kj we get, in matrix form
Sjχ ) ∫0 rjSχ(r) ) σ
(48)
( ) | ∂ j S (k) ∂(ik) χ
jgl-2
ik)0
j
j!
s)0
s! (j - s)!
) Kjχ + (-)j(Kjχ)T - ∑(-)s
s T (55) Kj-s χ (Kχ)
where we recall that
and Hχ is a Hilbert matrix. That is,
Hij,χ ) 1/(i + j)
(53)
Equation 51 is used in the equations for the boundary conditions discussed in the next subsection. 3.1. Boundary Conditions. The algebraic equations for the mnl coefficients bµν,i , are obtained inverting eq 30. For r > σ, this yields
mnl (r) cˆ µν
(46)
Λχ Λχ [Λχ] ) . . Λχ
Iχ 0 [I0] ) . . 0
1
1 qi,χ ∑ j + 1 i)1 (i - 1)! i + j
where we have defined
[] [ ]
25
j + 1 i)1 (i - 1)!
Kj,χ ) -
(52)
Here we use the definitions from our previous work5 and eq
Kj,χ ) ∫0 dr r ∑ qi,χ(r - σ ) ) i)1 i! i
(51)
and I0 is the appropriate unit matrix (6 by 6 in our present case).
where we have defined the moments of the factor distribution function 2lm
(50)
{I0 - nKχ - Λχ} ) [I0 - HBχ]-1{I0 - Λχ}
[(i - n + 1)K(i-n),χ - σi-n+1Λχ] (42)
j
]
The dimensions of the submatrices Bi,χ are different for each value of χ. Using this quantity, eq 42 yields the formal solution
lm
σ
i! B , ign (i - n + 1)! (n - 1)! i,χ
(49)
and the size of Hχ depends on the index χ of the irreducible representation. We can now write eq 42 in a condensed matrix form. To do so, we need to define a supermatrix Bχ with submatrix elements
Kjχ ) ∫0 rjQχ(r) σ
(56)
mn,j . and the matrix Sjχ has elements Sµνχ Since the invariant expansion coefficients of c(12), eq 11, for r > σ must be equal to those of eq 54, we require
A Simple Model of Water
J. Phys. Chem., Vol. 100, No. 4, 1996 1201
mnl µµν δjl )
2l + 1 l
j(j-1)aj
2πF
inf(m,n)
∑
χ)-inf(m,n)
(
b1 )
)
m n l mn,l-2 (-)χSµνχ χ -χ 0
lg2
mnl cˆ µν (0)
) finite
(58)
∞
33 (r) ) J22,0
2πF∫σ dr1 r1
Q′χ(0) + [QTχ ]′(0) ) Qχ(σ) QTχ (σ) - Qχ(0) QTχ (0) (59)
4. Molecules with Surface Adhesion of Octupolar Symmetry
u(1, 2) ) uHS(1, 2) + uD(1, 2) + uS(1, 2)
{
)
)b(33) 2
x
3 332 b22,3 x 2 7
(70)
332 ) 2πF∫σ dr1 r1-1hˆ 332 b22,2 22 (r1) ∞
(71)
(61) 332 Λ33 0 ) 4πFλ /x105
30 033 Λ03 0 ) -Λ0 ) -2πFλ /x7
(62) and
Λχ)0 )
2
10 µ Φ ˆ 112(Ω1Ω2ΩR) 3 R 3
(69)
and from eq 26
The dipolar part is
uD(1, 2) )
]}
where b0 is irrelevant and5
000 Λ00 0 ) 2πF
R12 < σ R12 > σ
)[
(
3 3 0 330 3 3 2 r2 hˆ 22 (r1) + 3 2-1 0 0 0 0 0 0 r 1
2 33 ) b0r + b(33) 2 r -Λ
where the hard-core contribution is
∞ 0
(68)
are to be Also η ) πFσ3/6. The parameters b1, b3, and b(33) 2 determined. The elements of the matrix Λ0 are related to the original sticky parameters:
The interaction potential is of the form
uHS(1, 2) )
{(
(60)
Equations 55, 59, and 60 are a set of coupled nonlinear mnl and represents the algebraic equations for the parameters bµν,i complete, formal analytic solution of the molecular OZ equation.
(67)
Similarly,
∞
Qχ(0) ) QTχ (0)
033 b02,3
2x7
033 b02,i ) 2πF∫σ dr1 r11-ihˆ 033 02 (r1)
which from eq 36 means that
Similarly, for l ) 1 (m + n ) odd) we need to cancel the 1/r2 term in eq 54 which from eq 36 means that the BW functions satisfy the symmetry condition
3
b3 ) -
033 b02,1
and from eq 26
(57) For l ) 0 (m + n ) even) we need to cancel the 1/r term in eq 54
5 2x7
(63)
12
and uS(1, 2) is the tetrahedral octupole surface adhesion. The correlation function h(1, 2) is
h(1, 2) ) -1 + λ000δ(R12 - σ-) +
[
Λ00 0
Λ03 0
(33) -Λ03 0 Λ0
(72)
]
(73)
which in our present work are treated as free parameters. In the full theory they will be derived from a closure relation, for example the PY closure,28 or others. In matrix form, we get H0 eq 33 and eq 21
H0(r) ) -B1 - 2B2r - 3B3R2 + Λ0δ(r - σ-)
ˆ 033 + λ303Φ ˆ 303]ξ3(R12 - σ-) + [λ033Φ
(74)
6
λ33lΦ ˆ 33lξl(R12 - σ-), ∑ l)0
R12 < σ (64)
where δ(x) is Dirac’s delta and ξl(x) is a singular function, related to δ(x). The explicit angular dependence of the rotational invariants2 has been omitted for simplicity. In the PY-MSA closure
c(1, 2) ) -βuD(1, 2), R12 > σ
(65)
From eqs 23-28 03 (r) ) J02,0
( ) { 0 3 3 2πF 0 0 0
[
]
3
1
) b1r + b3r - Λ 3
where5
}
033 ∫σ∞dr1 r1 5rr 3 - 3rr1 hˆ 033 02 (r1) - λ02 03
(66)
where the matrices B1, B2, and B3 are
[
]
[
]
[
-6η 0 b b3 0 0 B1 ) -b 1 ; B2 ) (33) ; B3 ) -b 0 0 b 1 3 0 2
]
(75) In the full theory they will be derived from a closure relation, for example the PY closure,28 or others. In the simple problem we need to consider only the χ ) 0. We will therefore drop this subindex in the notation. Notice that from eq 27
H0(r) ) -
∂J0(r) ∂r
(76)
With this result and using the procedure described in the previous section we solve eq 38 to get the coefficients of Q0(r) as functions of b1, b3, and b(33) as well as the sticky coef2 ficients Λ. The results are expressed by the relations (see also
1202 J. Phys. Chem., Vol. 100, No. 4, 1996 the Appendix)
[
] [ ]
B1 -B2 B3 B0 ) 2B2 -3B3 0 3B3 0 0
[ ]
q1;0 [q0] ) q2;0/2 q3;0/6
K0;0 [K0] ) K1;0 K2;0
[
(77)
(78)
1 - 4η 3η 0 1 + 2η 0 Aη ) [I - H3P] ) -3η -12/5η 3/2η 1
-1
{[I0] - n[K0] - [Λ0]} ) [I0 - H0B0] {[I0] - [Λ0]}
[
1/2 1/3 1/4 H3 ) 1/3 1/4 1/5 1/4 1/5 1/6 Similarly
(79) and [q0] can be obtained in the form
Aβ ) [I - H3Pβ] )
[q0] ) B0{[I0] - n[K0] - [Λ0]} ) B0[I0 - H0B0]-1{[I0] - [Λ0]} (80)
]
Aη0 ) [Aη]-1C0 )
[
[
(82)
]
(83)
0 -6η 0 0 P ) 12η 0 0 0 0
0
-b(33) 2 /2 -2/5b(33) 2
0
1+
2b(33) 2 /6
b(33) 2 /4
0 -b(33) 2 0 0
1
]
[M BT]-1 )
[
(81)
][
-Aη0 (Aη + C0Aβ0)-1 0 I AB0 I (Aβ + C0Aη0)-1 0
[
-b3/2 -b1/2 - 3/4b3(1 - 2η/5) b3(1 - η/4) 1 -b /3(1 + η/2) - 3b /5(1 - η/4) 3b /4 -b3/(1 + η/2) 3 1 3 ∆2 -b /4(1 + 4η/5) - b /2(1 - η/10)2 3b /5(1 + η/8) -b /4(1 + 4η/5) 3 3 1 3
β6 ) 1 - b(33) 2 /6
K0 - KT0 - K1KT0 + K0KT1 ) 0 03 30 00 30 03 33 00 30 03 33 - K1;0 - K1;0 K0;0 - K1;0 K0;0 + K0;0 K1;0 + K0;0 K1;0 ) 0 K1;0
(90) 2.
K0 + KT0 - K0KT0 ) 0 33 30 2 33 2 - [K0;0 ] - [K0;0 ] )0 (91) 2K0;0 3. and condition (60) which requires the symmetry of the BW function at the origin
]
(87)
] (88)
[
In the Appendix we outline the procedure. The parameters b1, b3, b(33) are obtained by solvating the 2 system of nonlinear algebraic equations formed by the three boundary conditions, eq 55 1.
(85)
Using the partitioned matrix method,29 we get
-b3/2 -b1/2 - 3/4b3(1 - b(33) b3(1 - b(33) 2 /15) 2 /24) 1 (33) (33) Aβ0 ) [Aβ] C0 ) 2 -b1/3(1 + b2 /12) - 3b3/5(1 - b2 /24) 3b3/4 -b3/3(1 + b(33) 2 /12) β6 (33) (33) (33) 2 -b1/4(1 + 2b(33) /15) b /2(1 b /60) 3b /5(1 + b /48) -b /4(1 + 2b 2 3 2 3 2 3 2 /15) -1
]
(84)
(86)
∆)1-η and
]
(33) 1 - 2b(33) 2 /3 b2 /2
Pβ ) 2b(33) 0 2 0 0
where where
]
0
MB ) [I0 - H0B0] We first rearrangement this matrix using a similarity transformation, which exchanges rows and columns (see Appendix):
[
[
where n is a matrix with diagonal blocks of the form nij ) δiji. The problem of finding [q0] in terms of b1 and b3 essentially reduces to that of the inverse of the matrix
[
]
b1/2 + 3b3/4 -b3 b3/2 C0 ) b1/3 + 3b3/5 -3b3/4 b3/3 b1/4 + b3/2 -3b3/5 b3/4
As was discussed in the earlier work, we get the matrix [K0] from
Aη C0 MBT ) -C 0 Aβ
Blum and Vericat
] (89)
Q(0) ) QT(0)
(92)
The physical branch of the solutions of these boundary conditions should satisfy the low-density asymptotic conditions. 4.1. Correlation Functions. The pair correlation functions are calculated following the method of earlier work.30,31 We recall from eq 4 and eq 19 that the radial functions hmnl(r) can 2 be written as functionals of the irreducible functions Hmn χ (r). mn All the irreducible functions, like H0 (r), can be arranged in matrix form according to
[
03 H00 0 H0 0 33 Hχ)0 ) H30 0 H0 0 0 0 H11 0
]
(93)
A Simple Model of Water
J. Phys. Chem., Vol. 100, No. 4, 1996 1203
We assume that the dipolar part (given by Wertheim’s MSA solution17) (with m, n ) 1, 1) decouples from the sphericaloctupolar one ((m, n) ) (0, 0); (3, 0); (0, 3) and (3, 3)). Only the upper-left block of the matrices in (eq 93) are then relevant. The irreducible functions Hmn 0 (r) can be written in terms of new functions Gmn 0 (r):
H0(r) ) G0(r) + E0(r)
(94)
where E0(r) is a matrix equal to H0(r) inside the hard core and is analytically continued outside:
E0(r) ) -B1 - 2B2r - 3B3r2 + Λ0δ(r - σ-)
(95)
G0(r) ) 0, r < σ-
(96)
so that
The pair correlation function matrix G0(r) verifies eq 96 for r < σ- and the integral equation30,31
G0(r) ) Q′0(r) + Λ0[Q0(r - σ) + Λ0] +
∫σrdr1 Q0(r-r1) G0(r1),
σ < r < 2σ (97)
and
G0(r) ) Q′0(r) + ∫σ dr1 Q0(r-r1) G0(r1), r > 2σ (98) r
Here Q0(τ) is given by eq 40. The invariant coefficients of the pair correlation functions are then computed using Perram’s method32 to obtain G0n 0 (τ), from eq 97, from where we get the invariant expansion eq 4. The intermolecular atom R-atom β pair correlation function is defined by the average
gRβ(|τRβ|) )
∫dbr 12 dΩ1 dΩ2 g(12) δ(br Rβ - br 12 - Bb R + Bb β)
(99)
where atom R is in molecule 1, atom β in molecule 2, b rRβ is the separation of the two atoms, b r12 is the separation of the centers of the two molecules, B bR is the vector for the position of atom R in molecule 1 and B bβ is the vector for the position of atom R in molecule 1 and B bβ is the vector for the position of atom β in molecule 2. The calculation of the atom-atom pair correlation function is identical to that of the neutron scattering cross section (see ref 2, section IV.A), and need not be repeated here. The result depends only on the χ ) 0 irreducible representation and is nm gRβ(r) ) ∑ ∑ ∑[(2m + 1)(2n + 1)]1/2F νµ;Rβ (r) × iR,jβ m,µ n,ν
C nν(bˆ iR) C µm(bˆ jβ) (100) where the sum is over all atoms iR, jβ of the same kind within each molecule, C nν(bˆ iR) is the Racah modified spherical harmonic34 of the angles bˆ iR ) θiR, φiR. The function χnm(r) is given by
(-)n bR nm F νµ;Rβ (r) ) ∫ dz ∫bβ dz H nm (r - z1 - z2) × 4bRbβ -bR 1 -bβ 2 νµ,χ)0 z1 z2 Pn P (101) bR m bβ
() ()
In the simplified model under consideration the three pair
Figure 1. (a, top) The oxygen-oxygen pair correlation function for 33 the case with Λ03 0 ) Λ0 ) 6 hard-core diameter σ ) 2.76 Å and -3 density F ) 0.0033 Å (F* ) 0.7035). (b, middle) The oxygenhydrogen pair correlation function for the same case. (c, bottom) The hydrogen-hydrogen pair correlation function for the same case.
correlation functions are (we drop unnecessary indices) 00 gOO(r) ) F OO (r)
(102)
00 03 (r) + x7[F OH (r)(C 32(bˆ H1) + C 32(bˆ H2)) gOH(r) ) F OH 30 (r)(C 32(bˆ H1) + C 32(bˆ H2))] (103) F HO 00 03 (r) + x7[F HH (r)(C 32(bˆ H1) gHH(r) ) F HH 30 (r)(C32(bˆ H1) + C 32(bˆ H2))] C 32(bˆ H2)) + F HH (33) (r) 7F HH
∑
i,j)1,2
C 32(bˆ Hi) C 32(bˆ Hj) (104)
Furthermore
bH ) 0.918 Å 3 (bˆ Hj) ) C (2
bO ) 0.0654 Å
x30 cos θ sin2 θ e(2iφ 4
The results are shown in Figure 1a-c for the sticky 33 parameters Λ03 0 ) Λ0 ) 6, and in Figure 2a-c for the sticky 03 33 parameters Λ0 Λ0 ) 8.
1204 J. Phys. Chem., Vol. 100, No. 4, 1996
Blum and Vericat the quadratic the hypernetted chain approximation (QHNC). In their work the coordination number is in good agreement with the experimental one, but the contact value of the averaged correlation functions is several times larger than the maximum of the first peak in the experimental oxygen-oxygen correlations. We have assumed that the coupling of the dipoles to the octupoles is weak, and we have neglected it here. We leave the solution with the full coupling case for the future. 5. Conclusion In this paper, we have simplified the analytic solution of the molecular Ornstein-Zernike equation for the general closure, eqs 7-11 presented in earlier work.3-6 Our results show that, in spite of the shortcomings of the sticky potential, the agreement with experimental structure functions is qualitatively correct. Acknowledgment. Support of this work by CONICET of Argentina, the U.S. National Science Foundation (Grants EPSCOR EHR-91-80725), and the Office of Naval Research is gratefully acknowledged. F.V. is member of CONICET. Appendix. Explicit Formulas for the Octupole-Density Case
(
)
For the simplest octupolar model we have from eq 77
-6η 0 0 3b3
b1
0 -b1
12η B0 ) 0 0 -3b3
0 0 2b(33) 2 3b3 0 0 0
0 0 (33) -b -b2 3 -3b3 0 0 0 0 0 0 0
( ) ( )
The Hilbert matrix is
1
Figure 2. (a, top) The oxygen-oxygen pair correlation function for 33 the case with Λ03 0 ) Λ0 ) 8 hard-core diameter σ ) 2.76 Å and density F ) 0.0033 Å-3 (F* ) 0.7035). (b, middle) The oxygenhydrogen pair correlation function for the same case. (c, bottom) The hydrogen-hydrogen pair correlation function for the same case.
In all of our calculations we have used the hard-core diameter σ ) 2.76 Å and density F 0.0033 Å-3 (F* ) 0.7035). We compare our results with the experimental pair correlations in water.22 The corresponding coordination numbers, defined by the integral up to the first minimum (about 3.5 Å, gives NCN ) 5.1 for both calculations. These numbers must be compared against the experimental value NCN ∼ 4.5 - 5.33 The coordination numbers appear to be rather good, in spite of the simplicity of the model. The agreement of the pair correlation functions is satisfactory, considering the fact that the sticky models have a jump in the pair correlations at r ) 2σ, which is clearly shown for the oxygen-oxygen correlations. This roughness is smoothened by the convolutions with the form factors in both the oxygen-hydrogen and hydrogen-hydrogen correlations, both of which exhibit the two peaks seen in the experiment. The oxygen-oxygen correlation function has a jump at r ) 2σ, but if one smoothens it out, then the shoulder at about 4.1 Å could be considered as being the one seen by Soper22 at 4.5 Å. The contact values appear to be satisfactory also. It is worth noting that Carnie and Patey20 have also applied integral equations to water. In their model, a fluid with a polarizable dipole-tetrahedral quadrupole interaction, they used
/2 0 1 / H0 ) 3 0 1 /4 0
0 1 /2 0 1 /3 0 1 /4
1
/3 0 1 /4 0 1 /5 0
0 1 /3 0 1 /4 0 1 /5
1
/4 0 1 /5 0 1 /6 0
0 1 /4 0 1 /5 0 1 /6
b3 0 0 0 0 0
(105)
(106)
so that we need the inverse of
MB ) [I0 - H0B0] Λ0
Λ(03)
-Λ Λ0
(03)
Λ0 )
(107)
Λ(33)
Λ(03)
-Λ(03) Λ(33) Λ0 Λ(03)
( )
(108)
-Λ(03) Λ(33)
1 0 0 n) 0 0 0
0 1 0 0 0 0
0 0 2 0 0 0
0 0 0 2 0 0
0 0 0 0 3 0
0 0 0 0 0 3
(109)
A Simple Model of Water
()
1 0 0 I0 ) 0 0 0
0 1 0 0 0 0
0 1 0 0 0 0
0 0 1 0 0 0
0 0 0 1 0 0
0 0 0 0 1 0
0 0 0 0 0 1
( ) Λ(03) -Λ0
1 - Λ(33)
Λ(03) -Λ0
-Λ(33)
Λ(03)
-Λ(33)
This matrix is then
-Λ(03)
M BT ) [I - HB]T ) T6[I - HB]TT6
(
which is then written as eq 81.
(113)
(112)
3η 1 + 2η 3/2η b3 3b3/5
0 0 1 -b3/2
0 0 0 1 0 0
0 1 0 0 0 0
0 0 0 0 1 0
0 0 1 0 0 0
0 0 0 0 0 1
(115)
1 0 0 T T6 ) 0 0 0
0 0 1 0 0 0
0 0 0 0 1 0
0 1 0 0 0 0
0 0 0 1 0 0
0 0 0 0 0 1
(116)
b1/2 + 3b3/4 -b3 b1/3 + 3b3/5 -3b3/4 b1/4 + b3/2 -3b3/5
b3/2 b3/3 b3/4
(33) 1 - 2b(33) 2 /3 b2 /2
0
-b3/3 -b(33) 2 -b3/4 2/5b(33) 2
References and Notes (1) Jepsen, D. W.; Friedman, H. L. J. Chem. Phys. 1963, 38, 846. (2) Blum, L.; Torruella, A. J. J. Chem. Phys. 1972, 56, 303. (3) Blum, L. J. Chem. Phys. 1972, 57, 1862. (4) Blum, L. J. Chem. Phys. 1973, 58, 3295. (5) Blum, L.; Cummings, P. T.; Bratko, D. J. Chem. Phys. 1990, 92, 3741. (6) Blum, L.; Vericat, F.; Bratko, D. J. Chem. Phys. 1995, 102, 1461. (7) Blum, L.; Vericat, F. Mol. Phys., in press. (8) Eisenberg, D.; Kauzmann, K. The Structure and Properties of Water; Oxford University Press: London, 1969. (9) Franks, F., Ed. Water, A ComprehensiVe Treatise; Plenum Press: New York, 1974; Vol. I. (10) Bernal, J. D.; Fowler, R. H. J. Chem. Phys. 1933, 1, 515. Weissmann, M.; Blum, L. Trans. Faraday Soc. 1968, 64, 2605. (11) Ben Naim, A.; Stillinger, F. H. J. Chem. Phys. 1969, 51, 900. (12) Stillinger, F. H.; Rahman, A. J. Chem. Phys. 1974, 61, 4973. (13) Bopp, P.; Jancso´, G.; Heinziger, K. Chem. Phys. Lett. 1983, 98, 129. (14) Jorgensen, W. L. J. Am. Chem. Soc. 1978, 100, 7824; 1981, 103, 335; J. Chem. Phys. 1982, 77, 4156. (15) Berendsen, H. J. C.; Postma, J. P. M.; van Gunsteren, W. F.; Hermans, J. In Intermolecular Forces; Pullman, B., Ed.; Reidel: Dordrecht, 1981; p 331. (16) Berendsen, H. J. C.; Grigera, J. R.; Straatsma, T. P. J. Phys. Chem. 1987, 91, 6269.
(114)
1 0 0 T6 ) 0 0 0 and the inverse
-b1/3 - 3b3/5 3b3/4 -b1/4 - b3/2
( ) ( )
where
(111)
-Λ(03)
1 - 4η -3η -12/5η M BT ) -b1/2 - 3b3/4
K003 K033 2K103 2K133 3K203 3K233
It will be convenient to rearrange these matrices by a similarity transformation
1 - Λ0 -Λ(03)
[I0] - [Λ0] )
K000 -K003 2K100 n[K0] ) -2K103 3K200 -3K203
(110)
( )
1 0 0 I0 ) 0 0 0
( )
J. Phys. Chem., Vol. 100, No. 4, 1996 1205
1 + 2b(33) 2 /6 0 b(33) 2 /4
1
)
(17) Wertheim, M. S. J. Chem. Phys. 1971, 55, 4291. (18) Blum, L.; Fawcett, R. W. J. Phys. Chem. 1993, 96, 408. (19) Carnie, S. L.; Chan, D. Y. C.; Walker, G. R. Mol. Phys. 1981, 43, 1115. (20) Carnie, S. L.; Patey, G. N. Mol. Phys. 1982, 47, 1129. (21) Bratko, D.; Blum, L.; Luzar, A. J. Chem. Phys. 1985, 83, 6367. (22) Soper, A. K.; Phillips, M. G. Chem. Phys. 1986, 107, 47. (23) Blum, L.; Degreve, L. J. Phys. Chem., in press. (24) Baxter, R. J. J. Chem. Phys. 1968, 49, 2770. (25) Wertheim, M. S. J. Stat. Phys. 1984, 35, 19; 1984, 35, 42; 1986, 42, 459, 447; J. Chem. Phys. 1980, 73, 1398; 1981, 74, 2466; 1986, 85, 2929. (26) Baxter, R. J. J. Chem. Phys. 1968, 49, 2770. (27) Wertheim, M. S. J. Math. Phys. 1964, 5, 643. (28) Cummings, P. T.; Blum, L. J. Chem. Phys. 1986, 84, 1833. (29) Joshi, A. W. Matrices and Tensors in Physics; Wiley Eastern Limited: New Delhi, 1984. (30) Høye, J. S.; Blum, L. J. Stat. Phys. 1977, 16, 399. (31) Vericat, F.; Blum, L. Mol. Phys. 1982, 45, 1067. (32) Perram, J. W. Mol. Phys. 1975, 30, 1505. (33) Sceats, M. G.; Stavola, M.; Rice, S. A. J. Chem. Phys. 1979, 59, 2254. (34) Edmonds, A. R. Angular Momentum in Mechanics; Princeton University Press: Princeton, 1960.
JP950817S
