A new family of bridge functions for electrolyte solutions
We present a new set of closures for restricted models of electrolyte solutions at the McMillan-Mayer level that improve upon the Hypernetted Chain prediction for the ion-ion pair correlation functions. The improvement is accomplished by proposing simple functional forms for the bridge functions...
Збережено в:
| Дата: | 2001 |
|---|---|
| Автори: | , |
| Формат: | Стаття |
| Мова: | English |
| Опубліковано: |
Інститут фізики конденсованих систем НАН України
2001
|
| Назва видання: | Condensed Matter Physics |
| Онлайн доступ: | http://dspace.nbuv.gov.ua/handle/123456789/120520 |
| Теги: |
Додати тег
Немає тегів, Будьте першим, хто поставить тег для цього запису!
|
| Назва журналу: | Digital Library of Periodicals of National Academy of Sciences of Ukraine |
| Цитувати: | A new family of bridge functions for electrolyte solutions / F.O. Raineri, G. Stell // Condensed Matter Physics. — 2001. — Т. 4, № 4(28). — С. 621-642. — Бібліогр.: 16 назв. — англ. |
Репозитарії
Digital Library of Periodicals of National Academy of Sciences of Ukraine| id |
irk-123456789-120520 |
|---|---|
| record_format |
dspace |
| spelling |
irk-123456789-1205202017-06-13T03:05:14Z A new family of bridge functions for electrolyte solutions Raineri, F.O. Stell, G. We present a new set of closures for restricted models of electrolyte solutions at the McMillan-Mayer level that improve upon the Hypernetted Chain prediction for the ion-ion pair correlation functions. The improvement is accomplished by proposing simple functional forms for the bridge functions and the specification of certain adjusting parameters according to several criteria. Under the new closures, and unlike the HNC case, the “sum” direct correlation function, which is crucial for determining the stability of the solution with respect to phase separation, remains finite at thermodynamic states along the spinodal and at the critical point. Для спрощених моделей розчинів електролітів ми представляємо новий набір замикань на рівні МакМіллана–Маєра, який покращує передбачення для парних кореляційних функцій іон-іон, отриманих у гіперланцюговому наближенні. Покращення здійснюється введенням простих функціональних форм для місткових функцій і деяких підгоночних параметрів, які відповідають різним критеріям. При нових умовах замикання, на противагу до гіперланцюгового наближення, пряма кореляційна функція “сума”, що є важливою для визначення стійкості розв’язку по відношенню до фазового розділення, залишається скінченою при термодинамічних станах уздовж спінодалі і в критичній точці. 2001 Article A new family of bridge functions for electrolyte solutions / F.O. Raineri, G. Stell // Condensed Matter Physics. — 2001. — Т. 4, № 4(28). — С. 621-642. — Бібліогр.: 16 назв. — англ. 1607-324X PACS: 61.20.Gy, 61.20.Qg, 64.60.Fr, 64.70.Ja, 64.75.+g DOI:10.5488/CMP.4.4.621 http://dspace.nbuv.gov.ua/handle/123456789/120520 en Condensed Matter Physics Інститут фізики конденсованих систем НАН України |
| institution |
Digital Library of Periodicals of National Academy of Sciences of Ukraine |
| collection |
DSpace DC |
| language |
English |
| description |
We present a new set of closures for restricted models of electrolyte solutions
at the McMillan-Mayer level that improve upon the Hypernetted Chain
prediction for the ion-ion pair correlation functions. The improvement is accomplished
by proposing simple functional forms for the bridge functions
and the specification of certain adjusting parameters according to several
criteria. Under the new closures, and unlike the HNC case, the “sum” direct
correlation function, which is crucial for determining the stability of the
solution with respect to phase separation, remains finite at thermodynamic
states along the spinodal and at the critical point. |
| format |
Article |
| author |
Raineri, F.O. Stell, G. |
| spellingShingle |
Raineri, F.O. Stell, G. A new family of bridge functions for electrolyte solutions Condensed Matter Physics |
| author_facet |
Raineri, F.O. Stell, G. |
| author_sort |
Raineri, F.O. |
| title |
A new family of bridge functions for electrolyte solutions |
| title_short |
A new family of bridge functions for electrolyte solutions |
| title_full |
A new family of bridge functions for electrolyte solutions |
| title_fullStr |
A new family of bridge functions for electrolyte solutions |
| title_full_unstemmed |
A new family of bridge functions for electrolyte solutions |
| title_sort |
new family of bridge functions for electrolyte solutions |
| publisher |
Інститут фізики конденсованих систем НАН України |
| publishDate |
2001 |
| url |
http://dspace.nbuv.gov.ua/handle/123456789/120520 |
| citation_txt |
A new family of bridge functions for
electrolyte solutions / F.O. Raineri, G. Stell // Condensed Matter Physics. — 2001. — Т. 4, № 4(28). — С. 621-642. — Бібліогр.: 16 назв. — англ. |
| series |
Condensed Matter Physics |
| work_keys_str_mv |
AT rainerifo anewfamilyofbridgefunctionsforelectrolytesolutions AT stellg anewfamilyofbridgefunctionsforelectrolytesolutions AT rainerifo newfamilyofbridgefunctionsforelectrolytesolutions AT stellg newfamilyofbridgefunctionsforelectrolytesolutions |
| first_indexed |
2025-07-08T18:01:06Z |
| last_indexed |
2025-07-08T18:01:06Z |
| _version_ |
1837102708800094208 |
| fulltext |
Condensed Matter Physics, 2001, Vol. 4, No. 4(28), pp. 621–642
A new family of bridge functions for
electrolyte solutions
F.O.Raineri, G.Stell
Department of Chemistry, State University of New York at Stony Brook,
Stony Brook, New York 11794–3400, USA
Received November 5, 2001
We present a new set of closures for restricted models of electrolyte solu-
tions at the McMillan-Mayer level that improve upon the Hypernetted Chain
prediction for the ion-ion pair correlation functions. The improvement is ac-
complished by proposing simple functional forms for the bridge functions
and the specification of certain adjusting parameters according to several
criteria. Under the new closures, and unlike the HNC case, the “sum” di-
rect correlation function, which is crucial for determining the stability of the
solution with respect to phase separation, remains finite at thermodynamic
states along the spinodal and at the critical point.
Key words: electrolyte solutions, restricted model, structure,
thermodynamics, bridge functions, thermodynamic consistency
PACS: 61.20.Gy, 61.20.Qg, 64.60.Fr, 64.70.Ja, 64.75.+g
1. Introduction
The accurate prediction of the structure and thermodynamic properties of elec-
trolyte solutions over a wide range of thermodynamic states has been the focus of
an enormous number of studies. In spite of these efforts, there are still a number of
issues that require more development, especially those concerning the description of
the structure and thermodynamic properties of an electrolyte solution close to criti-
cality (i. e., the critical point associated with the unmixing of an ionic solution into
two liquid phases of different electrolyte concentration). Another relevant issue is
the need of an improved description of the structure and thermodynamic properties
for aqueous solutions of 2–2 electrolytes (and also for electrolytes of higher charge)
at low concentrations.
The vast majority of the efforts have focussed on the calculation of the structure
and other properties of electrolyte solutions by statistical mechanical methods at
the McMillan-Mayer level [1]: integral equations and computer simulation. In the
first approach the ion-ion pair correlation functions gjl(r) are calculated by solving
an appropriate set of Ornstein-Zernike (OZ) equations, that is supplemented with
c© F.O.Raineri, G.Stell 621
F.O.Raineri, G.Stell
a set of additional relations between the indirect {hjl(r)} and the direct {cjl(r)}
correlation functions [2,3]. Because of their role in setting the mathematical problem
determinate, the additional relations are referred to as closure relations or integral
equation closures.
Among the closure relations proposed for ionic systems, the hypernetted chain
(HNC) closure [2,3] has probably been the closure most thoroughly investigated.
The HNC closure predicts the structural and thermodynamic properties of 1–1 elec-
trolytes quite accurately up to concentrations greater that 1mol dm−3. For 2–2 elec-
trolytes, however, this closure predicts a spurious peak in the like-sign (i. e., cation-
cation or anion-anion) correlation functions, which is absent in the correlation func-
tions extracted from computer simulation calculations [4]. At the same time, results
obtained from molecular dynamics computer simulations indicate that the height
(gmax
12 )MD of the peak in the cation-anion correlation function is considerably larger
than the corresponding quantity (gmax
12 )HNC obtained with the HNC closure [5,6]. A
somewhat less known problem is that the direct correlation functions c jl(r) under
the HNC closure become long-ranged as the state of the solution approaches the
critical point (or, more generally, the stability limit of the solution with respect to
phase separation) [7].
Recently we proposed an integral equation closure for general models of binary
electrolytes that gives direct correlation functions cjl(r) that remain of finite range
at every thermodynamic state [7]. The new closure is formulated by decomposing
the density fluctuation of each ionic species into contributions to two mutually or-
thogonal types of fluctuations: the charge density fluctuations and the neutral den-
sity fluctuations. The source in the ion-ion direct correlation functions responsible
for the long-range behaviour close to the critical point is eliminated in the closure
by appropriately choosing the bridge function. The particular form of the bridge
function, which is reminiscent of the Percus-Yevick (PY) type bridge function [see
equation (18) below], explains the name, HNC/PY, assigned to the closure. A short-
coming of the proposed HNC/PY closure is that it systematically gives (gmax
12 )HNC/PY
smaller than (gmax
12 )HNC and (gmax
12 )MD.
The purpose of this paper is to propose and investigate a new family of integral
equation closures for electrolyte solutions that are constructed using the HNC/PY
closure as the starting building block. The closures are obtained by simply extending
the bridge function associated with the HNC/PY closure with supplementary terms
of carefully selected form. In particular, the functional forms of the additional bridge
function terms are extracted from the asymptotic behaviour of the direct correlation
functions under the HNC/PY closure. Clearly this strategy guarantees the property
(shared with the HNC/PY closure) that the new closures remain of finite range close
to the critical point.
Although the HNC/PY closure was originally formulated for general solution
models of binary electrolytes at the McMillan-Mayer level [7], here we confine the
presentation to the simpler case of solutions of restricted or symmetric binary elec-
trolyte models (see the next section for details). Furthermore, with respect to the
performance of the new closures, in this work we consider only temperature and
622
Bridge functions for electrolyte solutions
concentration states where the electrolyte solution is stable.
The outline of the rest of the paper is as follows. In section 2 we discuss some
theoretical features of the HNC and the HNC/PY integral equation closures. Unlike
in the original presentation [7], here we use the language of the sum and difference
ion-ion correlation functions that is so convenient for problems involving symmet-
ric binary electrolytes. In this section we also introduce the new integral equation
closures HNC/λ1 and bD/λ1 by adding new terms to the bridge functions in the
HNC/PY closure. In section 3 we report calculations that test the performance of
the HNC/λ1 closure for models of aqueous solutions of a 2–2 electrolyte at 298.15K.
The results are compared with the results obtained from calculations with the HNC
and HNC/PY closures, as well as with the results obtained from molecular dynam-
ics computer simulations of the same solution models by Smith et al. [5]. Finally,
in section 4 we briefly discuss the results obtained under the HNC/λ1 closure, and
suggest alternatives for further improvement (the bD/λ1 closure and another exten-
sion, the HNC/λ2 closure). We also compare our approach with previous work on
the development of bridge functions appropriate for electrolyte solutions.
2. Theory
Our starting point is the well known McMillan-Mayer level expression [2,3]
hjl(r) = exp(−βūjl(r) + tjl(r) + bjl(r) ) − 1, j, l = 1, 2 (1)
for the ion-ion correlation functions, where the indices j and l refer to the ionic
species: 1 for the cations and 2 for the anions. As usual, β = (kBT )
−1 has the
meaning of the inverse of the temperature expressed in energy units (kB is the
Boltzmann constant).
We recall that the indirect correlation function hjl(r) is related to the ion-ion
radial distribution function gjl(r) [which is proportional to the probability that two
ions, one of species j and the other of species l, are separated by a distance r in the
ionic solution] according to the equation hjl(r) = gjl(r)− 1. The other functions in
equation (1) are ūjl(r), the solvent-averaged potential between two ions of species
j and l separated by a distance r at infinite dilution [1,2]; tjl(r) ≡ hjl(r) − cjl(r),
in which cjl(r) is the ion-ion direct correlation function; and finally the ion-ion
bridge function bjl(r) [2,3]. The sets of functions {hjl(r)} and {cjl(r)} are further
interrelated through (see below) the Ornstein-Zernike equations [2,3].
As mentioned in the introduction, in this study we focus on the case of a symmet-
ric (or restricted) model of an electrolyte solution. In such a model the electrolyte
dissociates into equal numbers of cations and anions. Furthermore, the cation and
anion species are entirely equivalent, except for their charges z1e and z2e, which are
of the same magnitude but of opposite signs (z1 = −z2; e is the charge of a proton).
Quite generally, then, we write the McMillan-Mayer level ion-ion solvent averaged
potentials ūjl(r) as the sum of a short-range contribution ū∗(r), that is the same
for every jl ionic pair-type, and a long-range Coulombic contribution whose sign
623
F.O.Raineri, G.Stell
depends on whether the labels j and l refer to the same or to different ionic species:
ūjl(r) = ū∗(r) + (−)j+l z
2e2
ǫ0r
, j, l = 1, 2 . (2)
It should be noted that in the McMillan-Mayer level description, the Coulomb in-
teraction between the charges is attenuated by the dielectric constant ǫ0 of the pure
solvent. A particular but important case of equation (2) is the restricted primitive
model, for which ū∗(r) is a hard sphere potential. In this study, however, the dis-
cussion will not be limited to that model, but instead it will consider more general
forms of ū∗(r).
For the restricted electrolyte models represented by equation (2) it is advanta-
geous to express any ion-ion function Fjl [such as any of the functions that occur
in equation (1)] in terms of the “sum” (subscript S) and difference (subscript D)
functions:
FS(r) ≡ [F11(r) + F12(r) ] / 2 , FD(r) ≡ [F11(r)− F12(r) ] / 2 . (3)
With this notation, equation (1) for jl =11, 12, 21, and 22 may be cast in the
more succinct form
hS(r) = e−βū∗(r)+tS (r)+bS(r) coshD(r) − 1 , (4)
hD(r) = e−βū∗(r)+tS (r)+bS(r) sinhD(r) , (5)
where we have introduced the auxiliary function
D(r) ≡ −βūD(r) + tD(r) + bD(r) . (6)
Notice that because ū∗(r) is the same function for every jl pair, we have used in
writing equations (4) and (5) that ūS(r) = ū∗(r), the short-range part of the ion-ion
potential energy of interaction. Correspondingly, we note that the nature of ūD(r)
is purely Coulombic: ūD(r) = z2e2/ǫ0r .
Equations (4) and (5) may also be rearranged to display the sum and the differ-
ence direct correlation functions explicitly:
cS(r) = e−βū∗(r)+tS (r)+bS(r) coshD(r) − 1 − tS(r) , (7)
cD(r) = e−βū∗(r)+tS (r)+bS(r) sinhD(r) − tD(r) . (8)
In terms of the Fourier transforms of the sum and difference functions, the
Ornstein-Zernike equations at the McMillan-Mayer level are
h̃S(k) = c̃S(k) + c̃S(k) ρt h̃S(k) , (9)
h̃D(k) = c̃D(k) + c̃D(k) ρt h̃D(k) , (10)
624
Bridge functions for electrolyte solutions
where ρt = ρ1 + ρ2 is the total number density of ions in the solution (for the
restricted electrolyte systems discussed here ρ1 = ρ2 applies). We recall that the
sum/difference decoupling suggested by equations (9) and (10) is only apparent, as
the set of functions {hS, cS} and {hD, cD} are still entangled by the relations (7)
and (8).
2.1. HNC closure
In the set of equations (7)–(10), which provide the means for the calculation
of the structure functions of the electrolyte solution, we still have to specify the
sum and difference bridge functions bS(r) and bD(r). To do so is one of the goals
of this study. In this section, however, we discuss the closure relations that result
from the simplest possible approximation for these functions, namely when both
bridge functions can be ignored (i. e., bS(r) = bD(r) = 0). This is the well-known
Hypernetted Chain (HNC) approximation, for which equations (4) and (5) simplify
to
hS(r) = e−βū∗(r)+tS (r) coshD0(r) − 1 , (11)
hD(r) = e−βū∗(r)+tS (r) sinhD0(r) , (12)
while the auxiliary function D(r) [cf. equation (6)] is given by
D0(r) = −βūD(r) + tD(r) . (13)
It is interesting to analyze here a shortcoming of the HNC closure that is not
widely appreciated. This may be seen most straightforwardly from equations (11)
and (12). These equations can be formally solved for the function D0(r):
D0(r) = arctanh{χ(r) } , (14)
in terms of the auxiliary function
χ(r) ≡
hD(r)
1 + hS(r)
. (15)
Equation (14) for D0(r) can then be substituted back into equation (11). Recall-
ing that in the resulting equation tS(r) = hS(r) − cS(r) in the argument of the
exponential, we can then solve for the sum combination cS(r) of the ion-ion direct
correlation functions. We obtain, for values of r large enough that we can make the
approximation exp(−βū∗(r)) ≃ 1,
cS(r) = ln
{
ehS(r)
1 + hS(r)
}
−
1
2
ln
{
1 − [χ(r) ]2
}
. (16)
The particular form of the second term in this equation follows straightforwardly
from the identity cosh arctanh(x) = 1/
√
(1− x2).
625
F.O.Raineri, G.Stell
The feature, revealed by equation (16), that cS(r) under the HNC closure is
expressed by a sum of two separate contributions, one that depends only on hS(r)
while the other depends only on χ(r), is the source of difficulties when the closure is
used to examine the material stability limit of the electrolyte solution. We recall [8]
that the stability of the solution against phase separation (unmixing) is determined
exclusively by the sum combination cS(r) of the direct correlation functions. (Briefly,
the solution is stable or metastable against phase separation as long as 1−ρ t c̃S(0) >
0, where c̃S(0) denotes the Fourier transform of cS(r) specialized at k = 0). On
the other hand, the signature of the states on the spinodal line (and, of course,
at the critical point) is that at those thermodynamic states the sum combination
hS(r) becomes infinitely long ranged (or, equivalently, h̃S(k) diverges as k → 0 at
these states). In contrast, the range of χ(r), which is proportional to the difference
combination hD(r), remains finite at the stability limit of the electrolyte solution
[9].
At large values of the interionic distance r, where both hS(r) and χ(r) are small,
we can approximate cS(r) in equation (16) by the first few terms of its power series
expansion with respect to hS(r) and χ(r):
cS(r) =
1
2
[hS(r)]
2 −
1
3
[hS(r)]
3 + · · ·
+
1
2
[χ(r) ]2 +
1
4
[χ(r) ]4 + · · · (17)
When the thermodynamic state approaches the stability limit and, correspondingly,
hS(r) becomes more and more long-ranged, equation (17) reveals that, under the
HNC closure, the sum combination cS(r) (due to the presence of terms with “pure”
[hS(r)]
n powers) also becomes infinitely long-ranged .
2.2. HNC/PY closure
As a first step in designing an integral equation relation that improves over the
HNC closure, it is natural to introduce into the closure a bridge function bS(r)
that eliminates from cS(r) the terms that contain pure [hS(r)]
n powers of the sum
combination function hS(r). This can be accomplished most simply by choosing
bS(r) of the form [7]
bS(r) = ln[ 1 + tS(r) ] − tS(r) . (18)
We still insist on the approximation bD(r) = 0. When these approximations are
introduced into the general equations (4) and (5) they become
hS(r) = e−βū∗(r) [ 1 + tS(r) ] coshD0(r) − 1 , (19)
hD(r) = e−βū∗(r) [ 1 + tS(r) ] sinhD0(r) , (20)
with D0(r) given by equation (13).
626
Bridge functions for electrolyte solutions
In the next section we show that the choice of equation (18) for the sum com-
bination bS(r) bridge function makes cS(r) under this closure short-ranged (which
is the correct behaviour) at the critical point or any other state along the spinodal
line of the electrolyte solution.
Equation (18) is the well-known form that the bridge function for a pure sub-
stance takes under the Percus-Yevick (PY) integral equation closure. In fact, the
PY flavour of this closure becomes immediately obvious in equation (19) if we ig-
nore the charge-fluctuation factor coshD0(r). On the other hand, through the choice
bD(r) = 0 for the difference bridge function combination, this closure retains all the
features of the HNC closure in what concerns the difference combination direct cor-
relation function cD(r). For these reasons we use the acronym HNC/PY to refer
to this closure. A version of this closure for general models (i. e., not restricted) of
electrolyte solutions has been reported [7].
As the calculations reported in section 3 show, the performance of the HNC/PY
closure for dilute solutions of 2–2 restricted electrolytes is somewhat disappointing,
with the HNC closure giving (broadly speaking) better results. The goal of this
paper is, using the HNC/PY closure as the foundation, to propose a new family of
improved integral equation closures for solutions of symmetric electrolytes.
2.3. Beyond the PY form for the sum bridge function bS(r)
To obtain improved results but, at the same time to preserve the feature that
cS(r) is of finite range at every thermodynamic state, we slightly modify the PY
form of the sum combination bS(r) of the bridge functions [cf. equation (18)]. We do
this by adding a corrector function λ(r) to the argument of the logarithmic term:
bS(r) = ln[ 1 + tS(r) + λ(r) ] − tS(r) . (21)
For generality we do not assume at this point that bD(r) = 0, but instead we leave
this function momentarily unspecified.
The general equations (4) and (5) take the form
hS(r) = e−βū∗(r) [ 1 + tS(r) + λ(r) ] coshD(r) − 1 , (22)
hD(r) = e−βū∗(r) [ 1 + tS(r) + λ(r) ] sinhD(r) , (23)
where D(r) was introduced in equation (6).
We follow now the discussion in section 2.1. We begin by solving equations (22)
and (23) for the auxiliary function D(r):
D(r) = arctanh [χ(r) ] , (24)
where χ(r) is given in equation (15).
At interionic distances such that exp(−βū∗(r)) ≃ 1, equations (22) and (24)
may be solved for the sum combination cS(r) of direct correlation functions, with
the result
cS(r) = λ(r) + ( 1 + hS(r) )
{
1 −
√
1− [χ(r) ]2
}
. (25)
627
F.O.Raineri, G.Stell
To derive this equation we have used the identity
sech arctanh [χ(r)] =
√
1− [χ(r) ]2 .
Equation (25) shows that the contribution λ(r) coming from the sum bridge
function bS(r) is separated from the remaining terms. Except for the multiplicative
factor (1 + hS(r)), these terms depend only on the function χD(r), which is short-
ranged at every thermodynamic state. Furthermore, the equation does not explicitly
involve the bridge function bD(r): this expression of cS(r) in terms of λ(r), hS(r)
and χ(r) is not affected by whichever approximation we choose for bD(r).
At large values of r, where χ(r) is small, we can approximate equation (25) for
cS(r) by
cS(r) = λ(r) + hD(r)
{
1
2
χ(r) +
1
8
[χ(r) ]3 + · · ·
}
, (26)
in which hS(r) does not appear explicitly
Leaving aside for a moment the function λ(r), equation (26) shows that at large
r the sum combination cS(r) involves terms from which pure powers [hS(r)]
n are
absent; such powers of hS(r) (that originate from the expansion of the denominator
of χ(r)) are always multiplied by powers [hD(r)]
m of the difference combination
function hD(r). The latter function is of finite range at every thermodynamic state
[9]. When λ(r) = 0, equation (26) coincides with the large-r behaviour of cS(r) under
the HNC/PY closure (cf. section 2.2). This demonstrates the statement made in that
section that cS(r) under the HNC/PY closure is of finite-range at the stability limit
of the ionic solution.
In view of the extreme complexity involved in the direct calculation of λ(r), here
we simply consider this function to be a functional of the functions χ(r) and hD(r)
(or, equivalently, given equation (15), a functional of hS(r) and hD(r)). Furthermore,
we propose that this functional has a form that closely resembles the first terms in
the asymptotic expansion of cS(r) under the HNC/PY closure [which is given by
equation (26) after making λ(r) = 0]. With these assumptions it follows that the
simplest candidate for λ(r) has the form
λ(r) =
1
2
BS hD(r)χ(r) , (27)
where the value of the bridge parameter BS needs to be specified.
Replacing this form of λ(r) back into equation (26) gives the asymptotic be-
haviour of the sum direct correlation function cS(r) at large r
cS(r) = (AS/2) hD(r)χ(r) + (1/8) hD(r) [χ(r) ]
3 + · · · , (28)
with the parameter AS is related to BS by AS = 1 +BS.
It is clear that this approach for constructing λ(r), that borrows from the asymp-
totic expansion of cS(r) under HNC/PY, guarantees that cS(r) under the new closure
[that incorporates λ(r)] remains of finite range at the stability limit of the electrolyte
solution. It is also worth noting that the strategy just outlined for the construction of
λ(r) is completely independent on which approximation we choose for the difference
bridge function bD(r).
628
Bridge functions for electrolyte solutions
2.4. Approximation for the difference bridge function bD(r)
As in the case of λ(r), we propose to construct an approximation for the difference
bridge function bD(r) by considering the first few terms of the asymptotic form of
cD(r) under the HNC/PY approximation.
From equations (6) and (24) we find for c∗D(r) ≡ cD(r)+βūD(r) (the short-range
part of the difference direct correlation function)
c∗D(r) = bD(r) + hD(r) −
1
2
ln
{
1 + χ(r)
1 − χ(r)
}
. (29)
The last term is, of course, equivalent to arctanh {χ(r) }. This expression is explicitly
linear in hD(r). Expanding the logarithmic term in powers of χ(r) we obtain
c∗D(r) = bD(r) + hD(r) − χ(r) −
1
3
[χ(r) ]3 −
1
5
[χ(r) ]5 + · · · . (30)
The second and third terms can be combined together, taking into account equa-
tion (15), so that the large-r asymptotic expansion of c∗D(r) becomes
c∗D(r) = bD(r) + hS(r)χ(r) −
1
3
[χ(r) ]3 + · · · . (31)
It is important to note that this result applies irrespective of any approximation
made to bS(r).
By taking bD(r) = 0 in equation (31) we obtain the asymptotic expansion of
cD(r) under any closure that approximates the difference bridge function to zero,
like the HNC or HNC/PY closures. We may generate an approximate expression
for bD(r) by proceeding as we did for constructing λ(r) in the last section; i. e., we
assume that bD(r) is proportional to the first term of the asymptotic expansion of
c∗D(r) under the HNC/PY closure. We then have
bD(r) = BD hS(r)χ(r) , (32)
in terms of the bridge parameter BD. With this approximation c∗D(r) [equation (31)]
behaves at large r as
c∗D(r) = AD hS(r)χ(r) −
1
3
[χ(r) ]3 + · · · , (33)
where AD = 1+BD. Clearly, this expansion does not contain terms with pure powers
[hS(r)]
n; such powers appear to be always multiplied by powers [hD(r)]
m that are
of finite range for every thermodynamic state.
Finally, it is straightforward to show from equations (28) and (33), together with
the expected small-k behaviour of the Fourier transforms of hS(r) and hD(r)
h̃S(k) = k0 h
(0)
S + k2 h
(2)
S + O(k4) , (34)
h̃D(k) = k0 h
(0)
D + k2 h
(2)
D + O(k4) , (35)
that the direct correlation functions cS(r) and cD(r) under the closures proposed in
this section have the form that is required for the Stillinger-Lovett second moment
condition [10] to be satisfied.
629
F.O.Raineri, G.Stell
3. Calculations
3.1. Closure relations
The development in section 2 motivated specific forms for the sum bS(r) and
difference bD(r) bridge functions [equations (21) and (27) for bS(r); equation (32)
for bD(r)]. Having this information, we may think of several approximate implemen-
tations of equations (7) and (8):
HNC closure. In this closure we have bS(r) = bD(r) = 0 at every r. The sum and
difference direct correlation functions cS(r) and cD(r) are given by relations
similar to equations (7) and (8), but with bS(r) absent from the argument of
the exponential, and with D0(r) in place of D(r) [cf. equations (13) and (6)].
HNC/PY closure. In this closure λ(r) = bD(r) = 0 at every r. Hence cS(r) and
cD(r) are given by expressions similar to equations (7) and (8), but with bS(r)
given by equation (18) and, again, with D0(r) in place of D(r).
HNC/λ1 closure. In this closure bS(r) in the argument of the exponentials in equa-
tions (7) and (8) is given by equation (21), with λ(r) represented by (27). The
auxiliary function D(r) again takes the form of D0(r) [i. e., in this closure we
again assume that bD(r) = 0 at every r].
bD/λ1 closure. In this closure cS(r) and cD(r) take the form of equations (7) and (8),
with bS(r) given by equation (21), λ(r) given by equation (27), and bD(r) given
by equation (32).
It is obvious that the implementation of closures HNC/λ1 and bD/λ1 require the
specification of the bridge parameters BS and BD. We address this important issue
in the case of the HNC/λ1 closure in the next subsection.
In what remains of section 3 we report numerical results for the structure and
thermodynamic properties of dilute solutions of a model 2–2 electrolyte as calculated
with the HNC, HNC/PY, and HNC/λ1 closures. We briefly comment on closure
bD/λ1 in section 4.
3.2. Numerical algorithm, electrolyte model, and calculations
To solve the Ornstein-Zernike integral equations [cf. equations (9) and (10)]
under any of the above closure relations we have used the familiar Picard iterative
approach. The technical difficulties associated with the long-range nature of the
difference direct correlation function cD(r) are avoided by means of Ng’s method
[11,12].
As a representative of a restricted model of an electrolyte solution, we consider
the model of a 2–2 electrolyte studied by Duh and Haymet [6]. In this model the
ions interact through a McMillan-Mayer level potential of the form of equation (2),
with ū∗(r) given by
ū∗(r) =
kBB
σ
( σ
r
)9
, (36)
630
Bridge functions for electrolyte solutions
where B = 5377(ze)2 Å K and σ = 2.8428 Å. In the Coulombic term ūD(r) we specify
z = 2 and ǫ0 = 78.358, the value of the dielectric constant of water at 298.15K.
The calculations reported here for this model electrolyte parallel those in the
study by Duh and Haymet [6]. We have solved each of the integral equation closures
summarized in the previous section for electrolyte solutions at several concentra-
tions in the range 0.001mol dm−3 6 cs 6 0.5625mol dm−3, where cs is the molar
concentration of the electrolyte.
In addition to the structure functions [the cation-cation g11(r) and the cation-
anion g12 radial distribution functions, obtained from gS(r) and hD(r) by inversion
of equations (3)], we have calculated the following thermodynamic properties:
(a) the reduced excess internal energy per ion u ≡ βU ex/Nt:
u = 2πβρt
{∫
∞
0
dr r2 ū∗(r) gS(r) +
z2e2
ǫ0
∫
∞
0
dr r hD(r)
}
. (37)
(b) the osmotic coefficient φ ≡ βπ/ρt of the solution:
φ = 1−
2πβρt
3
{∫
∞
0
dr r3
dū∗(r)
dr
gS(r) −
z2e2
ǫ0
∫
∞
0
dr r hD(r)
}
. (38)
(c) the generalized compressibility:
β
(
∂π
∂ρt
)
T
= 1 − 4πρt
∫
∞
0
dr r2 cS(r) . (39)
In the above equations π is the osmotic pressure of the solution. It may be shown
that the generalized compressibility equals
β
(
∂π
∂ρt
)
T
= βρs
(
∂µs
∂ρt
)
T
= 1 +
(
∂ ln γ±
∂ lnρs
)
T
, (40)
where µs and γ± are, respectively, the chemical potential and the mean activity
coefficient of the electrolyte, while ρs is the number density of the salt (ρs = ρt/ν,
where ν = ν1 + ν2 is the sum of the cation and anion stoichiometric coefficients).
It is important to remember that equations (37)–(39) correspond to different
thermodynamic routes [energy route for u, virial route for φ, and compressibility
route for β(∂π/∂ρt)T ], and when implemented using correlation functions gS(r) and
hD(r) calculated from approximate integral equation closures, their results are not
necessarily thermodinamically consistent.
3.3. Results
As a representative of the results obtained at each concentration examined, we
first discuss in some detail the results obtained under the HNC, HNC/PY, and
HNC/λ1 closures for the model electrolyte solution at concentration cs = 0.02mol
dm−3 at temperature 298.15K.
631
F.O.Raineri, G.Stell
0 5 10 15
r
0.0
0.5
1.0
g 1
1(
r)
(a)
HNC
HNC/PY
0 5 10 15
r
0
10
20
30
40
g 1
2(
r)
(b)
HNC
HNC/PY
Figure 1. Ion-ion radial distribution functions for an aqueous solution of a re-
stricted 2–2 electrolyte at cs = 0.02mol dm−3 and 298.15K. The correlation func-
tion are calculated under the integral equation closures HNC and HNC/PY. (a):
Like-sign (cation-cation or anion-anion) correlation function g11(r); (b): cation-
anion correlation function g12(r). r is in units of Å.
In figure 1 we report the pair correlation functions g11(r) and g12(r) at this
concentration as calculated under the HNC and the HNC/PY integral equation
closures. The first thing to notice is that (gmax
12 )MD = 48.0, the peak height for the
cation-anion correlation function g12(r) calculated frommolecular dynamics comput-
er simulation [5], is substantially larger than the corresponding HNC and HNC/PY
results, (gmax
12 )HNC = 37.6 and (gmax
12 )HNC/PY = 33.7. From this perspective, and in
spite of the fact that for this closure bS(r) 6= 0, the HNC/PY closure does not im-
prove over the HNC predictions. We notice, however, that there is no evidence of
the spurious peak at r ≃ 8.5 Å in the like-sign ion-ion correlation function g11(r)
calculated under the HNC/PY closure, contrary to the case for g11(r) under HNC
closure. The simulation study [5] indicates that such a peak is an artifact of the
approximate HNC closure.
With respect to the thermodynamic properties at this concentration, we find
uHNC/PY = −1.245 and φHNC/PY = 0.709, which should be compared with the HNC
results uHNC = −1.280 and φHNC = 0.713 and with the molecular dynamics results
[5,6] uMD = −1.413 and φMD = 0.706.
Thus, except for the case of the osmotic coefficient and the absence of the peak in
g11(r), at this concentration of the 2–2 electrolyte the structure and thermodynamic
properties obtained with the HNC/PY closure differ from the simulation results
by a larger extent than the predictions of the HNC closure. This pattern in the
performance comparison between HNC and HNC/PY is also observed, in different
632
Bridge functions for electrolyte solutions
0.0 0.2 0.4
B
S
-1.40
-1.35
-1.30
-1.25
-1.20
u
(a)
0.0 0.2 0.4
B
S
0.700
0.705
0.710
0.715
0.720
φ
(b)
0.0 0.2 0.4
B
S
30
40
50
60
70
80
P
ea
k
H
ei
gh
t
(c)
Figure 2. Dependence of thermodynamic and structural properties of an elec-
trolyte solution with the bridge parameter BS under closure HNC/λ1. The re-
sults shown correspond to a 2–2 aqueous electrolyte solution of concentration
cs = 0.02mol dm−3 at 298.15K. (a) Reduced excess internal energy per ion,
u, calculated with equation (37); (b) osmotic coefficient φ, calculated with equa-
tion (38); (c) peak height gmax
12 in the cation-anion pair correlation function g12(r).
degrees, at the other electrolyte concentrations. The differences between the two
closures disappear for concentrations above cs ≃ 0.5mol dm−3.
We now turn to the HNC/λ1 integral equation closure. The HNC/PY closure
discussed above is a particular case of HNC/λ1, and is obtained when we set to
zero the bridge parameter BS in equation (27). The complete specification of the
HNC/λ1 closure requires that we select the value of BS.
In figure 2 we illustrate the dependence of u, φ, and gmax
12 with BS at cs = 0.02mol
dm−3 and 298.15K. We observe that both u and gmax
12 are quite sensitive to the
value of BS, and that as we increase the value of this parameter both u and gmax
12
change, respectively, to more negative and more positive values, thus improving the
predictions of the closure when compared with the molecular dynamics data. On the
other hand, φHNC/PY is already in good agreement with φMD, and as BS is increased
the agreement of φ under this closure deteriorates. It should be noted, however,
that the dependence of φ with BS is much milder than the dependence of the other
properties.
More careful examination of the figure reveals u and gmax
12 change withBS towards
the target values (respectively uMD and (gmax
12 )MD) at different rates. For example,
at BS = 0.14 we find that (gmax
12 )HNC/λ1
= 47.49, which approximately matches the
molecular dynamics result (gmax
12 )MD = 48.0. However, at this particular value of BS
we find uHNC/λ1
= −1.302 which, although an improvement over both the HNC/PY
633
F.O.Raineri, G.Stell
0.0 0.2 0.4
B
S
0.60
0.65
0.70
0.75
0.80
G
en
er
al
iz
ed
C
om
pr
es
si
bi
lit
y
(a)
0.0 0.2 0.4
B
S
-0.8
-0.6
-0.4
-0.2
0.0
T
er
m
s
in
e
qu
at
io
n
(4
0)
(b)
Figure 3. Sensitivity of two thermodynamic consistency tests with the bridge
parameter BS in the HNC/λ1 integral equation closure. The results shown cor-
respond to a 2–2 aqueous electrolyte solution of concentration cs = 0.02mol
dm−3 at 298.15K. (a) Generalized compressibility test: β(∂π/∂ρt)T calculated
according to the compressibility route (upper curve) and the virial route (lower
curve). (b) Maxwel relation derived test: (energy route)-based left hand side of
equation (42) (lower curve) and (compressibility route)-based right hand side of
equation (42) (upper curve).
and HNC values, respectively uHNC/PY = −1.245 and uHNC = −1.280, differs from
the considerably more negative molecular dynamics result uMD = −1.413. With
regard to the osmotic coefficient, we find φHNC/λ1
= 0.7105 at BS = 0.14, which
should be compared with the simulation result φMD = 0.706. The HNC/λ1 result
is slightly better than the HNC prediction, φHNC = 0.713, but as mentioned above,
φHNC/PY = 0.709 is in a better agreement with simulation. Thus, at BS = 0.14,
when (gmax
12 )HNC/λ1
≃ (gmax
12 )MD, the reduced excess energy per ion uHNC/λ1
betters
uHNC, but it is still far from the target uMD.
Conversely, we observe from figure 2a that at BS = 0.4 the reduced excess
energy is uHNC/λ1
= −1.372, which is still short of matching uMD = −1.413. At this
value of BS, however, (g
max
12 )HNC/λ1
= 76 is far too large when compared with the
(gmax
12 )MD = 48.0.
To further illustrate the dependence with BS of the thermodynamic properties
calculated under the HNC/λ1 closure, figure 3 reports two independent consistency
tests for properties calculated according to different thermodynamic routes.
One test [figure 3a] addresses the thermodynamic consistency for the generalized
compressibility β(∂π/∂ρt)T under the HNC/λ1 closure when this property is calcu-
lated using either the compressibility route [equation (39)] or the virial route. The
634
Bridge functions for electrolyte solutions
latter route involves numerical computation of the indicated partial derivative, with
the function βπ at each density calculated with the relation βπ = ρtφ, in which the
osmotic coefficient is obtained with the help of the (virial route)-based equation (38).
The other thermodynamic consistency test considered follows from the Maxwell
relation between the reduced internal energy and the reduced chemical potential of
the electrolyte [8]:
ν
(
∂(ρt[3/2 + u])
∂ρt
)
β
=
(
∂(βµs)
∂β
)
ρt
. (41)
Notice that the left hand side of the relation requires the full reduced internal
energy per ion; the term 3/2 is just the value of the non-interacting contribution
βU id/Nt to the reduced internal energy. By taking the derivative of each side of
equation (41) with respect to the total ionic density ρt we derive, after taking into
account the first equality of equation (40),
ρ2t
∂2u
∂ρ2t
+ 2ρt
∂u
∂ρt
= −T
∂
∂T
{
β
(
∂π
∂ρt
)
T
}
. (42)
It should be noted that the partial derivatives on the left hand side of this equation
are carried at constant temperature, while those on the right hand side are carried at
constant density of the ions. To use this equation as a measure of the thermodynamic
consistency of closure HNC/λ1, we calculate each side of the equation following a
different thermodynamic route. Specifically, for the left hand side we calculate u at
several values of ρt using the energy route, equation (37). For the right hand side
we calculate β(∂π/∂ρt)T at several temperatures using the compressibility route,
equation (39). The results are shown in figure 3b.
Figure 3a shows that at BS = 0, which corresponds to the HNC/PY closure,
the difference between the values of the generalized compressibility β(∂π/∂ρ t)T cal-
culated under the compressibility and virial routes is the largest, with the value
calculated by the compressibility route larger than the one obtained through the
virial route. The thermodynamic inconsistency decreases, however, when BS is in-
creased, and the figure shows that at BS ≃ 0.395 the two routes in fact give the same
value for β(∂π/∂ρt)T . Although this is satisfying, this thermodynamic consistency
is achieved at a value of BS where gmax
12 is too large and, furthermore, for which the
shape of g11(r) differs considerably from the like-sign correlation function extracted
from molecular dynamics simulation (see below).
Figure 3b illustrates the performance of the HNC/λ1 closure with respect to the
second thermodynamic consistency test, equation (42). It is evident from the fig-
ure that under this closure the (energy route)-based left hand side of the relation is
different at every BS from the (compressibility route)-based right hand side. Further-
more, the magnitude of the thermodynamic inconsistency is practically independent
of the value of the bridge parameter BS.
Finally, in figure 4 we report the ion-ion correlation functions g11(r) and g12(r) for
the 2–2 electrolyte at concentration cs = 0.02mol dm−3. In addition to the HNC and
HNC/PY results already shown in figure 1, we also display the correlation functions
635
F.O.Raineri, G.Stell
0 5 10 15
r
0.0
0.5
1.0
g 1
1(
r)
(a)
HNC/λ1 0.14
HNC/λ1 0.083
HNC
HNC/PY
0 5 10 15
r
0
10
20
30
40
50
g 1
2(
r) (b) HNC/λ1 0.14
HNC/λ1 0.083
HNC
HNC/PY
Figure 4. Ion-ion radial distribution functions for an aqueous solution of a re-
stricted 2–2 electrolyte at cs = 0.02mol dm−3 and 298.15K. The correlation
function are calculated under the integral equation closures HNC, HNC/PY, and
HNC/λ1 for two values of the bridge parameter BS . These are BS = 0.14, when
(gmax
12 )HNC/λ1
= (gmax
12 )HNC, and BS = 0.083, when uHNC/λ1
= uHNC. (a): Like-
sign (cation-cation or anion-anion) correlation function g11(r); (b): cation-anion
correlation function g12(r). r is in units of Å.
calculated under HNC/λ1 using two different choices, BS = 0.14 and BS = 0.083,
of the bridge parameter.
As discussed above, when BS = 0.14 (at cs = 0.02mol dm−3) the HNC/λ1 clo-
sure gives a cation-anion correlation function g12(r) whose preak height (gmax
12 )HNC/λ1
approximately equals (gmax
12 )MD, the value observed in the molecular dynamics sim-
ulation. Figure 4b shows that for this value of BS the HNC/λ1 correlation function
g12(r) compares very favourably with the corresponding molecular dynamics result
[5,6]. (We refer the reader to references [5,6] for the relevant figures of the corre-
lation functions obtained by molecular dynamics). Unfortunately, figure 4a reveals
that for this value of BS the like-sign correlation function g11(r) develops the same
spurious peak at r ≃ 8.5 Å that plagues the HNC closure. We have found this same
behaviour of closure HNC/λ1 at electrolyte concentrations cs 6 0.02mol dm−3: for
values of BS such that (gmax
12 )HNC/λ1
≃ (gmax
12 )MD the shape of g11(r) under HNC/λ1
deteriorates to a level comparable to g11(r) under the HNC closure. For even larger
values of BS [as would be needed for achieving compressibility-virial thermodynamic
consistency in the computation of β(∂π/∂ρt)T ] the spurious peak is more intense,
and the quality of the predicted structure of the solution is not acceptable.
The results reported so far eliminate the possibility of seeking thermodynam-
ic consistency as a practical and unbiased approach for fixing the strength of the
636
Bridge functions for electrolyte solutions
0 5 10 15
r
0.00
0.10
0.20
0.30
0.40
0.50
g 1
1(
r)
(a)
HNC/λ1 0.047
HNC
HNC/PY
0 5 10 15
r
0
50
100
150
200
250
g 1
2(
r)
(b)
HNC/λ1 0.047
HNC
HNC/PY
Figure 5. Ion-ion radial distribution functions for an aqueous solution of a re-
stricted 2–2 electrolyte at cs = 0.001mol dm−3 and 298.15K. The correlation
function are calculated under the integral equation closures HNC, HNC/PY, and
HNC/λ1 with the bridge parameter BS = 0.047, for which uHNC/λ1
= uHNC. (a):
Like-sign (cation-cation or anion-anion) correlation function g11(r); (b): cation-
anion correlation function g12(r). r is in units of Å.
bridge parameter BS in the HNC/λ1 closure. An appealing alternative for estab-
lishing the value of BS that gives under HNC/λ1 improved results (relative to HNC
and HNC/PY) is suggested by the correlation functions gjl(r) reported in figure 4
using BS = 0.083. For this value of the bridge parameter (at cs = 0.02mol dm−3)
uHNC/λ1
≃ uHNC = −1.280. Inspection of figure 4a reveals that at this value of BS
there is only a minimal trace of the spurious peak in g11(r), with the overall shape of
the function in better agreement with the simulation result [5,6]. At the same time
figure 4b shows that g12(r), although not as good as the result at BS = 0.14, is still
an improvement over the predictions of HNC and HNC/PY [(gmax
12 )HNC/λ1
= 41.75
versus (gmax
12 )HNC = 37.6, (gmax
12 )HNC/PY = 33.7, and (gmax
12 )MD = 48.0]. For the os-
motic coefficient we find φHNC/λ1
= φHNC/PY = 0.709, which is in better agreement
with φMD = 0.706 than the HNC result φHNC = 0.713.
We have tested this strategy for fixing the BS parameter of closure HNC/λ1
for solutions of the 2–2 electrolyte in the concentration range 0.001 − 0.5625mol
dm−3. The results obtained are summarized in table 1. The second column in the
table reports the value of BS for which uHNC/λ1
= uHNC at each concentration. It
is interesting that this “optimun” value of the bridge parameter does not change
monotonically with the electrolyte concentration.
By following this procedure, the HNC/λ1 results for the reduced excess energy
u will be obviously equivalent to those of HNC; they are systematically smaller
637
F.O.Raineri, G.Stell
Table 1. Comparison of the predictions of the HNC/λ1 against the results by
molecular dynamics, HNC, and HNC/PY calculations for model aqueous solu-
tions of a 2–2 electrolyte at 298.15K. (a) cs is in units of mol dm−3. (b) For
each electrolyte concentration BS is selected such that uHNC/λ1
≃ uHNC. (c) At
cs = 0.5625mol dm−3, (gmax
12 )MD and φMD have not been reported; the data re-
ported in the table are integral equation results from Duh and Haymet [6] based
on bridge functions extracted from simulation data.
cs
a BS
b uMD
c uHNC/λ1
uHNC uHNC/PY
0.001 0.047 –0.469 –0.432 –0.4328 –0.4053
0.005 0.096 –0.9226 –0.8283 –0.8285 –0.784
0.020 0.083 –1.413 –1.280 –1.280 –1.245
0.0625 0.040 –1.834 –1.713 –1.713 –1.697
0.200 0.005 –2.255 –2.197 –2.197 –2.196
0.5625 0.0 –2.666 –2.648 –2.648 –2.648
cas BS
b φMD
c φHNC/λ1
φHNC φHNC/PY
0.001 0.047 0.898 0.889 0.890 0.892
0.005 0.096 0.808 0.800 0.802 0.802
0.020 0.083 0.706 0.709 0.713 0.709
0.0625 0.040 0.646 0.639 0.642 0.637
0.200 0.005 0.611 0.593 0.594 0.592
0.5625 0.0 0.638 0.615 0.615 0.615
cs
a BS
b (gmax
12 )MD
c (gmax
12 )HNC/λ1
(gmax
12 )HNC (gmax
12 )HNC/PY
0.001 0.047 234 231 186 147
0.005 0.096 112 104.7 85 69.6
0.020 0.083 48 41.75 37.6 33.8
0.0625 0.040 22.5 18.94 18.4 17.7
0.200 0.005 9.6 8.82 8.81 8.77
0.5625 0.0 5.02 4.73 4.73 4.73
638
Bridge functions for electrolyte solutions
in magnitude than the corresponding molecular dynamics results. For the osmotic
coefficient, HNC and HNC/PY are quite comparable, and the table shows that
these two closures perform marginally better against simulation than the HNC/λ 1
closure. On the other hand, the last part of the table reveals that HNC/λ1 performs
considerably better than HNC and HNC/PY with respect to the height gmax
12 of the
peak in the cation-anion radial distribution function. Especially encouraging are the
results at the smallest concentrations cs = 0.001 and 0.005mol dm−3. To illustrate
this point, in figure 5 we compare the gjl(r) at cs = 0.001mol dm−3 calculated under
closures HNC, HNC/PY, and HNC/λ1 with BS = 0.047. From table 1 we see that
(gmax
12 )HNC/λ1
= 231 compares very well with (gmax
12 )MD = 234. Inspection of figure 5b
reveals that the peak in g12(r) predicted by HNC/λ1 is slightly narrower than the
peak calculated under HNC (and it is also narrower than the peak extracted from
the molecular dynamics data [6]). With respect to the like-sign correlation function
g11(r), the HNC/λ1 closure gives, compared against the simulation results, better
results than HNC/PY. Unfortunately, for BS = 0.047, at which uHNC/λ1
= uHNC,
g11(r) under closure HNC/λ1 displays the spurious peak closeby to r ≃ 8 Å. The
size of the peak, however, is much smaller than the peak in g11(r) under the HNC
closure.
4. Discussion
Starting from the previously reported HNC/PY integral equation closure for
solutions of binary electrolytes [7], we have proposed a very simple strategy for
improving the closure performance concerning the structure and thermodynamic
properties of the solution. Basically, the improvement is accomplished by adding new
terms to the ion-ion bridge functions bjl(r) of the HNC/PY closure. The general idea
is to borrow from the first few terms in the asymptotic expansion of the HNC/PY
cjl(r) the specific functional form of the new terms in the bridge functions bjl(r).
We note, however, that although the original HNC/PY work was concerned with
general McMillan-Mayer level models of binary electrolytes, the present extension of
the HNC/PY closure has been implemented in the present work only for restricted
(or symmetric) models of electrolytes. In fact, by exploiting the symmetry of such
electrolyte models, the discussion was actually cast in the language of the sum and
difference functions introduced in equation (3).
The calculations presented in section 3.3 tested the predictions of the simplest
extension, the HNC/λ1 closure, for aqueous solutions of a restricted 2–2 electrolyte.
For this closure, only the sum combination bS(r) of ion-ion bridge functions is mod-
ified relative to the corresponding HNC/PY bridge function; the difference bridge
function bD(r) for this closure, like in the HNC and HNC/PY closures, is assumed
to vanish at every r. The explicit of expression for bS(r) is given by equations (21)
and (27), and contains an unspecified parameter BS. Initially we had hoped that the
optimal value of BS would be fixed by demanding thermodynamic consistency for
β(∂π/∂ρt)T when calculated through the compressibility and virial routes. While
at every concentration studied we succeeded finding a BS for which the thermody-
639
F.O.Raineri, G.Stell
namic consistency is realized, the calculated correlation functions g jl(r) at those BS
deviate substantially from the structure functions gjl(r) extracted from computer
simulations. We were therefore forced to adopt the less elegant (but quite practical)
criterion of choosing BS such that uHNC/λ1
= uHNC is satisfied.
The performance of the HNC/λ1 closure based on this criterion for fixing BS was
examined in section 3.3, where we concluded that, except for the osmotic coefficient,
this new integral equation provides a substantial improvement over the HNC/PY
and the HNC closures. Here we notice that, while still insisting on the approximation
bD(r) = 0, there is room for improving the closure relation by further extending bS(r)
with a 2-parameter λ(r) function:
λ(r) =
BS
2
hD(r)χ(r) +
B′
S
8
hD(r) [χ(r) ]
3 . (43)
The functional form of each term is borrowed from the respective term in the
asymptotic expansion of cS(r) under the HNC/PY closure [the second term in equa-
tion (26)]. We refer to this closure by the acronym HNC/λ2, as its implementation
requires the specification of the two bridge parameters BS and B′
S.
When considered from the viewpoint of the ion-ion bridge functions bjl(r), the
assumption bD(r) = 0 in the HNC/λn family of closures implies that under these clo-
sure relations b11(r) = b12(r) applies. Ichiye and Haymet [13] and Duh and Haymet
[6] have extracted the bridge functions for 2–2 electrolytes at several concentrations
from computer simulations, and their results indicate that b11(r) and b12(r) have
different sign, with b12(r) > 0. (See also Kjellander and Mitchell [14]). The bridge
functions in the Ionic-Percus Yevick (IPY) integral equation closures developed by
Ichiye and Haymet [6,13,15] specifically account for this feature, and are able to
account for the simulation results to a better degree than the HNC/λ1 closure de-
scribed in section 3.3. It should be mentioned that the IPY closures have no free
parameter. To our knowledge, the behaviour of these closures in the context of ther-
modynamic consistency between the various routes has not been examined. (The
exponential exp[b12(r)] in the IPY2 closure [15] contains a term with the factor
exp[−τ12(r)], where τ12(r) is an optimally designed t12(r) function. The presence of
this exp[−τ12(r)] factor suggests that the range of the corresponding direct correla-
tion function would become infinite as the electrolyte solution approaches the limit
of material stability).
Only when we allow bD(r) 6= 0 we can have b11(r) and b12(r) of different sign.
This leads us to the bD/λ1 integral equation closure introduced in section 3.1. This
closure requires the specification of the two bridge strength parameters BS and BD.
The criterions for fixing these parameters and the performance of this new integral
equation are the subject of ongoing research in our group.
We conclude with the mention that, for the dilute systems studied here, the
results of the HNC/PY closure are practically identical to the results calculated
with the Livshits-Martinov integral equation relation [16].
640
Bridge functions for electrolyte solutions
Acknowledgements
We dedicate this work to Jean-Pierre Badiali on the occasion of his 60th an-
niversary. 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.
References
1. McMillan W.G., Mayer J.E. The statistical thermodynamics of multicomponent sys-
tems. // J. Chem. Phys., 1945, vol. 13, No. 7, p. 276–305.
2. Friedman H.L. A Course in Statistical Mechanics. Englewood Cliffs, New Jersey,
Prentice-Hall, 1985.
3. Hansen J.-P., McDonald, I.R. Theory of Simple Liquids (2nd. edition). London, Aca-
demic Press, 1990.
4. Rossky P.J., Dudowicz J.B., Tembe B.L., Friedman H.L. Ionic association in model
2–2 electrolyte solutions. // J. Chem. Phys., 1980, vol. 73, No. 7, p. 3372–3383.
5. Smith D., Kalyuzhnyi Yu.V., Haymet A.D.J. Computer simulation of a model 2–2
electrolyte: multiple time-step molecular dynamics. // J. Chem. Phys., 1991, vol. 95,
No. 12, p. 9165–9171.
6. Duh D.-M., Haymet A.D.J. Integral equation theory for charged liquids: model 2–2
electrolytes and the bridge function. // J. Chem. Phys., 1992, vol. 97, No. 10, p. 7716–
7729.
7. Raineri F.O., Friedman H.L., Stell G. Single-ion contributions to activity coefficient
derivatives, second moment coefficients, and the liquid junction potential. // J. Solu-
tion Chem., 1999, vol. 28, No. 5, p. 463–488.
8. Xu H., Friedman H.L., Raineri F.O. Electrolyte solutions that unmix. Hydrophobic
ions in water. // J. Solution Chem., 1991, vol 20, No. 8, p. 739–773.
9. Stell G. Criticality and phase transitions in ionic liquids. // J. Stat. Phys., 1995,
vol 78, Nos. 1/2, p. 197–238.
10. Stillinger F., Lovett R. General restriction on the distribution of ions in electrolytes.
// J. Chem. Phys., 1968, vol. 49, No. 5, p. 1991–1994.
11. Ng K.-C. Hypernetted chain solutions for the classical one-component plasma up to
Γ=7000. // J. Chem. Phys., 1974, vol. 61, No. 7, p. 2680–2689.
12. Morris G.P., MacGowan D. Solution of the SSOZ equation for molecules of arbitrary
symmetry. // Mol. Phys., 1986, vol. 58, No. 4, p. 745–761.
13. Ichiye T., Haymet A.J.D. Integral equation theory of ionic solutions. // J. Chem.
Phys. 1990, vol. 93, No. 12, p. 8954–8962.
14. Kjellander R., Mitchell D.J. Dressed-ion theory for electrolyte solutions: a Debye-
Hückel-like reformulation of the exact theory for the primitive model. // J. Chem.
Phys., 1994, vol. 101, No. 1, p. 603–626.
15. Babu C.S., Ichiye T. New integral equation theory for primitive model ionic liquids:
from electrolytes to molten salts. // J. Chem. Phys., 1994, vol. 100, No. 12, p. 9147–
9155.
16. Livshitzs A.I., Martinov G.A. The statistical theory of electrolyte melts. // J. Chem.
Phys., 1989, vol. 90, No. 11, p. 6603–6609.
641
F.O.Raineri, G.Stell
Новий клас місткових функцій для розчинів
електролітів
Ф.О.Райнері, Дж.Стелл
Державний університет Нью-Йорка, Стоні Брук, Нью-Йорк, США
Отримано 5 листопада 2001 р.
Для спрощених моделей розчинів електролітів ми представляємо
новий набір замикань на рівні МакМіллана–Маєра, який покращує
передбачення для парних кореляційних функцій іон-іон, отриманих
у гіперланцюговому наближенні. Покращення здійснюється введен-
ням простих функціональних форм для місткових функцій і деяких
підгоночних параметрів, які відповідають різним критеріям. При но-
вих умовах замикання, на противагу до гіперланцюгового наближен-
ня, пряма кореляційна функція “сума”, що є важливою для визначен-
ня стійкості розв’язку по відношенню до фазового розділення, зали-
шається скінченою при термодинамічних станах уздовж спінодалі і в
критичній точці.
Ключові слова: розчини електролітів, спрощена модель,
структура, термодинаміка, місткові функції, термодинамічна
сумісність
PACS: 61.20.Gy, 61.20.Qg, 64.60.Fr, 64.70.Ja, 64.75.+g
642
|