Dielectrically Nontrivial Closures for the RISM Integral Equation - The

11880

J. Phys. Chem. B 2001, 105, 11880-11892

Dielectrically Nontrivial Closures for the RISM Integral Equation† Fernando O. Raineri* and George Stell‡ Department of Chemistry, State UniVersity of New York at Stony Brook, Stony Brook, New York 11794-3400 ReceiVed: June 5, 2001; In Final Form: July 24, 2001

We report three new closures to the reference interaction site method (RISM) integral equation. The formulation of the closures is motivated by an analogy with an alternative form of the Ornstein-Zernike hypernettedchain closure for a simple ionic mixture. Unlike other RISM integral equation closures proposed to date, the closures reported here predict nontrivial results for the dielectric constant 0 of model polar fluids. Interestingly, we find that, for models of polar fluids comprising molecules with only two nonequivalent interaction sites (such as dipolar dumbbells or the popular three-site representations of water), all of the new closures give simple analytical formulas for 0. These simple analytical results indicate that only one of the closures is based on a relation that can account for the dependence of 0 upon the thermodynamic state parameters of the liquid and the features of the molecular model (size and shape) separately, rather than in the ubiquitous dipolar strength parameter combination 4πFµ2/3kBT. For this promising closure we solve numerically the RISM integral equation and calculate the site-site correlation functions gjl(r) and the dielectric properties for dipolar dumbbell liquids, for the SPC/E model of water, and for a four-site interaction site model of dimethyl sulfoxide.

1. Introduction Andersen,1

Since its introduction in 1972 by Chandler and the reference interaction site method (RISM) has played an important role in the study of the structure and thermodynamic properties of molecular liquids.2-5 It is specifically adapted to deal with interaction site model (ISM) representations of fluids, in which the potential energy of interaction between two molecules is the sum of site-site terms, including Coulombic interactions between partial charges qj located at the molecular sites. In the RISM formalism the structure of the liquid is described in terms of the site-site correlation functions hjl(r), where the subscripts j and l label the interaction sites in two different molecules. The basic relation is the RISM integral equation

h˜ (k) ) ω ˜ (k)‚c˜ (k)‚ω ˜ (k) + ω ˜ (k)‚c˜ (k)‚G‚h˜ (k)

(1.1)

where ω ˜ (k), h˜ (k), and c˜ (k) are square matrices with elements equal to the Fourier transforms ω ˜ jl(k), h˜ jl(k), and c˜ jl(k) of, respectively, the site-site intramolecular correlation function ωjl(r), the site-site (indirect) pair correlation function hjl(r), and the site-site RISM direct correlation function cjl(r). Finally, G is the square diagonal matrix with elements equal to the number density of each interaction site (considered as distinguishable). Equation 1.1, which relates the two unknown matrices h˜ (k)and c˜ (k), must be supplemented by a closure relation; i.e., an independent interrelationship between the sets of correlation functions {hjl(r)} and {cjl(r)} in order to have a closed set of equations. In the original approximation of Chandler and Andersen, which was concerned with hard core ISM models, the closure is specified by the equations hjl(r) ) -1 for r < djl and cjl(r) ) -βujl(r) for r > djl, where djl is the contact distance between interaction sites j and l, ujl(r) is the potential energy of †

Part of the special issue “Howard Reiss Festschrift”. * Corresponding author. E-mail: [email protected]. ‡ E-mail: [email protected].

interaction between a site of type j and another site of type l in different molecules when the distance between the sites is r, and β ) (kBT)-1 is the inverse of the temperature expressed in energy units. A subsequent site-site hypernetted-chain (HNC) type closure, relevant to ISM models of more general form, was proposed by Hirata, Pettitt, and Rossky.6 Explicitly, the RISM-HNC [or extended-RISM (XRISM)] closure can be written as

hjl(r) ) exp[-βujl(r) + tjl(r)] - 1

(1.2)

tjl(r) ) hjl(r) - cjl(r)

(1.3)

For the molecular models of interest here, ujl(r) has the general form

ujl(r) ) u*jl(r) +

qjql r

(1.4)

where the non-Coulombic part u*jl(r) is typically taken to be of the Lennard-Jones type. It is important to notice that all of the functions in eqs 1.2 and 1.3 explicitly involve only the j and l interaction sites. Thus the RISM-HNC closure is a strictly local relation, involving only a single pair of sites. A characteristic feature of the HNC and most other closures used to date in conjunction with the RISM integral equation is the implication that, at large distance r, the RISM direct correlation functions cjl(r) are proportional to the corresponding potential energies of interaction between the interaction sites (just the Coulombic tail qjql/r at such large r)

cjl(r) ) -βqjql/r ≡ φjl(r) at large r

(1.5)

Intuitively, this large-r behavior would seem quite natural, since7 there are formally exact arguments which demonstrate that an analogous relation holds between the molecule-molecule direct correlation function c(1,2) and the molecular pair potential u(1,2)

10.1021/jp0121163 CCC: $20.00 © 2001 American Chemical Society Published on Web 09/27/2001

Closures for the RISM Integral Equation

J. Phys. Chem. B, Vol. 105, No. 47, 2001 11881

(the arguments 1 and 2 specify both location and orientation) at large intermolecular distance r

c(1,2) ) -βu(1,2) at large r

(1.6)

as long as one is not at a critical point.7 Unfortunately,8 sitesite closures that satisfy eq 1.5 give a “trivial” result11,12 for the dielectric constant 0 of the liquid, namely the ideal gas value

0 ) 1 + V

(1.7)

here expressed in terms of the dipolar strength parameter V of the liquid

V ≡ 3y ) 4πβµ2F/3

(1.8)

In this equation µ is the strength of the molecular dipole moment (the molecules we treat here are assumed to be nonpolarizable) while F is the number density of the solvent. This shortcoming of the HNC and other site-site closures prompted several efforts to correct or ameliorate the problem. These modifications are based on a result by Cummings and Stell13,14 for the special case of liquids comprising polar molecules with at most two inequiValent interaction sites. Cummings and Stell showed that for such systems a more complicated relation

cjl(r) ) A(V,0) φjl(r) at large r

(1.9)

holds between the RISM site-site cjl(r) and the corresponding site-site potential at large r. The scalar factor A(V,0) is determined by the dielectric constant 0 of the liquid and by the dipolar strength parameter V introduced in eq 1.8:

A(V,0) )

1 + 0(V - 1)

(1.10)

V(1 - 0)

Stimulated by the simplicity of eqs 1.9 and 1.10, Rossky, Pettitt, and Stell15 proposed slightly modifying the site-site HNC closure, eq 1.2, to be asymptotically compatible with eq 1.9 by letting

hjl(r) ) exp[-βu′jl(r;V,in 0 ) + tjl(r)] - 1

(1.11)

In this paper we investigate a new family of HNC-like closures to the RISM integral equation that give highly nontrivial and quantitatively useful results for the dielectric constant of ISM model fluids. The proposed closures preserve the form of eq 1.2, namely

hjl(r) ) exp[-βujl(r) + τjl(r)] - 1

and are distinguished by the particular form of τjl(r) that replaces the function tjl(r), eq 1.3, of the HNC closure. The outline of the paper is as follows. In section 2 we briefly recall the dielectric formulas for interaction site models. Section 3 presents the main ideas on which we base the particular choices for the τjl(r) function for the first two closures. In section 4 we implement the dielectric formulas of section 2 under these new closures and find that 0 can be calculated analytically for a particular class of ISM models. Although the results are encouraging, we point out that for this particular class of models closures I and II share with RISM-HNC an important limitation, namely that 0 is predicted to depend on the fluid thermodynamic state and the molecular parameters only through the dipolar strength parameter V (cf. eq 1.8). This feature motivates our departure in section 5 from our original strategy. There we consider an alternative closure relation, closure III, that results from the combination of the HNC tjl(r) and some projected terms extracted from the closure I τjl(r) function. We then show, based on numerical examples, that closure III has the potential for predicting realistic estimates for the dielectric constant of molecular fluids while still giving good results for the structure of the liquid [as represented by the site-site correlation functions gjl(r)]. We conclude the paper in section 6 with a brief summary of the results obtained. Several equations and techniques used in the text are discussed in the appendices. 2. Dielectric Formulas for ISM Models In this section we present the dielectric formulas, based on site-site correlation functions, from which we can calculate the dielectric constant 0 of a model fluid comprising rigid nonpolarizable molecules with p interaction sites. For convenience we adopt the matrix notation of eq 1.1. 1. The first dielectric formula11,12,18 allows the calculation of 0 from the site-site correlation functions hjl(r):

1 - 1 ) -4πβ(q|S2|q) 0

in which the original site-site potential energy of interaction ujl(r), cf. eq 1.4, is replaced by a modified site-site potential

