Local density of overlapping spheres near a hard wall: A density functional approach

25 December 1995

PHYSICS

LETTERS

A

Physics Letters A 209 ( 1995) 3 17-320

Local density of overlapping spheres near a hard wall: A density functional approach Andrij Trokhymchuk a lnstituto b DepArtAmenfO

de Fisica.

Universidad

a*‘, Douglas Henderson b, Stefan Sokolowski b,2 de @imica,

Authoma

UNAM,

Coyoach,

Metropa~ilanall~lapalapa,

DF 04510, Mexico Apdu. Postal

55-534,

1:tapalapa.

DF 09340.

Mexico

Received 25 May 1995; accepted for publication 19 October I995 Communicated

by A.R. Bishop

Abstract

A nonlocal density functional theory, developed by Zhou and Stell [J. Chem. Phys. 92 ( 1990) 55331 for simple nonuniform fluids is used to describe the local density of overlapping particles near a hard wall. A comparison of the results of this theory and those obtained from the singlet level approach with computer simulations is presented. This comparison shows that the density functional theory gives quite satisfactory results.

The structure

of simple nonuniform fluids can be by using a singlet level approach [ I] or a pair level theory [ 2-41. It is known, however, that the singlet integral equations, such as the PercusYevick equation (PY 1) or the hypernetted chain equation (HNC 1) are not able to account for surface phase transitions, such as wetting or layering transitions [ 51. The pair theories are able to describe wetting transitions [ 61, but their solution requires remarkable numerical work. As a result, the interpretation of wetting phenomena has relied mainly on density functional (DF) theories. The latter have been proposed in several versions and each of them provides a recipe for calculating the properties of inhomogeneous fluids (see Ref. [ 51 for a review). The commonly used approaches are based on a perturbational division of determined

’ Permanent

either

address: Institute for Physics of Condensed Mat-

ter, National Academy of Sciences of the Ukraine,

29001

1 Lviv,

Ukraine. 2 Permanent address: Computer Laboratory, Faculty of Chemistry, MCS University, 20031 Lubhn, Poland.

0.375.96Ot/95/$09.50

@

SSDIO375.9601(95)00835-7

the interparticle potential and on a mean-field type approximation of the contribution to the free energy functional arising from the attractive part of the interparticle forces. An alternative DF approach, based on the use of the Omstein-Zemike (OZ) equation for the wall-particle correlations and on a nonlocal DF closure for the direct wall-particle correlation function has been proposed by Zhou and Stell [ 71. In our previous works [ 8,9] we have performed a study of the structure of chemically associating fluids (CAPS) near solid surfaces by using mainly the Cummings-&e11 [lo] model of a bulk fluid and the singlet level approximations for the density profiles. A comparison with Monte Carlo (MC) simulations indicated [9] that, for dimerizing hard spheres near a hard wall, the PYl and HNCl approximations reproduce the MC data reasonably well, except for the region close to the wall. These approximations, however, seem to have a rather limited application in studies of the phase behaviour of nonuniform CAFs. Recently, Kierlik and Rosinberg [ 111 have applied

I995 Elsevier Science B.V. All rights reserved

318

A. Trokhymchuk

et al. /Physics

the DF theory to describe the structure of a nonuniform CAF, employing the Wertheim model [ 131. Their theory, however, is suitable for studying a system with adhesive associative forces. Looking for other possibilities to describe the density distribution of overlapping spheres, we have decided to test the DF theory [7] and this is the aim of this work. Theoretical predictions are compared with the results of MC simulation. The model we investigate here consists of hard spheres of diameter u = 1 with an intracore squarewell attraction, u(r) = co,

r<

L-iw,

= -_E,

L-iw
=C3,

L+iw
= 0,

r>

1,

(1)

where E characterizes the strength of the “associative” interactions, w is the width of the square well and L is the average bond length. Note that this model differs considerably from that used in Ref. [ lo]. In particular, geometrical constraints, which determine the type of products of association in the case of a two-component model, are no longer valid. Indeed, the model allows for the formation of n-mers. This model was considered for the first time by Cummings and Stell [ 131; in their work the reader can find a discussion on the allowed classes of n-mers formed for different values of L. In this work we restrict ourselves to the case L < i. In this case only dimers and n-mers built up of rigid, regular polygons are possible. In the case of the analogous model [ 101, the theoretical description is based on the procedure introduced by Baxter [ 141 (for a discussion of the Baxter model see Ref. [ I51 ) . According to this method a deep and narrow square well is replaced by a stickytype potential, defining a new Mayer function, f’ ( r), containing the a-function, f*(r)

= -1 + &S(r

= 0, The

- L),

r<

1,

r>

1.

(2)

parameter 7 is evaluated by equalizing the second virial coefficients for the initial potential ( 1) and the Mayer function (2), i.e. from 7 = L3exp(-@)/( 12L’w + w”) [ 101.

Letters A 209 (199.5) 317-320

Having defined the Baxter-type replica of the Hamiltonian, the structure and thermodynamic properties of the bulk CAF can be evaluated in a standard manner [ 7,8], using the integral equation approach. In this study, we apply the Percus-Yevick approximation to obtain the bulk correlation functions and the chemical potential. The latter quantity is computed as pa$“/apa = - J dr c(r), where ,u”, is the excess chemical potential over the ideal gas term, po is the bulk density and c(r) is the bulk direct correlation function (DCF). The knowledge of both, c(r) and pUex,is necessary to perform the density profile calculations. We are aware that this treatment of the bulk system introduces some inaccuracies, but this choice is just a compromise between accuracy and efficiency; we intend to return to this problem in future studies. The present, exploratory work is devoted to the study of a nonuniform system and the bulk fluid properties are the input to the theory. At the singlet level, the local density p( z ) can be determined by using the PYI or HNCI approximations. The OZ equation for the fluid-wall correlation functions is [ 1 ] hv(r)

= Gv(z> cc

+ 2Tpo

J

03

dz’ h,(z)

--oo

J

dr rc(r),

(3)

I;--ill

where cw(z ) and h,( z ) are the direct and total wallparticle correlation functions. Supplementing the last equation with the Percus-Yevick closure, c;y’(z)

=J%(z)[h,(z)

+ 11 +h,(z.)

- 1,

(4)

or with the hypemetted chain closure, C,(Z)HNC’ = -&J,(Z)

-In]&(z)

+ 11 + h,(z), (5)

we obtain the PY 1 and HNC 1 approximations, respectively. In the above, yw (z ) = exp[ -&.,( z)] and uw(z ) is the potential energy of the particle-wall interaction. The density profile is p( z ) = PO[ h,( z ) + 11. The DF theory [ 71 is based on the expansion of the corrections to the wall-particle DCF, c,(z), of the HNCl approximation. It has been proved that

A. Trokhymchuk et al. / Physics Letters A 209 (1995) 317-320

(6) i=l

The definition of the Sync’ terms can be found in Ref. [7]. Up to now, neither the averaged density, p( z ), nor the total correlation function, h,(z), have been defined. The Zhou and Stell approach [7] consists of the following assumption,

P(a) =

s

dr2 dzdW(lrl - r21),

(7)

where W(r) is a weight function; fi( z ) = po [ h,( z ) + 11. One of the approximations is W(r) = f(r)/ S dr f(r), where f(r) is the Mayer function. This choice is suitable for a repulsive interaction, not for a fluid with attractive forces, such as a CAF or even a Lennard-Jones fluid [ 71, because it can generate negative averaged densities, p( z ) . Thus, we have assumed that f(r) is just a Mayer function for a hard sphere interaction, i.e. f(r) = -1 for r 6 1 and f(r) = 0 otherwise. Other choices may improve the results and we intend to pursue this point systematically in future work. In this connection one can try to introduce selfconsistency into the scheme using a weight function such that the thermodynamics of the system matches the chemical potential input in a homogeneous system. When all SYNC’(z ) terms in Eq. (6) are neglected, Eqs. (3)) (5) and (6) form a closed set, which can be solved. This approximation is called [ 71 “the hydrostatic hypemetted chain approximation” (HHNC 1) . Our numerical procedure was identical to that described in the Appendix B of the second paper of Zhou and Stell [ 71. In order to test the accuracy of the HHNC 1 approximation, we have also performed MC simulations, using the same algorithm as in Refs. [ 9,161. We have simulated the fluid interacting via the potential ( 1) in a rectangular box XL x l2 x ZL of size XL = l!L = 10 and ZL = 15. The simulations were carried out for a number of particles N = 905 and also for N = 500. Periodic boundary conditions in the x and y directions have been applied and the system was closed at the

319

top and the bottom by hard reflecting walls. Ensemble averages were accumulated over at least 80 x lo6 configurations after equilibrating from random distributions for 30 x lo6 configurations. Density profiles were evaluated from histograms updated every N configurations with a grid size of 0.01. The acceptance ratio was of the order of 0.25-0.3. no potential models u(r) were studied: model A: L = 0.45, P& = 2.5215 (which corresponds to T = 0.03) and model B: L = 0.45, PE = 1.317 (7 = 0.10) and in both cases w = 0.1. A comparison of the predictions of the PY 1, HNC 1 and HHNC 1 equations with the MC results is given in Fig. 1. Before performing the surface calculations, we have checked if the density integration of the derivative of the chemical potential does not pass through a metastable region. All our calculations have been carried out far away from bulk instabilities. Figs. la and 1b show the results for models A and B, respectively. In each case the lower group of curves corresponds to a lower density, whereas the upper group of curves corresponds to a higher density. Note the presence of cusps at z = L, resulting from “associative” interaction in the pair potential u( r) . In the case of the MC data these cusps are somewhat smeared out. This is connected with the finite width of the square-well attraction, w = 0.1. As we have already noted, our previous calculations have indicated that the PY 1 and HNC 1 equations correctly reproduce the MC density profiles, except for the region close to the wall. Thus, one could expect that the HHNCl equation should improve the accuracy of the predictions just for small values of z. Indeed, an inspection of all curves given in the figure indicates that the DF theory describes better the behaviour of the density profile in the vicinity of the wall and the contact values, p(O), in particular. Obviously, the agreement of the HHNCl predictions is better for a lower density - the results of the MC simulations and HHNC 1 calculations almost completely coincide. In each case, the PY 1 approximation underestimates, whereas the HNCl approximation overestimates the contact value, p( z ). The HHNCl approximation leads to a better coincidence than the two remaining equations. Some discrepancies between the MC simulations and HHNCl profiles can be attributed either to the difference between the Hamiltonian used in the simulation and in the theoretical calculations,

320

A. Trokhymchuk et (11./Physics

Letter.? A 209 (1995) 317-320

the possibility of formation of several types of association products. One can thus expect that for L > i the weight function used may be unappropriated. Obviously, the approach applied here can be extended to other classes of models of CAFs. The application of Eqs. (6) requires the knowledge of the bulk DCF and the chemical potential. Their source, however, can be quite arbitrary, e.g. computer simulations. Thus, we can also consider the use of the Wertheim theory to obtain the bulk properties while solving the PYI and HNCI equations. Indeed, in the absence of an associative fluid-wall interaction, we do not need to use the Wertheim formalism for the corresponding fluid-wall correlation functions. A generalization of this procedure to the case of a two-component model of the Cummings and Stell type [ IO] is also possible. All these problems will be discussed in a future work. We wish to express our gratitude to the CONACYT of MCxico (Grant No. 4186-E9405 and el Fondo para CBtedras Patrimoniales de Excelencia) for financial support of the project.

1.5

1

References 0.5

III D. Henderson, ed., in: Fundamentals fluids (Dekker,

I21 D. 0 Fig.

I.

2

1

Comparison of the MC, HHNCI,

profiles. Results for: (a)

3

,5

I

density

(lower group of

[31 S.

(lower group of curves) and at po = 0.7 I5 (upper

group of curves). are for the HHNCI,

In (a) PYI

approximations; in (b)

curves at ; = 0 from the bottom are for the PY I, HHNCI HNCI

approximations.

the

Henderson (Dekker,

New York, 1992).

(1987)

1210.

[71 Y. Zhou and G. Stell, J. Chem. Phys. 92

( 1990)

5533.5544;

L. Blum and G. Stell, J. Stat. Phys. I5 ( 1976) 439.

and

In the case of lower densities for both

cases the simulation and the HHNCI

Sokolowski, J. Chem. Phys. 73 ( 1980) 3507.

161 E. Bruno, C. Caccamo and P. Tarazona, Phys. Rev. A 35

the curves at ; = 0 from the bottom and HNCl

( 1985)

[41 P Attard, J. Chem. Phys. 91 ( 1989) 3072. I51 R. Evans, in: Fundamentals of inhomogeneous fluids, ed. D.

curves) and at po = 0.741 (upper group of curves); (b) model B at po = 0.409

Henderson and M. Plischke, Proc. R. Sot. A 400

163.

HNC 1 and PY

model A at pt, = 0.403

of inhomogeneous

New York, 1992).

[81 0. Pizio, D. Henderson and S. Sokotowski, J. Phys. Chem. 99

results practically coincide.

( 1995)

2408.

[91 A. Trokhymchuk, D. Henderson and S. Sokolowski, Can. J. Phys., in press.

as well as to the uncertainties in the chemical potential and the bulk DCF. In general, we can state that the application of the DF-type theory to the case of a CAF is encouraging. This theory is able to predict the structure of the fluid considered here reasonably well and, thus, it is possible to apply this kind of approach to further studies on wetability of surfaces by chemically associating fluids. Assessing the validity of our results, one should remember that all calculations have been done for L < f. This condition restricts

[ 101 PT. Cummings and G. Stell, Mol. Phys. 51 (1987)

[Ill

( 1984)

253; 60

1315.

E. Kierlik

and M. Rosinberg, J. Chem. Phys. 100

( 1994)

1716.

1121 MS. Wertheim, J. Stat. Phys. 35

( 1984)

19.35;

459, 477; J. Chem. Phys. 85 ( 1986) 2929. PT. Cummings and G. Stell, Mol. Phys. 55

42

( 1985)

( 1986) 33.

R.J. Baxter, J. Chem. Phys. 49 ( 1968) 2770. G. Stell, J. Stat. Phys. 63 ( 1991 ) 1203. Yu.V.

Kalyuzhnyi,

G.

Stell,

Chapman and M.F. Holovko, 7939.

M.L.

Llano-Restrepo,

J. Chem. Phys. IO1

W.G.

( 1994)