Capillary condensation in pores with rough walls: A density functional approach

Journal of Colloid and Interface Science 313 (2007) 41–52 www.elsevier.com/locate/jcis

Capillary condensation in pores with rough walls: A density functional approach P. Bryk a , W. R˙zysko a , Al. Malijevsky b , S. Sokołowski a,∗ a Department for the Modeling of Physico-Chemical Processes, Maria Curie-Skłodowska University, 20-031 Lublin, Poland b E. Hála Laboratory of Thermodynamics, ICPF, Academy of Sciences, 165 02 Prague 6, Suchdol, Czech Republic

Received 18 January 2007; accepted 27 March 2007 Available online 27 April 2007

Abstract The effect of surface roughness of slit-like pore walls on the capillary condensation of a spherical particles and short chains is studied. The gas molecules interact with the substrate by a Lennard-Jones (9, 3) potential. The rough layer at each pore wall has a variable thickness and density and consists of a disordered quenched matrix of spherical particles. The system is described in the framework of a density functional approach and using computer simulations. The contribution due to attractive van der Waals interactions between adsorbate molecules is described by using first-order mean spherical approximation and mean-field approximation. © 2007 Elsevier Inc. All rights reserved. Keywords: Adsorption; Pore; Capillary condensation; Density functional theory; Quenched-annealed systems

1. Introduction Real adsorbing surfaces are usually rough, so a fluid upon adsorption encounters geometrically and/or energetically heterogeneous walls. Recently, the effects of heterogeneity on adsorption have been considered theoretically and employing simulation techniques for several model systems [1–19]. In particular, it has been shown that the surface roughness can influence their wettability [20–24]. To study this effect some of us proposed a model of a rough wall created by depositing a layer of molecules on an idealized substrate and next imposing a quenching procedure [23,24]. The quenched layer was then exposed to the adsorbing gas. Surface roughness in this model is a combined effect of the preparation procedure, the density of the deposited layer and its thickness. To describe the adsorption of a fluid a modified Fischer–Methfessel approach [25] that is closely related to the so-called “Mark-I” version of the Tarazona’s density functional theory (DFT) [26] was applied. This model was used to study the effect of surface roughness on the wetting transitions [24], to investigate adsorption of fluids in * Corresponding author.

E-mail address: [email protected] (S. Sokołowski). 0021-9797/$ – see front matter © 2007 Elsevier Inc. All rights reserved. doi:10.1016/j.jcis.2007.03.077

pores with rough walls and to examine the behavior of mixtures near semi-permeable membranes [27–29]. Recently, Schmidt and co-workers proposed a DFT to describe quenched-annealed mixtures [30–33]. Their approach is based on the replica trick and allows the treatment of situations where the quenched random matrix as well as the annealed fluid are inhomogeneous. Applications of the theory included investigation of the adsorption properties of hard spheres and model colloid–polymer mixtures in bulk matrices and at matrix surfaces, and the influence of the quenched disorder on phase transitions like fluid demixing, isotropic-nematic ordering and freezing. Particularly rich wetting behavior was found for colloid–polymer mixtures adsorbed against porous wall. For review of all these applications see Ref. [32]. The theory of Schmidt and co-workers treats both the annealed and the quenched species on the level of their one-body density distributions. Hence the complicated external potential that the quenched particles exert on the fluid never enters explicitly into the theoretical framework. This approach constitutes an enormous simplification as far as practical computations are concerned, but also suffers of several limitations, as it was noted by Lafuente and Cuesta [34]. These authors derived an alternative DFT for quenched-annealed mixtures from first principles

42

P. Bryk et al. / Journal of Colloid and Interface Science 313 (2007) 41–52

without resorting to the replica trick. They showed that the disorder-averaged free energy of the fluid is a functional of the average density profile of the fluid as well as the pair correlation of the fluid and matrix particles. However, an implementation of this theory is difficult and from this reason it is preferable to use the approach of Schmidt and co-workers. To employ theoretical methods for calculating adsorption isotherms a model for fluid–fluid and fluid–adsorbent interactions must be formulated first. The molecular parameters for fluid–fluid interaction are usually determined from bulk thermodynamic data, whereas solid–fluid potential parameters are obtained by fitting theoretical results against experimental adsorption isotherms. It appears, however, that in a vast majority of cases [18,19] the predicted isotherms are stepwise indicating the sequential formation of molecular layers, which, however, are not observed experimentally. In the case of some adsorbents, like graphitized carbon, improvement of the predictability of isotherms is achieved by modifying the fluid–fluid interaction parameters [18,19,35]. The altering of their values is justified by the importance of multi-body effects [35–37]. The success is of this procedure is possible with graphitized carbon black because of the fairly homogeneous nature of the surface, however, such an approach does not enjoy a similar success in the description of adsorption isotherm for other adsorbents, e.g., for silica gel [38–40]. One possible reason is energetic heterogeneity of adsorbing surfaces (chemical heterogeneity) surface roughness (geometrical heterogeneity), which is considered at the molecular level and attributed to a molecular feature of amorphous solid [18,19]. In the case of heterogeneous adsorbents the adsorbing potential is a complicated function that varies in three-dimensional space. Apart from the fact that usually only a little is known about the spatial dependence of the adsorbing potential for real substrates, numerical calculations (including simulations) would be extremely time-consuming. One possible route to describe effectively such systems is to replace the energy landscape by the density of “auxiliary matrix particles” landscape and to employ the methods developed for quenched-annealed mixtures [41]. According to this treatment a complicated adsorbing potential does not enter explicitly theoretical equations. The model developed in our previous works [23,24,27] seems to be appropriate to handle the effects of surface heterogeneity. Indeed, a similar model was used recently by Ravikovitch and Neimark [42], who applied DFT of Schmidt [30,32], and Reich and Schmidt [31] to model adsorption on geometrically heterogeneous surfaces and porous solids with rough pore walls. In this work we present DFT of adsorption of chain molecules in pores with rough walls. The theory combines the approach of Schmidt [32] with the theory of Yu and Wu [43] and Cao and Wu [44]. The latter theory can be successfully applied to describe adsorption of such important from experimental point of view [45–49] complex fluids as carbon dioxide, hydrocarbons, polymers and many others [50]. Moreover, this theory is numerically almost as simple as density functional theories of simple fluids [26]. Obviously, the proposed approach can be also used to study adsorption of spherically symmetric (atomic) molecules.