qjql r

/ in u′jl(r;V,in 0 ) ) ujl(r) + A(V,0 )

(1.13)

(2.1)

where S2(k) is the coefficient of k2 in the small-k expansion

(1.12)

in which the long-range Coulombic interaction has been rescaled. The parameter A(V,in 0 )is selected as in eq 1.10, with the argument 0 specialized at the targeted value in 0 of the dielectric constant of the fluid.15 The modified closure, eqs 1.11 and 1.3, is then guaranteed to produce site-site correlation functions hjl(r) that are consistent with the value of in 0 supplied as input to the calculation. This methodology is clearly not predictive with respect to 0; nevertheless it provides a robust method to calculate the structure of polar liquids and at the same time allows the detailed study of their wave vector dependent dielectric properties.16,17 It is important to remember, however, that this modification of the HNC closure is inspired by the Stell-Cummings result, eq 1.9, which only applies to liquids of polar ISM molecules with at most two inequivalent interaction sites.

S(k) ) k0S0 + k2S2 + ...

(2.2)

S(k) ) ω ˜ (k)‚G + G‚h˜ (k)‚G

(2.3)

of the matrix

with elements the site-site structure factors Sjl(k). In eq 2.1 |q) ) col(q1,q2,...,qp) is a column vector of dimension p whose elements are the partial charges at the molecular sites of the liquid molecules, while (q| is the transpose row vector. The matrix functions h˜ (k) and ω ˜ (k) in eq 2.3 have the small-k expansions

h˜ (k) ) k0h0 + k2h2 + ...

(2.4)

ω ˜ (k) ) k0ω0 + k2ω2 + ...

(2.5)

11882 J. Phys. Chem. B, Vol. 105, No. 47, 2001

Raineri and Stell

where the matrix coefficients ω0 and h0 have the simple dyadic form14

3. From the site-site RISM direct correlation functions cjl(r):

h0 ) h0|1)(1|

(2.6)

ω0 ) |1)(1|

1 - 1 ) -4πβ(q|{(ω2G)-1 - c-2}-1|q) 0

(2.7)

Here |1) is a column vector of dimension p with all its elements equal to unity. The electroneutrality of the molecules therefore implies the identities ω0|q) ) h0|q) ) |0) and (q|ω0 ) (q|h0 ) (0|. From eq 2.3 and the small-k expansions of ω0 and h0,we immediately conclude that S2 in eq 2.1 is given by the expression

S2 ) ω2‚G + G‚h2‚G

(2.8)

(2.9)

where S*2 is the coefficient of k2 in the small-k expansion

S*(k) ) k0S*0 + k2S*2 + ...

(2.10)

of the matrix of modified site-site structure factors (or sitesite hyperVertex functions)

S*(k) ) ω ˜ (k)‚G + G‚h˜ *(k)‚G

(2.11)

The modified h˜ /jl(k) in eq 2.11 are defined from the RISM-like integral equation (compare with eq 1.1)

h˜ *(k) )ω ˜ (k)‚c˜ *(k)‚ω ˜ (k) + ω ˜ (k)‚c˜ *(k)‚G‚h˜ *(k)

(2.12)

in terms of the modified RISM direct correlation functions

c˜ *(k) ) c˜ (k) - O˜(k)

(2.13)

The elements of the square matrix

O˜(k) ) (-4πβ/k2) |q)(q|

(2.14)

are the Fourier transforms of the functions φjl(r) defined in eq 1.5. The dyadic notation in eq 2.14 emphasizes that the rank of matrix φ˜ (k) is unity. The small-k expansion of h˜ *(k) is

h˜ *(k) ) k0h*0 + k2h*2 + O(k4)

(2.15)

with h* 0 ) h0. It then follows that S* 2 in eq 2.9 is given by

S*2 ) ω2‚G + G‚h* 2‚G

(2.16)

In view of the RISM equations (1.1) and (2.12), we can also express S2 and S* 2 in terms of the coefficients c-2 and c* -2 in the small-k expansions

Y ˜ (k) ) k-2Y-2 + k0Y0 + ...

(2.17)

where Y ) c or c*. Using eqs 2.3 and 2.11, we calculate

X2 ) [(ω2‚G)-1 - Y-2]-1

4. From the site-site modified RISM direct correlation functions c/jl(r):

0 - 1 ) 4πβ(q|{(ω2G)-1 - c*-2}-1|q)

(2.20)

From the last equation it is clear that 0 reduces to the ideal gas value when c*-2 ) 0. As mentioned in the Introduction, this behavior is typical of any RISM closure for which eq 1.5 applies.11,12 3. Development of Nontrivial Dielectric Closures

2. An alternative expression for 0 is19

0 - 1 ) -4πβ (q|S*2|q)

(2.19)

(2.18)

where the pair (X2,Y-2) is (S2,c-2) or (S* 2,c* -2). In this way we generate the alternative dielectric formulas:

As mentioned in the Introduction, our choices of τjl(r) functions in eq 1.13 is inspired in the unusual representation

tjl(r) ) hjl(r) - cjl(r) ) F -1{c˜ (k)‚S(k)‚c˜ (k)}jl (SIM) (3.1) for tjl(r) in a simple ionic mixture (SIM); i.e., a multicomponent liquid system of monatomic ions, such as mixture of molten salts. The second equality in eq 3.1 may be easily proved with the help of the Ornstein-Zernike equation for a SIM. In this equation the operator F -1 indicates the inverse Fourier transform in three dimensions. The square matrices c˜ (k) and S(k) have elements, respectively, c˜ jl(k) [the Fourier transform of the direct correlation function cjl(r)] and

Sjl(k) ) Fjδjl + FjFl h˜ jl(k) (SIM)

(3.2)

the partial structure factor between ionic species j and l in the SIM. Finally, Fj is the number density of ionic species j, while δjl is the Kro¨necker symbol. 3.1. New RISM Closure I. As the first alternative closure for the RISM integral equation, we propose to adopt the HNClike eq 1.13 together with

τjl(r) ) F -1{c˜ (k)‚S(k)‚c˜ (k)}jl (I)

(3.3)

This choice for τjl(r) is clearly fashioned after eq 3.1 for the SIM, but of course the symbols in eq 3.3 are interpreted differently. In particular, the subscripts j and l now refer to the interaction site types of polyatomic molecules, while c˜ (k) and S(k) represent, respectively, the matrices of site-site RISM direct correlation functions c˜jl(k) and the matrix of site-site structure factors (cf. eq 2.3). We now express the matrix τ(k) [whose elements are the Fourier transformed functions τjl(r)] in a form that can be directly contrasted with the Fourier transformed version ˜t(k) ) h˜ (k) - c˜ (k) of eq 1.3. Using eq 1.1, we derive

τ˜ (k) ) H ˜ (k) - c˜ (k) (I)

(3.4)

where H ˜ (k) is related to the matrix h˜ (k) of site-site correlation functions by the equation

H ˜ (k) ) Ω ˜ (k)‚h˜ (k)‚Ω ˜ (k) (I)

(3.5)

Matrix Ω ˜ (k) is, in turn, defined as the matrix inverse of the intramolecular correlation function matrix

Ω ˜ (k) ≡ [ω ˜ (k)]-1

(3.6)



The matrix product structure of H ˜ (k) gives rise to an important difference between τjl(r) of this integral equation closure20 and tjl(r) of the RISM-HNC closure, namely that τjl(r) is not only generated from hjl(r) and cjl(r), but also from the site-site correlation functions hmn(r) associated with other pairs mn of interaction sites. Thus, closure I is a nonlocal relation in r-space.20 The difference between H ˜ (k) and h˜ (k) [and, correspondingly, between τ˜ (k) and ˜t(k)] at small values of k is dramatic. In Appendix A we show that the matrix Ω ˜ (k) has the following behavior at small values of the wave vector k

Ω ˜ (k) ) k-2Ω-2 + k0Ω0 + O(k2)

(3.7)

and give explicit expressions for the matrix coefficient Ω-2. Equations 2.4 and 3.7 for, respectively, h˜ (k) and Ω ˜ (k), imply τ˜ (k) ) k-2τ-2 + k0τ0 + O(k2) when k is small. It then follows from eq 1.13 that the Fourier transform of the RISM direct correlation function matrix has the following small-k behavior

(3.8)

Of immediate significance for the calculation of the dielectric constant 0 is the expression of the c-2 matrix coefficient:

c-2 ) -4πβ|q)(q| + H-2 (I)

(3.9)

H-2 ) Ω-2‚h2‚Ω-2 (I)

(3.10)

where

