Adsorption of hard sphere fluid in porous material : A Monte Carlo simulation approach for pressure calculation ∗

We have performed Monte Carlo simulations in the canonical ensemble of a hard-sphere fluid adsorbed in microporous media. The pressure of the adsorbed fluid is calculated by using an original procedure that includes the calculations of the pressure tensor components during the simulation. In order to confirm the equivalence of bulk and adsorbed fluid pressures, we have exploited the mechanical condition of equilibrium and performed additional canonical Monte Carlo simulations in a super system “bulk fluid + adsorbed fluid”. When the configuration of a model porous media permits each of its particles to be in contact with adsorbed fluid particles, we found that these pressures are equal. Unlike the grand canonical Monte Carlo method, the proposed calculation approach can be used efficiently to obtain adsorption isotherms over a wide range of fluid densities and porosities of adsorbent.


Introduction
Investigation of adsorption, diffusion and reaction of fluids in porous materials is of much interest to the basic science and to many applications [1].It is of relevance for applications in membrane separations and heterogeneous catalysis by zeolites.From a simplified view, zeolite structure is characterized by the presence of a regular network of pores and cavities, whose diameter is commensurate with the molecular dimensions [2,3].Continuum systems with quenched disorder include, for example, porous and sintered porous materials, gels, amorphous substances, glaces, engineering composites, etc.A more detailed description of porous systems can be found in [4].
Computer simulation continues to play an important role in developing a full understanding of fluid adsorption in porous materials [5,6].Monte Carlo (MC) and molecular dynamics studies of the fluid-solid interface probe behaviour are inaccessible to experiments and provide essentially exact data for well-defined model potentials.In the absence of accurate information concerning the adsorbate-adsorbate and adsorbate-solid interactions, such information is of crucial importance to the theoretical predictions [6][7][8][9][10][11][12].
The most commonly used technique for the adsorption study has been the grand canonical Monte Carlo (GCMC) method [5][6][7][8][9].In grand canonical ensemble, the temperature and chemical potential are imposed and the number of particles is allowed to fluctuate during the simulation.However, the GCMC method is limited in its range of applicability to the systems where fluid and solid packing fractions are not too large.In a high packing fraction regime (or/and if the fluid system of interest contains relatively large molecules) the probability of particle creations and destructions becomes very low, which leads to poor convergence of the MC averages.For chain molecules, for example, almost all attempts of insertion result in an overlap with one of the atoms of the molecules in the system and, therefore, such a move has an extremely low probability of acceptance.Configurational-bias GCMC technique proposed recently [5,9], does not sufficiently improve the situation at high packing fractions.
On the other hand, to study the fluid adsorption in microporous solid (matrix) the canonical ensemble MC (NVT MC) simulation has been used on a supersystem which contains both the fluid confined in the porous solid and the bulk fluid phase outside it [11,12].However, in such a simulation there is an interface between the bulk phase and the matrix.In this region the fluid properties differ from both the bulk and adsorbed fluid properties of interest.Since in the simulation the systems are small compared with the experimental (hence, the interfacial region is a relatively large part of such a system), we have to simulate an unnecessarily large system to minimize the effect of this interfacial region.Such a simulation would, of course, be appropriate if the interest lay in this very region.
Therefore, it is clearly of interest to develop a simulation approach suitable for studying the adsorption under a wide variety of states of the bulk fluid and matrix.This is the primary aim of the present work.
Theoretical tools for the investigation of the properties of fluids and fluid mixtures adsorbed in disordered porous materials have been mainly adapted from the liquid-state statistical mechanics and include integral equation approach, pioneered by Madden and Glandt [13,14].These authors have presented exact Mayer cluster expansions for the correlation functions for the case when the matrix (quenched species subsystem) is generated by a quench from an equilibrium distribution, as well as for the case of arbitrary distribution of obstacles.However, their integral equations for the correlation functions appeared to be only approximate.This has been explicitly demonstrated by Given and Stell [15,16].By performing a detailed analysis of the cluster series, they have derived a correct set of equations and called them the replica Ornstein-Zernike equations.These approaches have been applied for calculation of structural and thermodynamic properties of a few molecular models of fluids confined in disordered porous materials [7,8,[17][18][19][20].Of particular importance is the adsorption isotherm which describes how the average density of the adsorbed fluid changes with the changes in pressure or chemical potential of the bulk fluid at constant temperature.Since the mentioned integral equation theories deal with fixed fluid density in the porous material, a relationship between the chemical potential and this averaged density is needed to calculate the adsorption isotherm.Madden [14] and Fanti et al. [21] have suggested that an appropriate route to this is to calculate the pressure of the fluid in the porous material via a virial theorem and then to obtain the chemical potential via integration of the isothermal Gibbs-Duhem equation.Their expression for the virial pressure consists of the sum of kinetic term, fluid-fluid and fluid-matrix interaction terms; matrix-matrix interaction term is absent since matrix particle is fixed.
However, it was shown later on that the corresponding adsorption isotherm calculated via integration of the isothermal Gibbs-Duhem equation was inconsistent with MC simulation results for a specific example [17].The origin of this discrepancy has been revealed by Rosinberg et al. [22] who derived the correct expression for the virial pressure which consists of the additional terms which are not presented in the equation of Madden [14].In a parallel study, Dong [23] has suggested that the expression of the mechanical pressure could be obtained by considering the average of the microscopic force balance between the fluid inside and outside the porous matrix.In contrast to previous results, Dong's expression is simpler and consists of the kinetic term and only fluid-fluid interaction term.Moreover, it was argued in [23] that this pressure should coincide with the thermodynamic pressure, although no formal proof was given.Finally, Kierlik et al. [24] by using the Green-Bogolubov method for the thermodynamic pressure have shown that applying the condition of mechanical equilibrium between two systems does not reveal information about the thermodynamics of adsorbed fluid systems.They concluded, that the application of the condition of mechanical equilibrium does not provide a route to the calculation of the thermodynamics of these systems.
As can be seen, there is a controversy in the pressure calculations for an adsorbed fluid.The simulation approach we propose in this paper, is based on the calculation of the pressure of the fluid adsorbed inside the porous material using the virial-like expression, similar to that suggested by Fanti et al. [21] and Madden [14].Since this pressure is equal to the bulk fluid pressure (we confirmed this by additional NVT MC simulation of supersystem "bulk fluid + fluid in matrix") our method directly yields the bulk pressure versus adsorbed fluid density.The simple model and the simulation method are described in the next section.In section 3, as an illustration, simulation results are presented.Concluding remarks are given in the last section.

Model and simulation approach
In order to present a new technique of pressure tensor calculation during NVT MC simulation we consider a simple case of hard sphere (HS) fluid particles (f ), interacting through the potential U f f (r), that is, The fluid particles are confined to a matrix of HS particles (m) interacting through the potential where σ and σ m are the diameters of fluid and matrix species, respectively.The cross-interaction, U f m , is given in the same form as in equations ( 1) and ( 2) but with characteristic length σ f m , σ f m = 0.5(σ + σ m ).We assume that σ = 1, and σ is taken as the length unit.
Estimation of the bulk fluid pressure in molecular simulation is usually accomplished via calculation of the virial [5].The bulk pressure of hard sphere fluid in MC simulations can be obtained through the radial distribution function, g(r), evaluated at contact distance, r = σ + , where ρ = N/V is the fluid number density.The bulk fluid pressure is obtained at the end of MC simulations by extrapolation of g(r) to the contact.As has been shown before [25][26][27], the bulk pressure of the fluid with discontinues (hard sphere-like) potential can be calculated during the NVT MC simulation run.Besides, the procedure of pressure tensor calculation of nonadditive hard sphere mixture adsorbed in a slit-like pore [25] and of square well fluid with liquid-vapor interface [27] has been recently presented.
The component P bulk xx of the bulk fluid pressure, P bulk , can be calculated using the mechanical definition [28] where x ij = x i − x j , x i is the coordinate of atom i in X direction.The force component between atoms i and j in the X direction is For the hard sphere potential [26], As has been shown, the derivative of the potential can be obtained during MC simulations [26] if the Dirak δ function is calculated numerically, where Θ(x) is the step function.The approximation in equation ( 7) is valid when ∆ σ → 0. Substitution of equations (5-7) into equation ( 4) gives the P bulk xx component of pressure tensor for the hard sphere fluid, The double sum in equation ( 8) is calculated for the distances between particles i and j that lie within the interval σ r < σ + ∆ σ .In a simulation P bulk xx is obtained at different values of ∆ σ to estimate its value when ∆ σ → 0. In this work we have taken ∆ σ equal to 0.005, 0.01 and 0.015.The procedure of pressure value estimation when ∆ σ → 0, as well as its dependence on the choice of ∆ σ have been discussed in [27].
When HS fluid is adsorbed into the HS matrix, in order to calculate total fluid pressure component, P T xx , the additional term should be added to equation (8), which represents the sum of x-components of the matrix-fluid forces.Such a definition of the pressure treats the matrix-fluid interaction as a contribution to the intermolecular forces in a system consisting of the fluid and solid, rather than as an external field acting on the fluid.Therefore, the total fluid pressure component is The components P T yy and P T zz are obtained by applying the same procedure as that of P T xx .In homogeneous fluids the three components of the pressure are equivalent and the total pressure is

Applications
To illustrate the use of the simulation approach described in this work, the HS fluid adsorption isotherms in different matrices have been calculated.
The model adsorbent (matrix) has been prepared by a random distribution of nonoverlapping hard-sphere particles of diameter σ * m = σ m + L, inside the cubic simulation cell with periodic boundary conditions.From 30 to 160 matrix particles have been used, depending on the values of particle diameters and number density, ρ m .Such number of matrix particles is believed to be large enough to form a statistically representative sample of a macroscopic porous medium while it is small enough to accommodate a manageable number of adsorbed fluid molecules.The matrix particles are rigidly fixed in their positions.In the case of disordered matrix, for some cases we have repeated our calculations for three different matrix particle configurations but this does not change the results to a marked degree (at least for the presented parameters of the model).For high matrix packing fraction, η m = π/6ρ m σ 3 m , the particles were allocated using the fcc crystalline array.The extended matrix particle diameter σ * m , (when L 0), allows to form the channels between the matrix spheres such that the fluid particles can pass through and touch each matrix particle.If L = 0, some of the matrix particles may be "hidden" from the contact with the fluid (especially at high matrix packing fraction).In this case, we have found that P T xx = P T yy = P T zz .Although L in the interval σ > L > 0 can also guarantee the nonzero probability of fluid contact with each solid particle (depending on σ m ), in our calculations, presented in the next section, the parameter L has been chosen to be equal σ.
An initial configuration for the fluid is prepared by inserting the f -particles into the available space.For the high fluid density we apply the algorithm of particle insertion and growing [5].
In order to verify that the fluid pressure in bulk phase is the same as the pressure of an adsorbed fluid we have performed NVT MC simulation on a supersystem ("bulk fluid + matrix fluid").In this case, the simulation cell was a parallelepiped in shape with the central part occupied by the matrix particles [12].Periodic boundary conditions were imposed in all three directions and neighbor listing was used.Each f -particle was chosen in turn and randomly displaced for N eq configurations to equilibrate the system.The displacement parameter was chosen to provide an acceptance ratio around 40%.Equilibration of the supersystem is more time consuming than the adsorbed fluid system in cubic cell and was performed from over N eq = 10 5 to over 4 × 10 5 trials to displace each fluid particle and an equal number for obtaining the ensemble averages.Care was taken so that equilibration was reached.We have evaluated the symmetry of the density profiles with respect to the matrix center and the stability of the averaged bulk and adsorbed fluid density from run to run.For the fluid densities that were considered, we have taken care to obtain a sufficiently wide "bulk" region for the bulk fluid and fluid inside the matrix.We have defined the "bulk" regions of the cell as that in which the density profile is constant up to fluctuations.In that regions we have calculated the pressure by using equations 8 and 10.
For each of the isotherms presented below, we have performed two NVT MC simulations of the supersystem "bulk fluid + fluid in matrix" at different bulk fluid densities.In these simulations we have confirmed that the pressure of bulk fluid (calculated from equation 8) and the pressure of fluid adsorbed in the matrix (calculated from equation 12) are the same.Note that the simulation of the supersystem "bulk fluid + matrix fluid" is 6-9 times more timeconsuming than our new simulation approach at the same conditions.
Initially, in figure 1 we present the adsorption isotherms of HS fluid into the HS disordered medium with packing fraction η m = 0.3077 and different sizes of matrix particles, σ m .In both cases considered here we use σ * m = σ m +L, (L = σ), in order to form the channels inside the matrix.In order to obtain the porous medium of high packing fraction, the matrix of σ m = 3 has been formed by allocation of the particles with σ * m = 4 using the fcc crystalline array.Such allocation allows us to form the matrix with initial packing fraction η * m = πρ m σ * m /6 ≈ 0.74.Thus, the matrix with σ m = 3 has the fcc structure and η m = 0.3077.Meanwhile, the second matrix with σ m = 4 has been formed by random insertion of HS particles with σ * m = 5 into a cubic simulation cell.
As expected, the adsorption of a hard sphere fluid is slightly higher into a matrix made of the large particles.The difference in adsorption increases for higher bulk pressure.In figure 2 we study the adsorption of HS fluid as a function of matrix packing fraction, η m .The fluid density is fixed, ρ = 0.42035.As in the previous case one can see that adsorption is higher in a matrix made of large particles.When the packing fraction decreases this difference disappears.The matrices with the highest η m has the fcc structure (due to the reasons described above) while the other matrices are disordered.
In order to investigate the effect of matrix configuration on fluid adsorption we performed the following computer experiment.We prepared two matrices with the same packing fraction η m = 0.3745 and particle size σ m = 4.The first one has the fcc structure (σ * = 5).While the second matrix contains sheetlike pores (as in carbonate rocks or clays [1]), i.e. the matrix particles are randomly allocated on the n squares in Y Z plane.The distance between neighbouring square planes in X direction is L x = 4.7 which allows the entrance of fluid particles inside the pores and their possible contact with each matrix particle.In other words, the x i coordinate of i matrix particle must be equal to n × L x , (n = 0, 1, 2, ...).
As can be seen in figure 3, the configuration of matrix HS particles does affect the adsorption isotherm of HS fluid.Namely, at low pressure the fluid is better adsorbed into the sheet-like structure porous material than in the fcc structure porous material.This is just opposite at high pressure.
It is worth noting, that at some cases of the highest adsorbed fluid density presented above the value of the total packing fraction in the simulation box, η T = η m + η, is around 0.55-0.6.As well known, it is very problematic or even impossible to use GC MC simulation at these values of packing fractions [5].

Conclusions
In this paper the adsorption of a hard sphere fluid in a hard sphere porous material has been studied applying a new approach to the pressure calculation during NVT MC simulation.This approach has two features which make it a useful tool for studying the adsorption in porous material.Firstly, one omits the calculation of bulk fluid pressure.Secondly, this method may be efficiently used over a wide range of conditions.The first advantage makes this method desirable for the investigation of adsorption from mixture and adsorption of chain molecules.Besides, the proposed method is considerably less timeconsuming than other MC methods that used to be applied for adsorption investigations.
We should stress that, the present NVT MC approach, like all simulation methods, has also got its own disadvantages.The method imposes certain restrictions on the model matrix particle configurations: the distance between matrix particles (size of channels) must permit the entrance of fluid particles.In other words, none of matrix particles can be "hidden" from the contact with fluid particle.Only then, there is an equivalence of bulk and adsorbed fluid pressures within our calculation approach.This might be a reason why the simulation investigation made by Vega, Kaminsky, and Monson [17] has shown that the virial expression of pressure given by Madden [14] fails to produce thermodynamically consistent results for the adsorption isotherm calculated from Gibbs-Duhem equation.
We are in the progress of extending the calculation reported here to a variety of other systems including fluid mixture, macromolecules, and more realistic model of the fluid-solid interactions.Our preliminary calculations for the dimer fluid adsorption in porous media give promising results [29] in a wide range of fluid densities and structure of matrix.

Figure 1 .
Figure 1.Adsorption isotherms of a hard sphere fluid in hard sphere matrices with the packing fraction η m = 0.3077.The matrix is made of the particles with the diameter σ m = 4 (stars) and σ m = 3 (circles).The accuracy of the calculations (error bar) does not extend the size of symbols.

Figure 2 .
Figure 2. Adsorption isotherms of a hard sphere fluid with density ρ = 0.42035 in a hard sphere matrices.The matrix is made of the particles with the diameter σ m = 4 (stars) and σ m = 3 (circles).The accuracy of the calculations (error bar) does not extend the size of symbols.

Figure 3 .
Figure 3. Adsorption isotherms of a hard sphere fluid in a matrix of hard spheres with diameter σ m = 4 and η m = 0.3745.The accuracy of the calculations (error bar) does not extend the size of symbols.