Previous applications of the density functional theory have been primarily focused on various weighted density approximations for the repulsive part of the Helmholtz energy functional and on the mean field approximation (MFA) for an attractive contribution. The major problem of the MFA is that it significantly underpredicts the critical temperature for uniform fluids; hence the corresponding interfacial properties maybe unphysical. Recently Tang [51] proposed a first-order mean-spherical approximation (FMSA) theory for a long ranged attraction. The FMSA theory, with the analytical expression of a radial distribution function or a direct correlation function as input, significantly improves the MFA for the long ranged attraction and is applicable to both bulk and inhomogeneous systems using a single set of molecular parameters. The DFT for simple nonuniform fluids based on the Fundamental Measure Theory [52–54] for hard-sphere free-energy contribution and on the FMSA for the attractive contribution is probably the most accurate approach to the theory of simple confined fluids [55–58]. Therefore the aim of this paper is also to extend the FMSA to the case of nonuniform systems with quenched component. In the case of adsorption of spherical symmetric molecules the theoretical predictions are compared with the results of computer simulations. The phase diagrams for model systems have been evaluated using hyper-parallel tempering technique [59]. 2. Theory The system under consideration is a fluid composed of chain molecules confined in a slit-like pore of the width H . The chains are built of M segments, each of diameter σss . We stress that setting M = 1 in all equations presented below the theory reduces to that for adsorption of spherically symmetric (atomic) adsorbates. The chain connectivity is assured by imposing the bonding potential [43], Vb (R) = M−1 j =1 vb (|rj +1 − rj |), where R ≡ (r1 , r2 , . . . , rM ) denotes a set of coordinates describing the positions of all segments and vb is the bonding potential between two adjacent segments. We assume that the segments are tangentially jointed, so that the binding potential is     M−1  exp −βVb (R) = δ |rj +1 − rj | − σss /4πσss2 .

(1)

j =1

We also take into account attractive van der Waals potential between the segments. This potential is the same for all the segments and is given by    4εss (σss /r)12 − (σss /r)6 , for r  rcut,ss , uss (r) = (2) 0, for r > rcut,ss . In the above rcut,ss is the cut-off distance, and εss is the energy parameter. Each segment j interacts with a single pore wall via Lennard-Jones (9,3) potential   vj (zj ) = εgs (z0 /zj )9 − (z0 /zj )3 , (3) where zj is the distance between the j th segment and the wall. The total external potential, Vext (R), imposed on a chain molecule is the sum of the potentials acting on individual segments,

P. Bryk et al. / Journal of Colloid and Interface Science 313 (2007) 41–52

 due to both pore walls Vext (R) = M j =1 [vj (zj ) + vj (H − zj )]. The parameters εgs , and z0 are identical for all segments. The last assumption can be easily removed. In order to model pores with rough walls, a layer consisting of a quenched matrix is placed at each wall [23,24,27]. We assume that the matrix is built of spherical molecules of diameter σmm . However, the presented theory can be also extended to the case of matrices made of chain particles. The interaction between matrix species and a segment of the chain is given by the potential ums (r), which can be either a hard-sphere potential  ∞, for r  σms , ums (r) = (4) 0, for r > σms , or a truncated Lennard-Jones potential    4εms (σms /r)12 − (σms /r)6 , for r  rcut,ms , ums (r) = (5) 0, for r > rcut,ms . Here σms = (σss + σmm ) and rcut,ms is the cut-off distance for the fluid–matrix interaction. The local density of the matrix is described by the function ρm (z). At the moment we do not specify the form of ρm (z), nor the process leading to the formation of the matrix distribution at the pore walls. Of course, if the matrix inhomogeneity is generated in response to an external potential, Vm (z), that acts on the matrix particles before the quench, then the matrix configuration is drawn from an equilibrium distribution according to the Hamiltonian of the matrix particles, in particular an equilibrium DFT can be applied for this purpose [30–32]. In general, however, in the DFT of quenched-annealed systems the distribution ρm (z) is an input into the theory and the step of matrix development is completely separated from the second step, which relies on the adsorption of a fluid. Following Schmidt [30–33] the replica trick is applied to write down the grand potential functional for the confined fluid       Ω ρm (z), ρ(R) = Fid ρ(R) + Fex ρm (z), ρ(R)     + dRρ(R) Vext ρ(R) − μ , (6) where ρ(R) and μ are the local density of chain particles and their chemical potential, respectively. Moreover, Fid (ρ(R)) is the configurational ideal free energy and the excess free energy functional, Fex (ρm (z), ρ(R)), describes the interparticle interactions of adsorbate particles with adsorbate particles and with the matrix particles. We are aware that this theory has a number of weak points that should be pointed out. Although the replica trick is a widely applied method of statistical physics, it introduces assumptions that are difficult to justify concerning the analytic continuation of the grand potential as a function of the number of replicas, and the replica symmetry or its breaking. Moreover, contrary to DFT of fluids, the formulation of the density functional theory of quenched-annealed systems by Eq. (6) makes difficult to derive the set of Ornstein–Zernike equations for the system from functional relations. For example, at present it is not at all clear what the meaning of the second derivatives of Fex (ρm (z), ρ(R)) is [34]. Nevertheless, we have decided to follow the ideas of Schmidt [30–32] because at the moment this

43

approach is the only one that can be implemented to the considered case in a relatively easy manner. Similarly as in DFT of simple nonuniform fluids we start with writing down the system free energy. The ideal free energy functional is known exactly [43]   βFid ρ(R) = β dRρ(R)Vb (R)     + dRρ(R) ln ρ(R) − 1 . (7) The excess free energy, Fex , can be split into the hard-sphere contribution, Fhs , which results from the hard-sphere repulsion between segments and between segments and matrix particles, the attractive contribution due to segment–segment interaction, Fatt and, eventually, due to segment–matrix interaction, and the contribution due to chain connectivity, Fc . According to the approach of Yu and Wu, which is based on the Fundamental Measure Theory [52–54] and on Wertheim [60,61] theory of association, the contributions Fhs and Fc are functionals of the weighted densities. In order to write down the definitions of the weighted densities we first introduce the average segment density profiles for chain molecules ρs (r) =

M

ρs,j (r) =

j =1

M

dRδ(r − rj )ρ(R),

(8)

j =1