The H-2 term is responsible for this new RISM closure I to give 0 different from the ideal gas value, eq 1.7. 3.2. New RISM Closure II. Based on the same idea guiding the construction of closure I, we now propose a second closure by separately dealing with contributions τjl*(r) and τjlQ(r) to τjl(r) that originate from the partition of the RISM direct correlation function as in eq 2.13. The new closure II corresponds to eq 1.13 with the following interpretation of the function τjl(r):

τjl(r) ) τ*jl(r) + τQjl (r) (II)

(3.11)

τ*jl(r) ) F -1{c˜ *(k)‚S*(k)‚c˜ *(k)}jl

(3.12)

)F

-1

{O˜(k)‚S (k)‚O˜(k)}jl Q

(3.13)

Here, in addition to S*(k) defined in eq 2.11, we have introduced the auxiliary structure factor

˜ (k)‚G + G‚Q ˜ (k)‚G ) [(ω ˜ (k)‚G) S (k) ) ω Q

-1

-1

- O˜(k)]

(3.14)

in which the square matrix Q ˜ (k) is implicitly defined by the RISM-like equation

Q ˜ (k) ) ω ˜ (k)‚O˜(k)‚ω ˜ (k) + ω ˜ (k)‚O˜(k)‚G‚Q ˜ (k) (3.15) in terms of matrix O˜(k) introduced in eq 2.14. At small values of the wave vector k matrix Q ˜ (k) has the behavior

Q ˜ (k) ) k2Q2 + O(k4)

(3.17)

H ˜ ′(k) ) Ω ˜ (k)‚[h˜ *(k) + Q ˜ (k)]‚Ω ˜ (k) (II)

(3.18)

and

Here it is important to note that the sum h˜ *(k) + Q ˜ (k) does not equal the full matrix of site-site correlation functions h˜ (k), so that H ˜ (k) * H ˜ ′(k). Using eqs 2.15 and 3.16, and following steps similar to those at the end of section 3.1, it is straightforward to show that under the new closure II the RISM site-site direct correlation function matrix c˜ (k) again displays the small-k behavior of eq 3.8 with coefficient matrix c-2 given by

c-2 ) -4πβ|q)(q| + H′-2 (II)

(3.19)

and H′-2 given by

H′-2 ) Ω-2‚[h2* + Q2]‚Ω-2 (II)

(3.20)

4. Dielectric Constant under Closures I and II

c˜ (k) ) k-2c-2 + k0c0 + O(k2)

τQjl (r)

τ˜ (k) ) H ˜ ′(k) - c˜ (k) (II)

(3.16)

Proceeding as in section 3.1, it is simple to show that matrix τ˜ (k) for the new RISM closure II is given by the pair of equations

It is noteworthy that just by using the functional form that the HNC closure assumes for a SIM mixture, eq 3.1, we were able to develop two new closures for the RISM integral equation that give (for any ISM model) highly nontrivial results for the coefficient matrix c* -2. Inspection of eqs 3.10 and 3.20, together with the fact that c*-2 ) H-2 and c*-2 ) H′-2 for, respectively, closures I and II, leads to the conclusion that c* -2 is not zero for these closures because of the role of the inverse intramolecular correlation function matrix Ω ˜ (k). We now consider the calculation of the dielectric constant 0 of polar nonpolarizable fluids under closures I and II. We begin by noticing that the two h-based dielectric fluctuation formulas, eqs 2.1 and 2.9, can be written more concisely as follows:

1 - 1 ) -V - B 0

(4.1)

0 - 1 ) V + B*

(4.2)

The first term on the right-hand side of these equations, which involves the dipolar strength parameter V (cf. eq 1.8), originates from the intramolecular contribution to the matrix coefficients S2 and S*2 (eqs 2.8 and 2.16, respectively)

V ) 4πβ(q|ω2‚G|q) ) 4πβF(q|ω2|q)

(4.3)

The B-terms in eqs 4.1 and 4.2 originate, on the other hand, from the intermolecular contributions to S2 and S* 2

B ≡ 4πβ(q|G‚h2‚G|q) ) 4πβF2(q|h2|q)

(4.4)

B* ≡ 4πβ(q|G‚h*2‚G|q) ) 4πβF2(q|h*2|q)

(4.5)

In general, the calculation of 0 with either eq 4.1 or 4.2 requires explicit knowledge of the full set of site-site correlation functions hjl(r). There is a particular class of ISM models, however, for which we can derive simple analytical formulas for 0 under the new closures I and II. This is the class of twosite-polar ISM models, namely multisite ISM models comprising only two inequiValent types of interaction sites. For these models all of the interaction sites of one type have a partial charge of positive sign, while all of the remaining interaction sites of the other type have a partial charge of negative sign.21 Familiar examples of ISM models that belong to the two-site-polar ISM class are dipolar dumbbells and the popular three-site-non-


Raineri and Stell

polarizable models of water (SPC, SPC/E, TIPS3). As we shall see, the determination of 0 for this class of models does not demand the explicit calculation of the intermolecular factors B and B*. 4.1. Dielectric Constant for ISM Models in the Two-SitePolar Class. What actually makes the analytical determination of 0 feasible for this class of models is the dyadic form of the coefficient matrix Ω-2. For the simplest example of a model liquid that pertains to this class, consider the case of a liquid of dipolar dumbbells, i.e., a liquid of diatomic molecules of bond length l and partial charges q and -q at the atomic sites. When we apply eq A.7 to this model to calculate the coefficient matrix Ω-2, we get

Ω-2 )

(

3 1 l2 1

1 -1

)

(4.6)

A more convenient but entirely equivalent representation of Ω-2 for this model is

Ω-2 ) (3/µ2)|q)(q|

(4.7)

where we have expressed the rank 1 matrix in terms of the dyadic constructed with the column vector |q), of elements q and -q. Furthermore, µ ≡ ql is the magnitude of the dipole moment of the diatomic molecule. Although specifically derived for a dipolar dumbbell model, the 2 × 2 matrix Ω-2 in the form of eq 4.7 also applies to any model that belongs to the two-site-polar ISM class when the symmetry-reduced form of the RISM integral equation is used.21 We therefore exploit this very simple expression for Ω-2 to derive the analytical formulas for 0 under the two new RISM integral equation closures. Equations 2.13 and 3.9 indicate that c*-2 ) H-2 applies under closure I. For the new closure II, on the other hand, eqs 2.13 and 3.19 imply c* -2 ) H′-2. The coefficient matrices H-2 and H′-2 were defined, respectively, by eqs 3.10 and 3.20. When we take into account eq 4.7, we find that c*-2 for models in the two-site class under closures I and II are given by

c*-2 ) 4πβ(B/V2)|q)(q| (I)

(4.8)

c*-2 ) 4πβ[(B* + Q2)/V2]|q)(q| (II)

(4.9)

where Q2 ≡ 2|q). In Appendix B we show that Q2 ) -V2/(1 + V). It is interesting to point out that both new closures predict that for models in the two-site-polar ISM class matrix c*-2 is proportional to the coefficient matrix O-2, in agreement with Cummings and Stell13,14 [cf. eq 1.9, which implies c*-2 ) (A - 1)O-2]. Using eqs 4.8 and 4.9 in the dielectric formula eq 2.20, we obtain 4πβF2(q|Q

(

0 - 1 ) V/ 1 -

(

0 - 1 ) V/ 1 -

B (I) V

)

(4.10)

B* V (II) + V 1+V

)

(4.11)

With these results we can now obtain closed formulas for the dielectric constant of liquids comprising molecules represented by ISM models that belong to the two-site-polar class. Closure I. Through the elimination of B from eqs 4.1 and 4.10, we derive the following analytical expression for 0 as a

Figure 1. Dependence of the dielectric constant 0 with the dipolar strength parameter V for model fluids with only two inequivalent interaction sites. HNC: 0 calculated with eq 1.7; Clos. I, 0 calculated with eq 4.12; Clos. II, 0 calculated with eq 4.14.

function of the dipolar strength parameter V:

0(V) )

1 V2 - 2 + 2V + V xV2 + 4V (I) 2 2V - 1

(4.12)

At first sight it appears that 0 has a pole at V ) 1/2, but this is only apparent, because the numerator also has a simple root at the same value of V. Equation 4.12 implies the following dependence of 0 with V at small values of the latter:

3 17 0(V) ) 1 + V - V3/2 + V2 - V5/2 + ... (I) (4.13) 2 2 Closure II. For the new closure II, after eliminating B* from eqs 4.2 and 4.11, we calculate

0(V) ) 1 + V +

V V + x5V2 + 4V (II) 2 1+V

(4.14)

In this case there are no zeros in the denominator of the third term for the physically relevant values of the dipolar strength parameter V. The dependence of 0 with the dipolar strength parameter is

1 1 3 0(V) ) 1 + V + V3/2 + V2 - V5/2 - ... (II) 2 4 16

(4.15)