where ρs,j (r) is the local density of the segment j of the chain. The weighted densities are then defined as   nα (r) = nα,s (r) + nα,m (r) = dr ρs (r )wα |r − r |, dss     + dr ρm r wα |r − r |, dmm , (9) where scalar (α = 0, 1, 2, 3) and vector (α = V 1, V 2) weight functions, wα (r, σx ), x = ss, ms are given in Ref. [52]. The quantities dss and dmm are the effective hard-sphere diameters of the reference hard-sphere fluid. Each of the components, Fα , α = hs or c is expressed as a  volume integral Fα = Φα (r) dr. The hard-sphere contribution is evaluated from [52–54]     Φhs (r) = Φ˜ hs r; {nα,s }, {nα,m } − Φ˜ hs r; {nα,s ≡ 0}, {nα,m } . (10)   Φ˜ hs r; {nα,s }, {nα,m } n 1 n 2 − nV 1 · n V 2 1 − n3  3 n3 + (1 − n3 )2 ln(1 − n3 ) + n32 1 − ξ 2 , 36πn23 (1 − n3 )2

= −n0 ln(1 − n3 ) +

(11)

where ξ(r) = |nV 2 (r)|/n2 (r). The contribution Φc is obtained from Wertheim’s first-order perturbation theory [60,61] Φc (r) =

  1−M n0,s ζ ln yhs (dss ) , M

(12)

44

P. Bryk et al. / Journal of Colloid and Interface Science 313 (2007) 41–52

where ζ = 1 − nV 2,s · nV 2,s /(n2,s )2 and yhs is the contact value of the hard-sphere radial distribution function for the particles of the annealed component in the mixture containing matrix and annealed particles yhs (dss ) =

1 n2 dss ζ (n2 dss )2 ζ + + . 1 − n3 4(1 − n3 )2 72(1 − n3 )3

(13)

The simplest treatment of attractive forces follows from the assumption that dss = σss and dmm = σmm and from the application of a MFA. Introducing the Weeks–Chandler–Andersen division [62] of the Lennard-Jones potentials (Eqs. (2) and (5))  −εx , for r  21/6 σx , x = ss, ms, u˜ x (r) = (14) ux (r), for r > 21/6 σx , x = ss, ms, we have Fatt =

1 2



+



  dr dr u˜ ss |r − r | ρs (r)ρs (r ) dr dr u˜ ms

 |r − r | ρs (r)ρm (r ).



(15)

If the matrix–fluid interactions are of hard-sphere type, the second term in Eq. (15) disappears. As we have already noted, the application of FMSA instead of MFA improves substantially accuracy of theoretical predictions, compared with MFA. However, it is possible only when the fluid–matrix interactions are given by the hard-sphere potential. All the details concerning FMSA can be found in [55–58]. In the case of FMSA the effective hard-sphere diameter dss is used. If the density-dependence of dss is neglected, its temperature-dependence can be approximated by the expression [63] σss   dss = dr 1 − exp −βuss (r) .

(16)

0

The diameter dss is used so often that it is very tempting to employ an approximation by a simple analytical expression. One of the most widely exploited expression is that due to Cotterman et al. [64] dss =

1 + 0.2977T ∗ 1 + 0.33163T ∗ + 0.00104771(T ∗ )2

,

(17)

where T ∗ = kT /εss . A collection of other expressions proposed to approximate Eq. (16) can be found in Silva et al. [65] paper. Sensitivity of thermodynamics of bulk fluids to the Barker– Henderson diameter, which covers both its approximation in calculation and improvement in rationality was discussed by Tang [66]. In the case of nonuniform systems, however, a simplifying assumption that dss is just σss does not change the overall picture of surface phase transition occurring in the system. Obviously, dmm ≡ σmm , by definition. The Helmholtz free energy functional due to the long ranged attraction is then given by   1 dr dr c˜ss |r − r | ρs (r)ρs (r ), Fatt = − (18) 2

where css (|r − r |) is the attractive part of the direct correlation function, which can be expressed as [55–58] ⎧ Y ∗ c (T1 , z1 , ds s, r) ⎪ ⎪ ⎨ − cY (T2∗ , z2 , ds s, r), r < dss , css (r) = (19) ⎪ dss < r < σss , ⎪ ⎩ 0, r > σss , uss (r), where z1 = 2.9637/σss , z2 = 14.0167/σss , T1∗ = T ∗ dss /k1 , T2∗ = T ∗ dss /k2 , k1 = k0 exp[(σss − dss )z1 ] k2 = k0 exp[(σss − dss )z2 ] and k0 = 2.1714σss . The first-order Yukawa direct correlation function, cY , in Eq. (19) is given in appendix of Ref. [67]. The equilibrium configuration of fluid molecules satisfies the Euler–Lagrange equation   ρ(R) = exp βμ − βVb (R) − βv(R) − βΛ(R) , (20) where μ is the chemical potential and Λ(R) = δFex /δρ(R), where Fex = Fhs + Fatt + Fc is the effective potential field due to intra- and intermolecular interactions. The excess free energy functional Fex depends on the average segment density ρs , thus

δFex λi (r), = δρ(R) M

i=1

λi (r) =

δFex . δρs (r)

(21)

In this work we consider systems that are inhomogeneous in only one (say z) direction. Then one obtains (cf. Refs. [43, 68,69])   ρs,i (z) = exp(βμ) exp −βλi (z) Gi (z)GM+1−i (z), (22) where the functions Gi (z) are determined from the recurrence relation   θ (σ − |z − z |) Gi−1 (z ), (23) Gi (z) = dz exp −βλi−1 (z) 2σ for i = 2, 3, . . . , M and with G1 (z) ≡ 1. The confined system is in an equilibrium with the bulk fluid. The bulk fluid is characterized by the same value of the chemical potential μ as the confined system and its properties are described by the equations given above assuming that the external field Vext is zero and that the segment densities ρs,j and the total segment density ρs are constant in the entire bulk system and equal to ρbs,j and ρb,s , respectively. Obviously, the matrix concentration is zero in the bulk system. The numerical methods for the minimization of the grand potential are identical with those described in our previous papers [68,69]. In particular, we introduce the appropriate propagator functions, as described in details in our previous works. For the sake of brevity these equations are not presented here, but we refer the reader to original works [43,68,69]. Finally, we introduce the reduced quantities. All the en∗ = ε /ε , ε ∗ = ergy parameters are expressed in εss units, εgs gs ss ms εms /εss . Similarly, the diameter σss is used as unit of length. ∗ = σ /σ , H ∗ = H /σ , z∗ = z/σ and the reduced Then σms ms ss ss ss ∗ (z) are defined as ρ ∗ (z) = ρ (z)σ 3 , ρ ∗ densities ρs,j s,j ss s,j bs,j = ∗ = ρ (σ )3 , etc. ρbs,j (σss )3 , ρbs bs ss

P. Bryk et al. / Journal of Colloid and Interface Science 313 (2007) 41–52

45

3. Results and discussion As we have noted, setting M = 1 in the relevant equations given above, the theory reduces to that for atomic adsorbates. The aim of our model calculations performed for systems involving spherically symmetric molecules is two fold. First, we would like to compare the results of MFA and FMSA with computer simulations. Secondly, we intend to compare two cases: the case of matrix molecules interacting with adsorbate atoms via repulsive (hard-sphere) interactions and via Lennard-Jones (12, 6) potential. The MFA that was used in a vast majority of papers aiming at the application of DFT to experimental data, is perhaps the simplest possible that it is able to capture the behavior of nonuniform systems with repulsive and attractive forces at quantitative level. However, this approximation should not be treated as qualitatively correct. An attractive alternative for MFA provides FMSA. The latter theory, with analytical expression for a direct correlation function as an input, which makes it numerically as simple as MFA, significantly improves the MFA for the long ranged attraction [55–58]. At the beginning of our discussion we compare the density functional predictions with the results of grand canonical Monte Carlo simulations. Because the simulation method is quite standard, we omit all the technical details and only note that the phase diagrams for the confined systems were obtained employing hyper-parallel tempering technique [59]. The surface roughness was modelled assuming that the matrix distribution is  ρm0 , for 0  z  zm and H − zm  z  H, ρm (z) = (24) 0, otherwise, where zm is the matrix thickness at each pore wall. In the calculations below we have assumed that zm = σss . The diameter of matrix molecules was the same as the Lennard-Jones diam∗ = 1. The matrix–fluid interactions eter of fluid particles, σmm were of hard-sphere type and the cut-off of the fluid–fluid po∗ = 2.5. The pore width was kept constant and tential was rcut,ss ∗ equal H = 5 and the fluid–pore wall potential parameters were ∗ = 9 and z∗ = 0.5. εgs 0 Fig. 1 shows the density profiles obtained from computer simulations and from both theories at supercritical temperature, T ∗ = 1.5. For all the systems in Fig. 1 the bulk (reference) sys∗ = 0.6958, which corresponds to the value tem density was ρb,s of the chemical potential from the bulk system Grand Ensemble Monte Carlo simulations equal to μ/kT = −0.5. It is not surprising that the FMSA improves the agreement with simulations in comparison with the data obtained from the mean-field theory. This is particularly well seen for the system without matrix. The improvement is also evident when one compares the dependence of the average adsorbate densities in the pore, H ρ =

ρs (z) dz/H,

Fig. 1. Density profiles from simulations (symbols), from MFA (solid lines) ∗ = 0.6958, and from FMSA (remaining lines) inside the pore of H ∗ = 5 at ρb,s ∗ = 0.5, 0.2959, at T = 1.5. The consecutive curves from bottom are for ρm0 0.09694 and 0 (i.e., for the pore without the quenched particles). The values of all remaining parameters of the model are given in the text.

Fig. 2. Average fluid densities in the pore versus the bulk density. The consecutive curves from bottom are for ρm0 = 0.5, 0.2959, 0.09694 and 0 (the pore without the quenched particles). The results of computer simulations are given as points, theoretical curves resulting from MFA and from FMSA are given as (· · ·) and (—), respectively. Parameters as the profiles in Fig. 1 and the nomenclature is the same as in Fig. 1.

(25)

0

plotted versus the bulk density, see Fig. 2. The presence of quenched molecules causes significant changes in the structure

of adsorbed fluid. Due to pore walls “screening” by the matrix molecules, the height of the local density peak adjacent to the wall is lowered at higher matrix density ρm0 .

46

P. Bryk et al. / Journal of Colloid and Interface Science 313 (2007) 41–52

Fig. 3. Phase diagrams for Lennard-Jones fluid confined in the pore of H ∗ = 5 ∗ = 0.5 and without the matrix particles. with deposited matrix of density ρm0 Points are the results of the Grand Canonical Monte Carlo simulations, (—, · · ·) are the results of FMSA, whereas (- - -, -·-·) the results of MFA.

Fig. 3 compares the phase diagrams in the plane: the average density of the confined fluid, the temperature reduced by the appropriate bulk critical temperature (the bulk critical temperature resulting from Monte Carlo simulations is Tc∗ = 1.1876 [70]). We realize that the FMSA predicts the phase behavior of the confined system without hard-sphere matrix and with the matrix ∗ = 0.5 much better than the MFA does. This of the density ρm0 result confirms previous findings [57,58,67,71,72] and suggests that when a density functional theory is applied to describe experimental data (cf. [38–40,42]), i.e. when the accuracy at a qualitative level is required, one should use rather FMSA than MFA. Obviously, for model calculations, i.e. when accuracy at a quantitative level suffices, the MFA is good enough to get a quantitative insight into the system behavior. The results reported below are entirely based on the MFA. This is also because we would like to compare the results for the system with hard-sphere and with Lennard-Jones fluid–matrix interactions and the application of an analogue of FMSA to the latter system would be difficult. First, we consider adsorption is a wide pore of H ∗ = 30, ∗ = 12 at a low temperature of T ∗ = 0.75, and assume that εgs ∗ and z0 = 0.5. The matrix density at each pore wall is described ∗ = 2. The matrix is composed of hard by Eq. (24), but now zm spheres of diameter σmm = σss and the fluid–matrix interactions are of hard-sphere type. For the pore without matrix the capillary condensation is preceded by a series of three layering transitions at each wall, cf. Fig. 4, where in part a we have plotted the excess adsorption isotherms Γ ∗ = Γ σss2 , H Γ = 0



 ρs (z) − ρb,s dz,

(26)

(a)