Equations 4.12 and 4.14 apply to any polar liquid whose molecules belong to the two-site-polar class of ISM models. According to these equations, under closures I and II, 0 is exclusively determined by the single parameter V. This prediction, however, is at odds with calculations of 0 for liquids of symmetric dipolar dumbbells by computer simulation,23,24 which reveal that two such liquids with molecules with dipole moments of the same magnitude µ but different bond lengths have at temperature T and density F different dielectric constants. In Figure 1 we represent 0 (calculated with eqs 4.12 and 4.14) as a function of the dipolar strength parameter V. For reference we also include in the figure the line 0 ) 1 + V that corresponds to the result under the HNC closure. We observe that 0 predicted by closure I is always smaller than the ideal gas reference behavior (i.e., 0 under RISM-HNC), which again is at odds with simulation results for liquids that belong to the two-site-polar ISM class. In contrast, 0 predicted by closure II is larger than the ideal gas reference behavior, and this difference increases with increasing V.



To put the results for 0 by the different closures in perspective, we consider the specific example of the extended simple point charge (SPC/E) model of water at room temperature.25 The SPC/E model has three interaction sites, but due to the equivalence of the hydrogen sites, the model belongs to the two-site-polar ISM class. At 298 K the dipolar strength parameter of the liquid is V ) 18.8. The values of 0 predicted by RISM-HNC (eq 1.7), closure I (eq 4.12), and closure II (eq 4.14) are, respectively, 19.8, 10.6, and 49.1. The dielectric constant calculated by molecular dynamics simulation for the model is 70.26 Although still quite far from the simulation result, 0 predicted by closure II is certainly the best prediction among the three closures. The remaining question, of course, concerns the structure functions predicted by the two new integral equation closures. We have attempted solving both integral equations for fluids of symmetric dipolar dumbbells, but encountered in the process some technical difficulties. Rather than delving into these difficulties here, we turn instead to the development of another new integral equation closure that attempts to overcome the shortcomings of closures I and II discussed in this section.

τPjl (r) is a variant of the τjl(r) function, eqs 3.4 and 3.5, of our new closure I:

τPjl (r) ) H′′jl(r) - cjl(r) where H′′jl(r) ) F -1{H ˜ ′′(k)}jl and

1 H ˜ ′′(k) ) [P‚H ˜ (k) + H ˜ (k)‚P] 2

(5.4)

Here H ˜ (k) is the matrix (3.5) of closure I, while P is the idempotent matrix (projection operator)

P ) |q)(q|q)-1(q|

(5.5)

We use the symmetric form (1/2)[P‚H ˜ +H ˜ ‚P] rather than P‚ H ˜ ‚P so that the rank of the resulting coefficient matrix c*-2 will be greater than 1 for ISM models with three or more types of inequivalent interaction sites. With the help of eqs 1.3 and 5.4 we can write the hybrid τθjl(r) function in the more convenient form

τθjl(r) ) tjl(r) + θ[H′′jl(r) - hjl(r)] (III)

5. Closure III 5.1. General Formulation. In view of the shortcomings of closures I and II, we now propose a third closure that may be viewed as a hybrid between the RISM-HNC closure and a modified version of the new closure I. This closure is designed with the following guidelines in mind: 1. We require that for models with three or more inequivalent interaction sites (i.e., models of polar liquids that do not belong to the two-site-polar ISM class) the rank of the coefficient matrix c* -2 is greater than 1. That this must be so is easily seen from the general expression

c*-2 ) Ω-2‚h*2‚Ω-2‚[I + G‚h*2‚Ω -2]-1

(5.1)

that follows from the RISM-like equation (2.12) together with eqs 2.15, 2.6, and A.8, and the fact that h*0 ) h0. Only for models that belong to the two-site-polar ISM class the rank of c* -2 is 1. 2. It is reasonable to expect that the modifications to RISMHNC required to generate a nonvanishing coefficient matrix c* -2 should affect mainly the large-r behavior of the correlation functions hjl(r). The rationale for this expectation is the robust performance of the RISM-HNC closure to predict the shortrange structure in liquids and solutions. Therefore, in crafting an improved integral equation closure, we choose to modify the RISM-HNC closure mostly at large values of r, which is the regime that matters for the determination of the dielectric constant. 3. Finally, we need a closure that incorporates naturally the subtle sensitivity of the dielectric constant with the size and shape of the molecular model, such as the specific dependence23,24 of 0 with the bond length for symmetric dipolar dumbbell model liquids. Based on these general guidelines, we now propose a new integral equation closure similar to eq 1.13, but for the τjl(r) function we choose the hybrid form

τθjl(r)

(5.3)

) (1 - θ)tjl(r) +

θτPjl (r)

(III)

(5.2)

Here tjl(r) is the RISM-HNC function given in eq 1.3, while

(5.6)

With the results of section 3.1 we can show that c˜ (k) predicted by the new closure follows at small k the behavior displayed in eq 3.8, with the following expression for the c-2 matrix coefficient:

c-2 ) -4πβ|q)(q| + θH′′-2 (III)

(5.7)

From eq 5.4 we have

1 H′′-2 ) [P‚H-2 + H-2‚P] (III) 2

(5.8)

The reader should notice the explicit appearance of the mixing parameter θ in the expression of c*-2 (the second term on the right-hand side of the last equation). 5.2. Dielectric Constant for two-site ISM Models. As was the case with closures I and II, we can obtain a simple analytical formula for 0 under closure III for liquids comprising molecules represented by ISM models that belong to the two-site-polar ISM class. Equations 5.5 and 4.7 for, respectively, P and Ω-2 imply the simplifications P‚Ω-2 ) Ω-2‚P ) Ω-2. Equations 2.13, 3.10, 5.7, and 5.8 then give

c*-2 ) θH-2

(5.9)

which differs from the result for closure I only by the presence of the weight parameter θ. At this point we can repeat step by step the derivation in section 4.1 to get

(

B V

0 - 1 ) V/ 1 - θ

)

(5.10)

Combining this result with eq 4.1 to eliminate the intermolecular factor B, we finally derive

[

0(θ,V) ) 1 - V -

]

(V - θ + 1 - xV2 + 2Vθ + 2V + θ2 - 2θ + 1)V 2θ

-1

(5.11)

which explicitly exhibits the dependence of the dielectric


Raineri and Stell

Figure 2. Dependence of the dielectric constant 0 for closure III with the weight parameter θ at several values of the dipolar strength parameter V for model fluids with only two inequivalent interaction sites. The curves have been calculated with eq 5.11.

Figure 4. Oxygen-oxygen site-site radial distribution function for SPC/E water at 298 K. The solid line corresponds to gOO(r) calculated with the RISM-closure III. The dashed line corresponds to gOO(r) calculated with RISM-HNC. r is in units of Å.

Figure 3. Dielectric constant 0 according to closure III for SPC/E water at 298 K as a function of the mixing parameter θ. The curve is calculated with eq 5.11 using V ) 18.8.

room temperature, and how the calculated structure of the liquid compares with the structural predictions of RISM-HNC. To answer these important questions, we have solved the new integral equation numerically for the SPC/E model water at 298 K. The site-site correlation functions hjl(r) and the dielectric constant of the model liquid were determined for a range of values of the mixing parameter θ, starting from θ ) 0 (which corresponds to RISM-HNC) to as low as θ ) -0.8. Technical details of the numerical methodology are briefly described in Appendix C. In Figure 4 we report results for the oxygenoxygen site-site correlation function gOO(r) as a function of the separation r. The figure compares gOO(r) calculated under RISM-HNC with the same correlation function calculated under RISM-closure III with θ ) -0.77, for which 0 ) 70. It is apparent from the figure that only relatively minor changes are brought on gOO(r) by the addition of the term θ[H′′jl(r) - hjl(r)] to the HNC tjl(r) function (cf. eq 5.6) in the closure relation. This feature, we recall, was one of the criterions that guided the development of function closure III (see item 2 just above eq 5.4). The position of the first peak of gOO(r) is the same for both closures, and there are only very minor differences in the height and width of the peak. There are, however, small qualitative differences in the predictions of both closures for r beyond the first peak, especially the development of a shoulder at r = 4.2 Å in the case of closure III. The location of this peak is close to the position of the second maximun of gOO(r) calculated by computer simulation. In Figure 5 we report the wave vector-dependent longitudinal dielectric function L(k) for SPC/E water at 298 K calculated both under HNC and under the new closure III with θ ) -0.77. The dielectric function is calculated with the equation16

constant 0 of a two-site-polar class liquid on the dipolar strength parameter V (cf. eq 4.3), and also on the weight parameter θ. The dependence of 0 with θ is quite interesting, as it opens the possibility for closure III to capture the separate influence of the molecular size and shape on the dielectric constant of a polar liquid. To visualize the predictions of eq 5.11 in Figure 2 we represent 0 as a function of θ for three values of the dipolar strength parameter V. Clearly, for any θ the function 0(θ,V) increases as V increases. For V fixed, on the other hand, 0(θ,V) is a decreasing function of θ. In particular 0(θ)1,V), which coincides with 0 under new closure I (cf. eq 4.12), is considerably smaller than 0(θ)0,V), the value of the dielectric constant predicted by the HNC closure (cf. eq 1.7). This feature illustrates the fact, already clear in Figure 1, that for model liquids belonging to the two-site-polar class 0 under RISMHNC is always larger than 0 calculated under the new closure I. Of considerably more interest, however, is the pronounced sensitivity of 0 with θ as the latter becomes more negative. In Figure 3 we report the dependence of 0 with θ (calculated with eq 5.11) when the value of the dipolar strength parameter is fixed at V ) 18.8, which is appropriate for the SPC/E model of water at 298 K. The figure shows, intriguingly, that 0 predicted by closure III when θ = -0.77 is approximately 70, which is the value of the dielectric constant of SPC/E water calculated by computer simulation.26 Naturally, this raises the question of whether we can actually solve the RISM-closure III integral equation in the range -0.8 e θ e -0.7 for SPC/E water at

1 -4πβ (q|S(k)|q) -1) L(k) k2

(5.12)

where the combination (q|S(k)|q) is the polarization chargecharge structure factor of the liquid. The figure shows that at large k (small distance) 1/L(k) is the same for both integral equation approximations. This is not surprising, as the features of the longitudinal dielectric function at large k are dominated ˜ (k)|q) on the by the intramolecular contribution (-4πβF/k2)(q|ω right-hand side of eq 5.12. On the scale of the figure, the difference between {1/L(k)}k)0 ) 1/0 for the two approximations (1/19.8 for RISM-HNC and 1/70 for closure III) cannot be appreciated. The discrepancy between the theories is most manifest around k ) km, the location of the minimum of 1/L(k). Interestingly, the value of km according to both theories



Figure 5. Longitudinal dielectric function, eq 5.12, for SPC/E water at 298 K. The solid line corresponds to L(k) calculated with the RISMclosure III. The dashed line corresponds to L(k) calculated with RISMHNC. k is in units of Å-1.

is the same. We refer the reader to ref 16 for a detailed discussion of the meaning and the implications of L(k) being negative at intermediate values of k. 5.3. General Multisite Solvents under Closure III. We have applied the numerical method used in the last section (see Appendix C) to calculate the structure functions and the dielectric constant for several polar liquids with molecules with more than two inequivalent types of interaction sites (e.g., methanol, dimethyl sulfoxide, dichloromethane). In each case we encountered no problems in solving the RISM-closure III integral equation with a value of θ for which the correlation functions and the dielectric formula (2.1) give 0 equal to the simulated dielectric constant for the liquid (or, for dichloromethane, the experimental 0). Encouragingly, the calculated site-site correlation functions are very similar to the gjl(r) calculated with RISM-HNC, and are also reasonably consistent with the liquid structure generated by computer simulation methods. The coefficient matrix Ω-2 of ISM molecules with more than two inequivalent interaction sites is not related in any way to the Coulomb coefficient matrix O-2. For such liquids, then, closure III says that (cf. eq 5.7)

c*-2 ) θH′′-2

(5.13)

where the rank of matrix H′′-2 (defined by eq 5.8) is greater than 1. The consequence of this is that relation 1.9 does not hold for polar liquids that do not belong to the two-site-polar class under closure III.27 More specifically, the analytical formula, eq 5.11, is not valid for models outside the two-sitepolar ISM class; the determination of 0 for these models requires the calculation of the intermolecular factor B (cf. eq 4.1), which, in turn, demands detailed knowledge of the sitesite correlation functions. Nevertherless, it is interesting to ask how much 0 as a function of the mixing parameter θ departs, for an ISM polar liquid with more than two inequivalent interaction sites, from 0(θ,V) given by eq 5.11. Figure 6 illustrates such a comparison for dimethyl sulfoxide (DMSO) at 298 K. For the calculations in the figure the DMSO molecules are represented with the foursite ISM model of Liu et al.28 (each of the methyl groups is represented by one interaction site). Represented in Figure 6 is the calculated dependence of 0 with θ as calculated (i) with dielectric formula (2.1) and RISMclosure III site-site correlation functions hjl(r) (the upper curve) and (ii) with eq 5.11 (the lower curve). In both cases the dipolar

Figure 6. Dielectric constant 0 according to closure III for dimethyl sulfoxide at 298 K as a function of the mixing parameter θ. The upper curve is calculated with the dielectric formula (2.1) and the site-site correlation functions hjl(r) obtained from the numerical solution of the RISM-closure III integral equation. The lower curve is generated with eq 5.11 with the dipolar strength parameter V ) 18.2 appropriate for the model.

strength parameter is V ) 18.2. The figure shows that for this DMSO model both estimates of 0 are very similar when θ is in the range -0.25 e θ e 0.0. However, for θ < -0.25 the difference between these two estimates of 0 grows at a considerable rate as θ becomes more and more negative. For example, when θ ) -0.475 the dielectric constant according to the numerical solution of RISM-closure III is 0 ) 43.4, which is close to the value 0 ) 44 determined for the Liu et al. model of DMSO from a molecular dynamics computer simulation.29 At the same value of θ, the two-site-polar calculation of 0 with eq 5.11 gives 0 ) 34. The implication of Figure 6 is that the lack of a simple connection between c*-2 and O-2 for ISM models with more than two inequivalent interaction sites (a property of the RISM direct correlation function matrix c˜ *(k) respected by closure III) is quite important for the determination of the dielectric constant. 5.4. A Priori Determination of θ. So far we have applied closure III by presetting the value of θ such that the site-site correlation functions hjl(r) produce, when replaced in one of the dielectric formulas, a desired result of the dielectric constant (for example, the value of 0 calculated for the liquid by computer simulation). Of course, in order for the closure to be useful as a predictive route for the dielectric constant of a liquid, we must also specify a methodology for fixing θ without recourse to external input. A first approach to this problem is to exploit the analytical formula for 0, eq 5.11, valid only for ISM models that belong to the two-site-polar class. We make the reasonable assumption that θ is an analytic function of the dipolar strength parameter V of the liquid, and that for small values of V it can be approximated by the Taylor expansion

θ(V) ) θ0 + θ1V + θ2V2 + ...

(5.14)

The coefficients θ0, θ1, etc. may depend separately on the thermodynamic state parameters T and F, as well as on the particular features (size, shape, and charge distribution) of the molecules. Equation 5.11 for 0(θ,V) and eq 5.14 for θ imply

0(V) ) 1 + V + R2V2 + R3V3 + R4V4 + ...

(5.15)

at small values of V. Assuming that θ0 < 0 (as suggested by the results of section 5.2), we derive the following expressions


Raineri and Stell

for the coefficients Rn, n ) 2, 3, 4:

R2 ) R3 ) R4 )

θ0 θ0 - 1

θ03 - θ02 - θ0θ1 - θ0 + θ1 (θ0 - 1)3

(5.16)

(5.17)

θ05 - 2θ04 - g3θ03 + g2θ02 + g1θ0 - g0 (θ0 - 1)5

(5.18)

liquid is d ) 3.876 Å and the reduced density is F* ) Fd3 ) 1.048 Å-3. With this value for F* the parameter Γ in eq 5.19 is Γ ) 1.790. With this information we can easily estimate the coefficients Rn of eq 5.15. Comparison with eq 5.19 allows us to identify R2 ) 1/3, R3 ) -1/9, and R4 ) (Γ/27 - 1/9) ) -0.0448. Using eqs 5.16-5.18 and these values for the Rn, we determine the coefficients θn of eq 5.14: θ0 ) -1/2, θ1 ) 1/6, and θ2 ) -0.0936. Notice that θ0 ) -0.5 is the value of θ that we would have obtained had we assumed θ independent of the dipolar strength parameter V. The naive use of eq 5.14

θ(V) ) -0.5 + 0.167V - 0.0936V2

In the last equation

with V ) 4.47 gives θ ) -1.63, which is far too negative (in fact, with this value of θ eq 5.11 gives negative 0). Presumably this failure occurs because the expansions in powers of V are only valid for small V. We have generated [1/1] and [0/2] Pade´ approximants for θ(V) based on expansion (5.21):

g3 ) θ2 + 2θ1 + 1 g2 ) 3 + 3θ2 + θ12 + 6θ1 g1 ) 1 - 3θ1 - 2θ12 - 3θ2 g0 ) θ1 - θ2 - θ12

(5.19)