(b) Fig. 4. (a) Excess adsorption isotherm inside the pore of the width H ∗ = 30 ∗ = 0.005 (· · ·) and without matrix (solid line), with the matrix of density ρm0 ∗ ρm0 = 0.05 (- - -). Vertical parts of the isotherms correspond to equilibrium transitions and no loops of hysteresis are shown here. (b) (—, - - -, -·-·) Are re∗ = 0.001), the second (ρ ∗ = 0.0023) spectively the profiles after the first (ρb,s b ∗ = 0.0025) layering transition for the pore of H = 30 withand the third (ρb,s out matrix. (-·-·) Is the profile inside the pore with the matrix characterize by ∗ = 0.05 at ρ ∗ = 0.001, i.e. at the same bulk state point at which the profile ρm0 b,s plotted as solid line was calculated. (· · ·) Corresponds to the state point after the ∗ = 0.005. layering transition (ρb∗ = 0.0025) in the pore with the matrix of ρm0 The values of all remaining model parameters are given in the text.

P. Bryk et al. / Journal of Colloid and Interface Science 313 (2007) 41–52

versus the bulk density and in part b—examples of the density profiles. ∗ = A layer of hard-sphere matrix of very low density, ρm0 0.05, at the pore walls suppresses all the layering transitions and only capillary condensation preserves. Even the presence ∗ = 0.005 causes that the first and of the matrix of density ρm0 the third layering transitions vanish. The structure of adsorbed fluid, evaluated at the same bulk density for the pore with and without matrix is different. However, density profiles for different bulk densities, but at a fixed value of Γ , are similar (cf. Fig. 4b). The effect of suppressing layering transitions has been already observed and its occurrence led Ravikovitch and Neimark [42] to the application of a similar model for the description of experimental adsorption isotherms in geometrically heterogeneous pores. Let us now consider a quite narrow pore of the width H ∗ = 5 at T ∗ = 0.8. The density of the deposited matrix layer is given ∗ = 1 and ρ ∗ = 0.5. In Fig. 5a we comby Eq. (24) with zm m0 pare the plots of the average densities of adsorbed fluid vs the configurational chemical potential, μ in the pore with and without matrix. The fluid–pore wall interactions are characterized ∗ = 9 and z∗ = 0.5 and the fluid–matrix by the parameters εgs 0 potential is again of hard-sphere type. Fig. 5a displays the density profiles for dense and rarefied adsorbed phases. Note that the profiles for the system with the matrix have been multiplied by 20. The differences between both systems are striking. Whereas for the system without matrix we observe the capillary condensation, the presence of the matrix suppresses the condensation and the transition in the system is the capillary evaporation, i.e. the phase transition takes place at the chemical potential higher than at the bulk liquid–vapor transition. The location of the bulk transition is marked in Fig. 5a by vertical dotted line. Note that in addition to the equilibrium transition points for confined fluids, we have also included the hysteresis loops. The capillary evaporation occurs when the fluid–pore walls interaction is repulsive (or weakly attractive). The presence of matrix in the nearest vicinity of the pore walls causes that the adsorbate particles are moved towards the pore interior and thus their interaction with the walls becomes weaker. Obviously, the density of the matrix particles must be high enough in order to reduce the fluid–wall potential energy and to affect the fluid phase behavior. The occurrence of capillary evaporation at the ∗ = 0.5 means that a similar phenomenon will matrix density ρm0 be also observed at higher matrix densities. The density profiles in the system with the matrix exhibit only very small maxima at the pore walls for all bulk densities up to the bulk liquid–vapor coexistence. However, for the system without matrix the capillary condensation leads to a dense adsorbed phase that shows a pronounced layered structure. The layering transitions, if exist, appear at reduced temperatures lower than 0.8. Examples of the isotherms and profiles displayed in Figs. 4 and 5 indicate that the presence of matrix can substantially change the phase behavior of the confined fluid. In order to explore this problem we have evaluated phase diagrams for sev-

47

(a)

(b) Fig. 5. (a) Average densities inside the pore of the width H ∗ = 5 at T ∗ = 0.8 versus the configurational chemical potential. (—) Is for the system without ma∗ =1 trix, (- - -) for the pore with the matrix having the distribution (24) with zm ∗ = 0.5. The diameter of the matrix particles is the same as the diameter and ρm0 of adsorbate molecules and the interactions between matrix and fluid molecules are of hard-sphere type. The values of the parameters of the gas–solid potential ∗ = 9 and z∗ = 0.5. Vertical (-·-·) lines denote the equilibrium transitions. are εgs 0 (· · ·) Vertical line indicates the chemical potential at the bulk liquid–vapor transition. (b) Shows the density profiles for the systems from part (a). (—) Have been evaluated for the pore without matrix and (- - -) with the same matrix as in (a). There are two set of profiles for each system. Lower curves are at the ∗ = 0.0015, upper—for ρ ∗ = 0.0056. Note that the profiles bulk density ρb,s b,s for the system with the matrix have been multiplied by 20.

48

P. Bryk et al. / Journal of Colloid and Interface Science 313 (2007) 41–52

Fig. 6. Phase diagrams in the average density–temperature (left panel) and the chemical potential–temperature (right panel) planes for the fluid inside the pore of the width H ∗ = 5 with hard-sphere matrix deposited at each pore wall. The matrix particles have the same diameter as fluid molecules and their distribution ∗ = 0 (no matrix) (- - -), 0.1 (-·-·-), is described by Eq. (24) with zm = 1 and ρm0 0.19 (—), 0.5 (-··-··) line. (· · ·) In right panel shows the bulk phase diagram. The ∗ = 9 and z∗ = 0.5. Fluid–matrix parameters of the gas–solid potential are: εgs 0 interactions are of hard-sphere type.

eral model systems with hard-sphere fluid–matrix interactions, see Figs. 6 and 7. Left panels present the phase diagrams in the average density inside the pore-temperature plane, whereas right panels—the diagrams in the temperature—chemical potential plane. Fig. 6 shows the results for narrow pore, H ∗ = 5, and for ∗ εgs = 9 and z0∗ = 0.5. The matrix distribution is described by ∗ = 1, whereas the values of ρ ∗ are given in Eq. (24) with zm m0 the figure caption. In the right panel we have also included the bulk fluid phase diagram. For the system without the matrix ∗ = 0.1 we observe capillary condensation—the reland for ρm0 evant phase diagrams are shifted towards lower values of the chemical potential compared with those for bulk fluid. The crit∗ = 0.1 is only slightly lower than for ical temperature for ρm0 ∗ = 0.5 only the capillary the system without matrix. For ρm0 evaporation occurs within the investigated range of tempera∗ = 0.19 we observe a crossover between tures. However, for ρm0 capillary condensation and capillary evaporation regime. The entire μ/kT vs T ∗ curve lies very close to that for the bulk fluid (except for the fact that the critical temperature is much lower than for the bulk system), but at reduced temperatures lower than ≈0.88 the system undergoes capillary evaporation, whereas at higher temperatures we observe capillary condensation. The crossover is illustrated in the inset to the right panel and its occurrence is connected with the change of the ability of adsorbate to wet the pore wall. The presence of hard-sphere matrix inhibits wetting. Therefore, the investigation of wetting transitions at a single wall would be helpful in a precise determination of the crossover between capillary condensation and capillary evaporation regimes.

Fig. 7. Phase diagrams in the average density–temperature (left panel) and the chemical potential–temperature (right panel) planes for the fluid inside the pore of the width H ∗ = 8 with hard-sphere matrix deposited at each pore wall. The matrix particles distribution is described by Eq. (24) with zm = 1 and ∗ = 0 (no matrix) (- - -), 0.05 (-·-·) and 0.2 (-··-··) line. In all thees cases ρm0 the matrix particles have the same diameter as fluid molecules. Solid line is for ∗ = 0.003125, but for the matrix diameter four times larger than the diamρm0 eter of adsorbate molecules. Dotted line in right panel shows the bulk phase ∗ = 12 and z∗ = 0.5. diagram. The parameters of the gas–solid potential are: εgs 0 Fluid–matrix interactions are of hard-sphere type.

An increase of the fluid–wall energy increases ability of a liquid to wet the wall covered with matrix particles. All the capillary condensation diagrams presented in Fig. 7 are obtained ∗ = 12. Similarly as previously the matrix is built of hard for εgs spheres, distributed according to Eq. (24) and the fluid–matrix interaction is of hard-sphere type. Except for one case when ∗ = 2, the size of the matrix particles is the same as the size σmm of adsorbate molecules. The pore width is H ∗ = 8. The layering transition within the layer adjacent to the pore wall does not appear only for the system involving bigger matrix particles. For all remaining systems we observe the occurrence of the layering and its critical temperature and the triple point temperature between layering and capillary condensation shift down with an increase of the matrix density. The capillary condensation critical temperature depends on the size of the matrix particles. For example, the critical temperature is ∗ = 1 than for the systems with the matrix parlower when σmm ∗ = 2. This can be explained by the fact ticles of the size σmm that in the latter systems the capillary condensation is preceded by the layering transition and the capillary condensation takes place in effectively narrower “pore.” We now consider the case when the matrix–fluid interactions are given by the Lennard-Jones (12, 6) potential. For the sake of brevity only some selected results are presented below. Previous calculations (cf. Fig. 5a) have indicated that for hardsphere matrix–fluid interaction the capillary condensation may be superseded by the capillary evaporation. An attraction between matrix and fluid molecules diminishes the value of the

P. Bryk et al. / Journal of Colloid and Interface Science 313 (2007) 41–52

(a)

49

(b)

(c) Fig. 8. Adsorption isotherms (a) and the phase diagrams (b and c) inside the pore of H ∗ = 8 with walls covered by matrix interacting with fluid particles via the ∗ = 1. The isotherms in part (a) were calculated for ρ ∗ = 0.2, ε ∗ = 0.5 (-·-·) and Lennard-Jones potential. The matrix distribution is given by Eq. (24) with zm ms m0 ∗ = 0.5, ε ∗ = 1.5 (—). The loops of hysteresis are plotted here and the equilibrium transition is shown as dotted lines. Long dotted vertical line denotes the locus ρm0 ms ∗ = 0.2, ε ∗ = 0.5 (- - -), ρ ∗ = 0.5, ε ∗ = 1.5 (—), ρ ∗ = 0.5, ε ∗ = 0.5 of bulk transition point. The nomenclature of the lines in (b) and (c) is as follows: ρm0 ms ms ms m0 m0 ∗ = 0.2, ε ∗ = 1.5 (· · ·). The fluid–solid interactions are characterized by ε ∗ = 12 and z∗ = −0.5. (-··-··) In (c) denotes the bulk phase diagram. (-·-·) and ρm0 ms gs 0

chemical potential at the transition point and thus causes that capillary condensation phenomenon is privileged, cf. Fig. 8a. Moreover, distribution of the matrix particles causes “softness” of the pore walls and a “spatial distribution” of the attraction. Consequently, layering transition disappears (or is shifted towards much lower temperatures), see Figs. 8b and 8c.

The inclusion of attractive matrix–fluid forces is more efficient in deteriorating the layering transition than the effect of an increase of the matrix density in the case of solely repulsive matrix–fluid interactions. The same observation was the foundation of the approach of Ustinov et al. [19] that was used to describe experimental adsorption on geometrically heteroge-

50

P. Bryk et al. / Journal of Colloid and Interface Science 313 (2007) 41–52

(a)

(b)

(c) ∗ = 12 and z∗ = 0.5. Hard-sphere matrix is of the thickness z∗ = 2 and distributed according to Eq. (24). The Fig. 9. Adsorption of trimers in pore of H ∗ = 8, εgs m 0 size of matrix particles is the same as the size of timers’ segments. (a) Shows examples of adsorption isotherms at T ∗ = 0.8. (· · ·) Vertical line shows the chemical ∗ = 0.05, whereas (-··-··) potential for the bulk phase transition. (—) Is the isotherm for the pore without matrix, (- - -) is the isotherm with the matrix of density ρm0 ∗ 4 line abbreviates the matrix density ρm0 = 0.5. The latter isotherm has been multiplied by 10 for the values of the chemical potentials lower than the chemical potential at the bulk liquid–vapor transition (marked by (· · ·) vertical line). (b) Left panel shows the total segment profiles ρs∗ (z) at the bulk liquid–vapor coexistence ∗ = 0.5. Right panel (- - -) and for the rarefied phase at the evaporation transition, μ/kT = −9.95, calculated along the adsorption isotherm from (a), obtained for ρm0 ∗ (z), i = 1 and i = 2 for the dense phase at the capillary evaporation point, μ/kT = −9.8, calculated shows the end (—) and the middle-segment (- - -) profiles, ρsi ∗ = 0.5. (c) Displays phase diagrams in the average segment density–temperature (left panel) and in the chemical potential–temperature planes. (—) Is for for ρm0 ∗ = 0.1, (- - -) for ρ ∗ = 0.3, (-·-·) for ρ ∗ = 0.5, and finally, (-··-··) for ρ ∗ = 0.05, but nor z∗ = 2. (· · ·) In right panel is the bulk phase diagram. ρm0 m m0 m0 m0

neous surfaces. In order to derive a density functional approach, Ustinov et al. introduced the concept of “an effective density.” This density takes into account the presence of solid adsorbent particles. In the case of spherically symmetric molecules their approach is formally quite similar to the theory of Schmidt an co-workers. However, despite formal similarity the theory of

Ustinov et al. [19] leads to some inconsistencies, namely the excess free energy, as defined by Ustinov et al. is nonzero for zero adsorbate local density, because the second term in their counterpart of Eq. (10) is omitted. This leads to violation of the Gibbs adsorption equation. Obviously, functional derivative of the second term in Eq. (10) with respect to fluid local density is

P. Bryk et al. / Journal of Colloid and Interface Science 313 (2007) 41–52

zero and the occurrence of the term Φ˜ hs (r; {nα,s ≡ 0}, {nα,m }) in the free energy functional does not change the resulting density profile equation. Finally, we present examples of the application of the outlined theory to the case of adsorption of tangentially jointed trimers. The segment–matrix interaction is assumed to be of hard-sphere type and the matrix distribution describes Eq. (24) ∗ = 2. Fig. 9a shows examples of adsorption isotherms with zm ∗ = 12 and z∗ = 0.5 at T ∗ = 0.8. The evaluated for H ∗ = 8, εgs 0 plotted curves correspond to three characteristic cases. Solid line shows the results for the pore without any matrix. This curve exhibits two discontinuities, the first of which (at lower chemical potential) is associated with the layering transition, i.e. with the filling of the layers adjacent to the pore walls, the second—corresponds to the capillary condensation. Similarly as for adsorption of spherical molecules the presence of ∗ = 0.05, suppresses the layering and a low-density matrix, ρm0 ∗ = 0.5, causes that the presence of a high-density matrix, ρm0 the phenomenon of capillary condensation is replaced by the capillary evaporation. The data from Fig. 9a prove that quantitatively the behavior of short chains in pores with rough walls is similar as the behavior of spherical adsorbates (cf. Figs. 4a and 5a). Fig. 9a shows examples of the density profiles. Left panel ∗ = 0.5 displays total density profiles, ρs∗ (z), obtained for ρm0 at the chemical potential corresponding to bulk liquid–vapor transition point and for the rarefied adsorbed phase at the capillary evaporation point. Although the latter profile has been calculated at a higher value of μ/kT , the height of its maximum is much lower than the maximum of the former profile. It means that up to the capillary evaporation point the interior of the pore contacting with a bulk liquid phase is “more dry” than the pore being in contact with a gaseous phase at chemical potentials close to the bulk liquid–vapor transition. Right panel of Fig. 9c shows the profiles of end and middle segments for dense phase at the capillary condensation point. Here the ∗ = 0.05. Although the profiles ρ ∗ (z) and matrix density is ρm0 s1 ∗ ρs2 (z) differ slightly, especially close to the pore centre, there is no favorable orientation of the trimers inside the pore. Finally, Fig. 9c displays phase diagrams for several model systems depicted in the figure caption. The influence of the matrix on the capillary phase diagrams is quantitatively similar as in the case ∗ = 0.3 and 0.5 we of spherically symmetric adsorbates. For ρm0 observe capillary evaporation, in all remaining cases—capillary condensation. The theory presented in this work can be also applied to study adsorption of chain molecules in pores with rough walls. In particular, the theory can be generalized to the cases when consecutive segments of the chain molecule are of different size and differently interact with matrix particles. Also, one can imagine a model in which the matrix is built of quenched chain particles. Formally, the presented approach can be also extended to the case of a polydisperse matrix by employing the methods from Ref. [73]. The assumption of matrix polydispersity would probably lead to a more realistic model of geometrically heterogeneous adsorbents.

51

4. Conclusions Let us briefly summarize the results. The presence of a quenched component influences the phase behavior of the confined fluid. In the case of spherically symmetric adsorbates the main effect caused by the presence of a matrix relies on abolishing (or lowering the critical temperatures) the layering transitions. From the calculations performed it results that a version of density functional approach applied to the description of experimental data should be based on FMSA rather than on MFA. The former approximation significantly improves the accuracy of evaluated capillary phase diagrams. In general, answering the question how to model the surface roughens in order to observe desired effects is not simple. In general, the distribution of matrix particles does not need to be a stepwise function as in the present study. The size and the interaction of matrix particles with adsorbate molecules can also be different and their change may lead to different effects, e.g. one can observe a capillary evaporation instead of capillary condensation. The presence of matrix–fluid attractive forces more “effectively removes” the layering transitions and thus leads to “smoother” isotherms comparing with those for the system with harshly repulsive matrix–fluid interactions. Acknowledgments This work was supported by EU as TOK contract 509249. The research of Al.M. has been also partially supported by the Ministry of Education, Youth, and Sports of the Czech Republic under Project No. LC 512 (Center for Biomolecules and Complex Molecular Systems) and by the Grant Agency of the Czech Republic under Project No. 203/06/P432. References [1] W.A. Bakaev, W.A. Steele, Langmuir 8 (1992) 1372. [2] G. Chmiel, K. Karykowski, A. Patrykiejew, W. R˙zysko, S. Sokołowski, Mol. Phys. 81 (1994) 691. [3] J.L. Ricardo, W.A. Steele, A.J. Ramirez Cuesta, G. Zgrablich, Langmuir 13 (1997) 1064. [4] B. Kuchta, P. Llewellyn, R. Denoyel, L. Firlej, J. Low Temp. Phys. 29 (2003) 880. [5] E.A. Ustinov, D.D. Do, Langmuir 20 (2004) 3791. [6] P. Kowalczyk, K. Kaneko, A.P. Terzyk, H. Tanaka, H. Kanoh, P.A. Gauden, Carbon 42 (2004) 1913. [7] T. Miyata, A. Endo, T. Yamamoto, Mol. Simul. 30 (2004) 353. [8] D.D. Do, H.D. Do, Mol. Simul. 31 (2005) 651. [9] T.X. Nguyen, S.K. Bhatia, D. Nicholson, Langmuir 21 (2005) 3187. [10] D.D. Do, H.D. Do, J. Phys. Chem. B 110 (2006) 17531. [11] X.D. Liu, X.C. Lu, Q.F. Hou, K. Yang, Colloids Surf. A 281 (2006) 51. [12] G. Calleja, B. Coto, A.M. Morales-Cas, Appl. Surf. Sci. 252 (2006) 4345. [13] A. Wongkoblap, D.D. Do, J. Colloid Interface Sci. 297 (2006) 204. [14] C.G. Sonwane, C.W. Jones, P.J. Ludovice, J. Phys. Chem. B 109 (2005) 23395. [15] T.X. Nguyen, S.K. Bhatia, J. Phys. Chem. B 108 (2004) 14032. [16] L.D. Gelb, Mol. Phys. 100 (2002) 2049. [17] B. Coasne, R.J.-M. Pellenq, J. Chem. Phys. 120 (2004) 2913. [18] E.A. Ustinov, D.D. Do, M. Jaroniec, Appl. Surf. Sci. 252 (2005) 548. [19] E.A. Ustinov, D.D. Do, M. Jaroniec, Langmuir 22 (2006) 6238. [20] J.L. Harden, D. Andelman, Langmuir 8 (1992) 2547.

52

P. Bryk et al. / Journal of Colloid and Interface Science 313 (2007) 41–52

[21] R.R. Netz, D. Andelman, Phys. Rev. E 55 (1997) 687. [22] K. Grabowski, A. Patrykiejew, S. Sokołowski, E.V. Albano, A. de Virgiliis, Surf. Sci. 448 (2000) 11. [23] P. Bryk, D. Henderson, S. Sokołowski, J. Chem. Phys. 110 (1999) 15017. [24] P. Bryk, D. Henderson, S. Sokołowski, Langmuir 15 (1999) 6026. [25] J. Fischer, M. Methfessel, Phys. Rev. A 32 (1980) 2836. [26] R. Evans, in: D. Henderson (Ed.), Fundamentals of Nonuniform Fluids, Dekker, New York, 1992, chap. IV. [27] Yu. Duda, S. Sokołowski, P. Bryk, O. Pizio, J. Phys. Chem. B 102 (1998) 5490. [28] O. Pizio, S. Sokołowski, Phys. Rev. E 58 (1998) 2652. [29] P. Bryk, O. Pizio, S. Sokołowski, Mol. Phys. 95 (1998) 311. [30] M. Schmidt, Phys. Rev. E 66 (2002) 041108. [31] H. Reich, M. Schmidt, J. Stat. Phys. 116 (2004) 1683. [32] M. Schmidt, J. Phys. Condens. Matter 17 (2005) S3481. [33] A.J. Archer, M. Schmidt, R. Evans, Phys. Rev. E 73 (2006) 011506. [34] L. Lafuente, J.A. Cuesta, Phys. Rev. E 74 (2006) 041502. [35] D.D. Do, H.D. Do, Adsorpt. Sci. Technol. 23 (2005) 267. [36] O. Sinanoglu, K.S. Pitzer, J. Chem. Phys. 32 (1960) 1279. [37] E.A. Ustinov, D.D. Do, Part. Part. Syst. Char. J. 21 (2004) 112. [38] P.I. Ravikovitch, A.V. Neimark, Colloids Surf. A 187–188 (2001) 11. [39] P.I. Ravikovitch, A.V. Neimark, J. Phys. Chem. B 105 (2001) 6817. [40] P.I. Ravikovitch, A. Vishnyakov, A.V. Neimark, Phys. Rev. E 64 (2001) 011602. [41] W. R˙zysko, S. Sokołowski, O. Pizio, Z. Sokołowska, J. Colloid Interface Sci. 219 (1999) 184. [42] P.I. Ravikovitch, A.V. Neimark, Langmuir 22 (2006) 11171. [43] Y.-X. Yu, J.Z. Wu, J. Chem. Phys. 117 (2002) 2368. [44] D. Cao, J. Wu, J. Chem. Phys. 121 (2004) 4210. [45] M. Maeda, M.N. Kohonen, H.K. Christenson, J. Phys. Chem. B 105 (2001) 5906.

[46] D.D. Do, H.D. Do, Chem. Eng. Sci. 60 (2005) 1977. [47] J.W. Jiang, J. Phys. Chem. B 110 (2006) 8670. [48] M.B. Sweatman, N. Quirke, W. Zhu, R. Kapteijn, Mol. Simul. 32 (2006) 513. [49] T.J.H. Vlugt, R. Krishna, B. Smit, J. Chem. Phys. B 103 (1999) 1102. [50] D. Fu, J. Wu, Ind. Eng. Chem. Res. 44 (2005) 1120. [51] Y.P. Tang, J. Chem. Phys. 118 (2003) 4140. [52] Y. Rosenfeld, Phys. Rev. Lett. 63 (1989) 980. [53] R. Roth, R. Evans, A. Lang, G. Kahl, J. Phys. Condens. Matter 14 (2002) 12063. [54] Y.X. Yu, J.Z. Wu, J. Chem. Phys. 117 (2002) 10165. [55] Y.P. Tang, J.Z. Wu, Phys. Rev. E 70 (2004) 011201. [56] Y.P. Tang, J.Z. Wu, J. Chem. Phys. 119 (2003) 7388. [57] Z.D. Li, D.P. Cao, J.Z. Wu, J. Chem. Phys. 122 (2005) 224701. [58] D. Fu, J. Chem. Phys. 124 (2006) 164701. [59] Q. Yan, J.J. de Pablo, J. Chem. Phys. 111 (1999) 9509. [60] M.S. Wertheim, J. Stat. Phys. 42 (1986) 459. [61] M.S. Wertheim, J. Stat. Phys. 42 (1986) 477. [62] J.D. Weeks, D. Chandler, H.C. Andersen, J. Chem. Phys. 54 (1971) 5237. [63] J.A. Barker, D. Henderson, J. Chem. Phys. 47 (1967) 4714. [64] R.L. Cotterman, B.J. Schwarz, J.M. Prausnitz, AIChE J. 32 (1986) 1787. [65] C.M. Silva, H. Liu, E.A. Macedo, Ind. Eng. Chem. Res. 37 (1998) 221. [66] Y.P. Tang, J. Chem. Phys. 116 (2002) 6694. [67] Y.P. Tang, J. Chem. Phys. 119 (2005) 204704. ˙ [68] P. Bryk, K. Bucior, S. Sokołowski, G. Zukocinski, J. Phys. Chem. B 109 (2005) 2977. [69] P. Bryk, O. Pizio, S. Sokołowski, J. Chem. Phys. 122 (2005) 194904. [70] N.B. Wilding, Phys. Rev. E 52 (1995) 602. [71] Y. Tang, J. Chem. Phys. 121 (2004) 10605. [72] J. Mi, Y. Tang, C. Zhong, Y.-G. Li, J. Chem. Phys. 124 (2006) 144709. [73] O. Pizio, A. Patrykiejew, S. Sokołowski, Mol. Phys. 99 (2001) 57.