where Γ ) [0.0444/(F*)2 + 7/4], F* ) Fd3, and d is the diameter of the dipolar hard spheres. We test the procedure by calculating θ and 0 for a fluid of symmetric dipolar dumbbells (SDD) at T ) 72 K and number density F ) 0.018 Å-3. The Lennard-Jones parameters in the u/jl(r) term of eq 1.4 are σ ) 3.341 Å and /kB ) 44.0 K, and are the same for both interaction sites. The dumbbell bond length is l ) 1.3364 Å, so that l /σ ) 0.4. The partial charges at the interaction sites (of the same magnitude at both molecular sites, but of opposite sign) are chosen such that the dipolar strength parameter is V ) 4.47 at the given T, F, and l . To be able to use the Rn coefficients from eq 5.19 in the determination of θ, we must decide which DHS model corresponds to our SDD model. We suggest to choose the DHS diameter d so that the hard sphere volume is equal to the volume of a SDD molecule. The condition is33

d ) ds[1 + 1.5(l /ds) - 0.5(l /ds)3]1/3

-0.5 - 0.1143V 1.0 + 0.5619V

(5.22)

-0.5 1.0 + 0.3333V - 0.0857V2

(5.23)

θ(V) )

Equations 5.16-5.18 suggest that we could use information of the dependence of 0 with V at small values of V to determine the coefficients θn, n ) 0, 1, 2, of eq 5.14. Such a dependence of 0 with V has been studied in detail for fluids of dipolar hard spheres (DHS),30,31 but not, to our knowledge, for interaction site model fluids.32 Because of this unfortunate lack of information, we are forced to apply eqs 5.16-5.18 based on coefficients Rn appropriate for DHS fluids. Clearly this approach will be most reliable when used to determine θ of two-site-polar ISM liquids that deviate from spherical shape only slightly. A particularly interesting representation for 0 of a DHS fluid as a function of V at small V may be extracted from the work of de Leeuw et al.30

0 ) 1 + V + V2/3 - V3/9 + (Γ/27 - 1/9) V4 + ...

(5.21)

(5.20)

where ds is the diameter of the hard dumbbell interaction site. We estimate ds from the Lennard-Jones parameters of the SDD model and the temperature with the method of Hsu, Chandler, and Lowden.34 For the specific parameters of our SDD model we get ds ) 3.336 Å. With this value for ds we calculate with eq 5.20 that the hard sphere diameter of the corresponding DHS

θ(V) )

At V ) 4.47 these equations give θ ) -0.288 (Pade´ [1/1]) and θ ) -0.517 (Pade´ [0/2]). Use of these two predictions for the mixing parameter θ in eq 5.11 gives 0 ) 6.83 (Pade´ [1/1]) and 0 ) 8.64 (Pade´ [0/2]). These numbers should be compared with 0 ) 9.4 calculated by molecular dynamics35 and the RISMHNC or ideal gas value 0 ) 5.47. The site-site correlation functions calculated with closure III for both values of θ are very similar to the correlation functions calculated under the RISM-HNC closure. The methodology outlined above is limited to ISM models that belong to the two-site-polar class and that, furthermore, have a molecular core shape that is nearly spherical. At present we continue investigating alternative ways for the determination of θ that can be applied to arbitrary interaction site models. Our primary focus is the determination of θ by demanding the consistency between results calculated with RISM-closure III following two different thermodynamic routes. 6. Summary We have presented three new dielectrically nontrivial closures for the RISM integral equation which are reminiscent of the familiar RISM-HNC closure (cf. eq 1.13). The closures are distinguished by the choice of the function τjl(r) that replaces the HNC function tjl(r). The particular choices of τjl(r) in closures I-III are inspired by an unusual (although equivalent) representation (eq 3.1) of the HNC tjl(r) in the case of a simple ionic mixture. We justify this strategy by the following considerations: 1. The HNC closure (in the context of the Ornstein-Zernike equation) complies with the Stillinger-Lovett second moment condition36 for ionic systems. 2. The fact that 0 of a polar liquid is determined by the second moment of the site-site correlation correlation functions hjl(r) (cf. eq 2.1). Like the molecular site-site closures recently introduced in the context of the RISM polaron theory for the hydrated electron37 and of the polymer RISM theory,38 closures I-III are nonlocal in r-space. Thus, the expression of these molecular

Closures for the RISM Integral Equation closures not only involves a particular pair of site-site functions {hjl(r), cjl(r)}, but also incorporates correlation functions hmn(r) associated with other pairs mn of interaction sites. The key feature in the τjl(r) functions of the new closures I-III that is responsible for the nontrivial predictions of 0 is the inclusion of the inverse intramolecular correlation function matrix Ω ˜ (k) (eq 3.6), which for many ISM models has the behavior shown in eq 3.7 at small k. This is in marked contrast with the other molecular closures for the RISM integral equation reported in the literature37,38 that explicitly incorporate the intramolecular correlation function matrix ω ˜ (k) instead of the inverse matrix Ω ˜ (k). Although the molecular closures discussed in refs 37 and 38 improve over some of the limitations of the site-site closures for the RISM integral equation, they appear problematic when applied to the polar ISM models (which have unshielded Coulombic interactions between interaction sites) considered in this study.39 We illustrate this in Appendix D for the Laria-Wu-Chandler closure,37 where we show that for small values of the wave vector k the matrix of RISM direct correlation functions c˜ (k) is dominated by a term proportional to k-4. A similar behavior can be demonstrated for the molecular RISM-HNC closure. The prediction by these closures that for polar liquids c˜ (k) ) k-4c-4 + ... is problematic; such behavior implies that for large values of the distance r the site-site direct correlation functions cjl(r) are directly proportional to r. The three closures give simple analytical formulas for 0 for polar liquids represented by ISM models with only two inequivalent interaction sites. Although limited to a restricted class of models, these formulas pinpoint a conceptually important limitation of closures I and II, namely their failure to capture the dependence of 0 on the thermodynamic state parameters and molecular features separately from the combination of variables in the dipolar strength parameter V. This limitation of closures I and II (as well as other technical difficulties) prompted us to formulate τjl(r) in closure III as a hybrid between the tjl(r) function of RISM-HNC and a slight variant τPjl (r) of the τjl(r) function of closure I. Thus, closure III incorporates a new parameter θ that gauges the extent tjl(r) and τPjl (r) contribute to the τjl(r) of closure III. It is through this parameter that closure III is able to incorporate the subtle dependence of 0 on the molecular features (size, shape, and charge distribution) specific to each ISM model. In section 5.4 we discussed a first attempt to determine θ for a two-site-polar model, so that closure III can be used to predict the dielectric constant of the liquid. We are currently exploring alternative strategies for the determination of θ, based on the requirement of consistency of the closure for two different thermodynamic routes, that we expect will be more general (i.e., applicable to arbitrary multisite polar liquids) and robust. Acknowledgment. We dedicate this work to Howard Reiss, who has long been concerned with the interplay between the geometric and electrostatic features of molecular modeling. We are grateful to Professor Martin Neumann, of Vienna University, for communicating to us his result for the dielectric constant of the symmetric dipolar dumbbell model liquid discussed in section 5.4. F.O.R. and G.S. gratefully acknowledge the support of the Division of Chemical Sciences, Office of Basic Energy Sciences, Office of Energy Research, U.S. Department of Energy. Appendix A: The Matrix Ω ˜ (k) In this appendix we derive general conclusions on the form of the inverse intramolecular correlation function matrix Ω ˜ (k)

J. Phys. Chem. B, Vol. 105, No. 47, 2001 11889 ˜ (k) is an ≡ [ω ˜ (k)]-1. For a molecule with p interaction sites ω p × p matrix with elements ω ˜ jl(k) ) δjl + (1 - δjl)j0(kl jl), where j0 is the spherical Bessel function of zero order. Given that j0(0) ) 1, for our purposes we can write ω ˜ (k) in the more convenient form

ω ˜ (k) ) Λ ˜ (k) + |1)(1|

(A.1)

where |1)(1| is the p × p dyadic matrix with all its elements equal to 1, while Λ ˜ (k) is an p × p matrix of elements

Λ ˜ jl(k) ) (1 - δjl) × [j0(kl jl) -1]

(A.2)

We can easily calculate the inverse of ω ˜ (k) from eq A.1 using the Sherman-Morrison formula41

{A + |x)(y|}-1 ) A-1 -

A-1|x)(y|A-1 1 + (y|A-1|x)

(A.3)

which is obviously valid as long as matrix A is invertible and (y|A-1|x) * -1. Applying eq A.3 to calculate [ω ˜ (k)]-1, we get

Ω ˜ (k) ) [Λ ˜ (k)]-1 -

˜ (k)]-1 [Λ ˜ (k)]-1|1)(1|[Λ 1 + (1|[Λ ˜ (k)]-1|1)

(A.4)

From this equation we may derive an explicit expression for the coefficient Ω-2 in the expansion of eq 3.7. For this purpose we note that

Λ ˜ (k) ) k2ω2 + O(k4)

(A.5)

which in turn implies (for ω2 invertible)

[Λ ˜ (k)]-1 ) k-2ω2-1 + O(k0)

(A.6)

When we substitute this expression into eq A.4, we find that the coefficient of the power k-2 in eq 3.7 is

Ω-2 ) ω2-1 -

ω2-1|1)(1|ω2-1

(A.7)

(1|ω2-1|1)

From this expression we infer the important relations

Ω-2|1) ) |0), (1|Ω-2 ) (0|

(A.8)

A more thorough examination of the above derivation is required to identify some problematic cases. Taking into account that [Λ ˜ (k)]-1 ) A ˜ (k)/D ˜ (k), where D ˜ (k) and A ˜ (k) are, respectively, the determinant and the adjoint of Λ ˜ (k), eq A.4 can be expressed as

Ω ˜ (k) )

{

A ˜ (k)|1)(1|A ˜ (k) 1 A ˜ (k) D ˜ (k) D ˜ (k) + (1|A ˜ (k)|1)

}

(A.9)

For molecules with p interaction sites we have the general expansions

D ˜ (k) ) k2pD2p + k2(p+1)D2(p+1) + k2(p+2)D2(p+2) + ... (A.10) and

A ˜ (k) ) k2(p-1)A2(p-1) + k2pA2p + k2(p+1)A2(p+1) + ... (A.11) When these expansions are substituted into eq A.9, we recover


Raineri and Stell

eq 3.7 with an expression for Ω-2 which is equivalent to eq A.7:

Ω-2 )

{

}

A2(p-1)|1)(1|A2(p-1) 1 A2(p-1) D2p (1|A2(p-1)|1)

(A.12)

In some cases (for instance, linear ISM models with more than two interaction sites, like a three-site model of acetonitrile) the condition

(1|A2(p-1)|1) ) 0

(A.13)

applies. For these cases it is easy to see that instead of eq 3.7 we have at small k

iteration the second moment

hjl,2 ) (1/2)(d2h˜ jl(k)/dk2)k)0

(C.1)

of the site-site correlation functions. This is required for the implementation of the short-range/long-range separation of matrix H ˜ ′′(k) according to Ng’s method. We have found that to calculate numerically the required second moments, the approximation eq C.5 below, which is based on the representation of the analytical h˜ jl(k) functions in terms of Cardinal functions,44 is quite satisfactory. To derive this approximation, we assume that h˜ jl(-k) ) h˜ jl(k). We can then represent44 h˜ jl(k) in the form N

∑

h˜ jl(k) = CN(h˜ jl,∆)(k) )

h˜ jl(n∆) S(n, ∆)(k) (C.2)

n)-N

Ω ˜ (k) ) k-4Ω-4 + k-2Ω-2 + k0Ω0 + ...

(A.14)

where44

S(n,∆)(k) )

with

Ω-4 ) -

A2(p-1)|1)(1|A2(p-1) D2p[D2p + (1|A2p|1)]

(A.15)

The integral equation closures proposed in this work fail when applied to liquids comprising molecules for which eq A.13 applies.

˜ (k)‚O˜(k)‚ω ˜ (k) Q ˜ (k) ) (I - ω ˜ (k)‚O˜(k)‚G)-1‚ω

(B.1)

where I is the p × p unit matrix and G ) FI. Given the dyadic form of O˜(k), cf. eq 2.14, we can calculate the inverse in eq B.1 with the help of the Sherman-Morrison formula (eq A.3). After some manipulations we obtain

(q|Q ˜ (k)|q) )

{ } () dk2

From eq 3.15 we derive

(q|ω ˜ (k)‚O˜(k)‚ω ˜ (k)|q) 1 + (4πβF/k2)(q|ω ˜ (k)|q)

(B.2)

Using again the dyadic expression of O˜(k) and the small-k expansion

(q|ω ˜ (k)|q) ) k2(µ2/3) + O(k4)

(B.3)

=-

k)0

2 π h˜ jl(k)0)

∆

3

-

4

N

∑

(-)nh˜ jl(n∆)

∆2n)1

n2

(C.4)

We notice that the compressibility relation14 imposes the requirement that all of the zero-moment coefficients h˜ jl(k)0) be equal to each other. When we recall that the matrix Ω-2 satisfies the conditions displayed in eq A.8, we see that to calculate the second moment coefficients hjl,2 needed for the computation of H′′-2 (cf. eq 5.8) we may safely ignore the first term in eq C.4, and use the simpler expression

hjl,2 ) -

2

N

∑

∆2n)1

(-)nh˜ jl(n∆) n2

(C.5)

Appendix D: Laria-Wu-Chandler Closure The Laria-Wu-Chandler (LWC) closure,37 for large separation r between interaction sites j and l in two different molecules, may be stated in the form

˜ (k)‚c˜ (k)‚ω ˜ (k)}jl ) F -1{ω ˜ (k)‚O˜(k)‚ω ˜ (k)}jl + Ljl(r) F -1{ω (D.1)

we finally derive

Q2 ≡ 4πβF2(q|Q2|q) )

(C.3)

with ∆ the k-grid spacing used in the numerical procedure. In eq C.2 h˜ jl(n∆) stands for the value of h˜ jl(k) at k ) n∆. The function CN(h˜ jl,∆)(k) defined by the second equality in eq C.2 is analytic in k. We can estimate44 d2h˜ jl(k)/dk2 by simply taking the second derivative of CN(h˜ jl,∆)(k) with respect to k. Exploiting the assumed even character of h˜ jl(k), we calculate

d2h˜ jl(k)

Appendix B: The Coefficient Q2

sin(π(k - n∆)/∆) π(k - n∆)/∆

-V2 1+V

(B.4)

which was used in section 4.1. Appendix C: Computational Details To solve the RISM integral equation under closure III, we have used the Picard iterative approach. Technical difficulties associated with the long-range nature of the cjl(r) and τjl(r) functions were avoided by means of Ng’s method.42,43 A unique issue arises in the solution of the RISM-closure III integral equation; namely the necessity of calculating at each

where we have used the explicit k-space notation F -1{ω ˜ (k)‚x˜ (k)‚ ω ˜ (k)}jl (with x either c or O) to indicate the jl matrix element of the r-space convolution between matrices ω, x, and ω. The function Ljl(r) in eq D.1 is defined by

(-)n Ljl(r) ) hjl(r) - ln gjl(r) ) [hjl(r)]n ng2 n

∑

(D.2)

In k-space eq D.1 becomes

ω ˜ (k)‚c˜ (k)‚ω ˜ (k) ) ω ˜ (k)‚O˜(k)‚ω ˜ (k) + L ˜ (k)

(D.3)

By left- and right-multiplying each term of eq D.3 with the inverse of the intramolecular correlation function matrix Ω ˜ (k)



(cf. eq 3.6), we obtain

c˜ (k) ) O˜(k) + Ω ˜ (k)‚L ˜ (k)‚Ω ˜ (k)

(D.4)

The small-k behavior of matrices O˜(k) and Ω ˜ (k) is shown in eqs 2.14 and 3.7: in both cases it is dominated by a term ˜ (k) behaves at small proportional to k-2. In contrast, matrix L values of the wave vector k as

L ˜ (k) ) k0L0 + k2L2 + O(k4)

(D.5)

For our purposes, the most important feature of this expansion is that the matrix coefficient L0, with elements

(L0)jl ) 4π

∫0

∞

dr r Ljl(r) ) 4π 2

∑ ng2

(-)n n

∫0∞dr r2[hjl(r)]n (D.6)

is not of the rank 1 form: L0 ) (constant) × |1)(1|. Therefore, Ω-2 ‚L0 * 0 and L0‚Ω-2 * 0. Taking into account the alluded small-k behavior of matrices O˜(k), Ω ˜ (k), and L ˜ (k), it follows from eq D.4 that under the LWC closure

c˜ (k) )k-4c-4 + k-2c-2 + k0c0 + O(k)

(D.7)

with the following expressions for the c-4 and c-2 matrix coefficients

c-4 ) Ω-2‚L0‚Ω-2 c-2 ) O-2 + Ω-2‚L0‚Ω0 + Ω0‚L0‚ Ω-2 + Ω-2‚L2‚Ω-2 (D.8) That any site-site correlation function c˜ jl(k) is dominated by a k-4 contribution implies that in r-space cjl(r) is directly proportional to r.40 References and Notes (1) Chandler, D.; Andersen, H. C. J. Chem. Phys. 1972, 57, 1930. (2) Chandler, D. Annu. ReV. Phys. Chem. 1978, 29, 441. (3) Chandler, D. In The Liquid State of Matter: Fluids Simple and Complex, Montroll, E., Lebowitz, J. L., Eds.; North-Holland: Amsterdam, 1982. (4) Gray, C. G.; Gubbins, K. E. Theory of Molecular Fluids; Oxford University Press: New York, 1984; Vol. I. (5) Monson, P. A.; Morris, G. P. AdV. Chem. Phys. 1990, 50, 451. (6) Hirata, F.; Rossky, P. J. J. Chem. Phys. Lett. 1981, 83, 329. Hirata, F.; Pettitt, B. M.; Rossky, P. J. J. Chem. Phys. 1982, 77, 509. Rossky, P. J. Annu. ReV. Phys. Chem. 1985, 36, 321. (7) Stell, G. In Statistical Mechanics, Part A: Equilibrium Techniques; Berne, B. J., Eds.; Plenum Press: New York, 1977, and references therein. (8) As discussed in detail by Stell and Zhou9 and by Zhou and Stell,10 the difference between eqs 1.5 and 1.9 can be understood by noting that interaction site molecules can be thought of as representing the complete association limit of associating atoms centered at the interaction sites. For a fluid of associating particles, the full-association limit is a highly singular state, in much the same sense that the critical point is a highly singular thermodynamic state. In both cases, what one means by “large r” is an r that is significantly larger than a characteristic correlation length that goes to infinity at the singular state (as discussed in these references). Just as closures consistent with eq 1.6 yield mean-field critical behavior that is typically in substantial quantitative error, closures that satisfy eq 1.5 give a “trivial” result11,12 for the dielectric constant of the liquid, namely the ideal gas value given by the right-hand side of eq 1.7. (9) Stell, G.; Zhou, Y. J. Chem. Phys. 1989, 91, 4861. See also: Stell, G.; Zhou, Y. Fluid Phase Equilib. 1992, 79, 1. (10) Zhou, Y.; Stell, G J. Chem. Phys. 1993, 98, 5777. (11) Høye, J. S.; Stell, G. J. Chem. Phys. 1976, 65, 18. (12) Chandler, D. J. Chem. Phys. 1977, 67, 1113. (13) Cummings, P. J.; Stell, G. Mol. Phys. 1981, 44, 529. (14) Cummings; P. J.; Stell, G. Mol. Phys. 1982, 46, 383. (15) (a) Rossky, P. J.; Pettitt, B. M.; Stell, G. Mol. Phys. 1983, 50, 1263. (b) Yu, H.-A.; Roux, B.; Karplus, M. J. Chem. Phys. 1990, 92, 5020.

(16) Raineri, F. O.; Resat, H.; Friedman, H. L. J. Chem. Phys. 1992, 96, 3068. (17) Raineri, F. O.; Perng, B.-C.; Friedman, H. L. Electrochim. Acta 1997, 42, 2749. (18) Equation 2.1 was first reported in ref 11. Subsequent derivations of this equation are given in: Høye, J. S.; Stell, G. J. Chem. Phys. 1977, 67, 1776. The result is generalized to polar, polarizable site-site models in: Stell, G.; Patey, G. N.; Høye, J. S. AdV. Chem. Phys. 1981, 48, 183. Derivation of these results is also given by Chandler in ref 12. (19) Equation 2.9 was first derived, in somewhat different notation, in ref 11. It was also reported in ref 12. (20) This observation also extends to the τjl(r) functions in the integral equation closures II and III that are introduced in this work. (21) It is possible to reduce the dimension of the correlation function matrices in the RISM integral eq 1.1 and other equations by recognizing the chemical equivalence of interaction sites. This may be seen most simply by introducing the matrix w ˜ (k) of symmetry-reduced intramolecular correlation functions, with elements

1 1 w˜ JL(k) ) s ∑ ω ˜ (k) ) s ∑ ω ˜ (k) ) w ˜ LJ(k) nJ j∈J jl nL l∈L jl In this formula the indices J and L refer to the inequivalent types of molecular interaction sites, while nJ and nL are the corresponding symmetry numbers of sites of type J and L in the molecule. The range {j ∈ J} in the summation symbol indicates that the sum extends over all the equivalent interaction sites that belong to the J-type of interaction site. With this transformation the original RISM equation may be expressed in symmetryreduced form as follows: h˜ (k) ) w ˜ (k)‚C ˜ (k)‚w ˜ (k) + w ˜ (k)‚C ˜ (k)‚G‚h˜ (k), where the new symmetry-reduced matrix C ˜ (k) of RISM direct correlation functions has elements C ˜ JL(k) ) nJc˜ JL(k)nL. Each matrix in this equation has dimension P × P, where P is the number of inequivalent types of interaction sites in the ISM model. The density matrix G is diagonal, and each of its nonzero elements equals the number density of the liquid. All the simplifications alluded to in the text for certain classes of interaction site models follow from the systematic use of the symmetry-reduced formalism. For clarity, however, the equations in the text correspond to the RISM formalism in which all interaction sites are considered explicitly. (22) Notice that eq A.7 indicates that, for ISM models that do not belong to the two-site class, the rank of matrix Ω-2 is greater than 1. (23) Morris, G. P. Mol. Phys. 1982, 47, 833. (24) Morris, G. P.; Cummings, P. T. Mol. Phys. 1982, 45, 1099. (25) Berendsen, H. J. C.; Grigera, J. R.; Straatsma, T. P. J. Phys. Chem. 1987, 91, 6269. (26) Ho¨chtl, P.; Boresch, S.; Bitomsky, W.; Steinhauser, O. J. Chem. Phys. 1998, 109, 4927. (27) Actually, for models outside the two-site-polar ISM class eq 1.9 is not satisfied either by the correct RISM direct correlation function matrix or by the direct correlation function matrix generated by any of the new closures proposed in this work. This is to be contrasted with the modified RISM-HNC closure, eqs 1.11 and 1.12, which indeed assumes the validity of eq 1.9 for eVery ISM model of polar liquids. (28) Liu, H.; Mu¨ller-Plathe, F.; van Gunsteren, W. J. Am. Chem. Soc. 1993, 117, 4363. (29) Skaf, M. S. In Electron and Ion Transfer in Condensed Media; Kornyshev, A. A., Tosi, M., Ulstrup, J., Eds.; World Scientific: Singapore, 1997. (30) de Leeuw, S. W.; Perram, J. W.; Quirke, N.; Smith, E. R. Mol. Phys. 1981, 42, 1197. (31) Tani, A.; Henderson, D.; Barker, J. A.; Hetch, C. E. Mol. Phys. 1983, 48, 863. (32) Fluids with spherically symmetric short-range interactions but more complicated charge distributions have also been investigated; see for example: Goldman, S.; Joslin, C. J. Phys. Chem. 1993, 97, 12349. (33) See section 5.5.1 of ref 4. (34) Hsu, C. S.; Chandler, D.; Lowden, L. J. Chem. Phys. 1976, 14, 213. (35) Neumann, M. Private communication. (36) Stillinger, F. H.; Lovett, R. J. Chem. Phys. 1968, 49, 1991. (37) Laria, D.; Wu, D.; Chandler, D. J. Chem. Phys. 1991, 95, 4444. (38) Schweizer, K. S.; Yethiraj, A. J. Chem. Phys. 1993, 98, 9053. (39) The Laria-Wu-Chandler closure has been applied to the study of the structure functions of polyelectrolyte solutions (where only the polyelectrolyte molecules are considered explicitly) with excellent results. See: Shew, C.-Y.; Yethiraj, A. J. Chem. Phys. 1997, 106, 5706. Notice, however, that the site-site interactions for such system involve a shielded Coulombic contribution of the form qjql exp(-κr)/r, as opposed to the pure Coulombic term in eq 1.4. To our knowledge, the molecular closures introduced in refs 37 and 38 have not been used to calculate the structure functions of the class of ISM polar liquids [characterized by the site-site potentials (1.4)] considered in the present study.

11892 J. Phys. Chem. B, Vol. 105, No. 47, 2001 (40) Friedman, B. Lectures on Applications-Oriented Mathematics; John Wiley and Sons: New York, 1991; Chapter 1. (41) Lancaster, P.; Tismenetsky, M. The Theory of Matrices, 2nd ed.; Academic Press: Boston, 1985.

Raineri and Stell (42) Ng, K.-C. J. Chem. Phys. 1974, 61, 2680. (43) Morris, G. P.; MacGowan, D. Mol. Phys. 1986, 58, 745. (44) Stenger, F. Numerical Methods Based on Sinc and Analytic Functions; Springer: New York, 1993.

Dielectrically Nontrivial Closures for the RISM Integral Equation - The

Recommend Documents