Models Based on Statistical Thermodynamics

Included in this group are methods that model the adsorption system in terms of forces acting between individual molecules.

Theoretical Background

Traditional adsorption theories attempt to describe experimental adsorption isotherms with an isotherm equation containing a small number of parameters. At a minimum, these parameters include the extent of the surface, such as the monolayer capacity (Qm), and the molar intensity of the gas-surface interaction, such as the Langmuir “K” constant or the BET “C” constant. In some equations, additional parameters take into account the lateral interaction of adsorbed molecules with each other. Other theories, such as the Dubinin-Astakhov approach, also include parameters for the effect of adsorbent porosity.

Instead of this classical kinetic or phenomenological approach, we can use a molecular-based statistical thermodynamic theory that allows us to relate the adsorption isotherm to the microscopic properties of the system: the fluid-fluid and fluid-solid interaction energy parameters, the pore size, the pore geometry, and the temperature.

The following example is given so that you may understand how such a theory is constructed:

A clean sample of a solid material containing slit-shaped pores of a single width is placed in an evacuated space. It is kept at a fixed temperature as a known quantity of pure argon gas is admitted into the space surrounding the sample. The pressure within the space is recorded over time. In this situation, the pressure falls rapidly from its initial value and gradually approaches a steady reading, called the equilibrium pressure. The amount adsorbed corresponds to the quantity of gas effectively removed from the gas phase by the solid surface. A graph that plots amount adsorbed versus equilibrium pressure is called an adsorption isotherm.

Under such conditions, the argon atoms that randomly enter the pore space feel the presence of the solid surface as the action of an external attractive force (the dispersion forces or Van der Waal’s forces) and spend more time near the surface. As a result, the space near the surface acquires a greater average density of argon atoms than regions farther removed.

If the equilibrium distribution of the gas atoms near the surface could be described as a function of pressure and the molecular properties of the components of the system, then a model could be constructed for the adsorption isotherm for the system. Modern physical chemistry provides several ways to calculate this distribution. All these methods are based on the fundamental thermodynamic law that such a system adopts a configuration of minimum free energy at equilibrium. Also needed is a description of the pairwise interaction energy between atoms, U(s), commonly given by a Lennard-Jones potential:

where

ε = a characteristic energy of the adsorptive,
σ = the diameter of the adsorptive molecule, and
s = the separation distance.

Molecular Simulation Methods

Two simulation techniques are commonly used to determine the distribution of gas molecules in a system in equilibrium: the molecular dynamics method and the Monte Carlo method. Both of these are used as reference methods because their results are considered exact.

Molecular Dynamics Method

In the molecular dynamics method, the position and velocity of individual gas particles are calculated over time at very short intervals. This method takes into account both the forces acting between the gas particles themselves and those acting between the gas particles and the atoms of the simulated surface. As the simulated particles collide with each other and with the surface, the average concentration of particles in the space near the surface is calculated; this calculation yields the amount of gas adsorbed.

This method can be thought of as a way to determine the chronological record of the movement of each particle in the system using time steps of 10-14 seconds. Although the mathematics are simple, the number of calculations required for a system of even a few hundred particles is astronomical and challenges even the fastest computers.

Monte Carlo Method

In the Monte Carlo method, determination of the system equilibrium distribution begins with an assumption (which may be only approximate) about the initial configuration of particles in the system. The system is “equilibrated” through a process of randomly selecting one particle and conditionally moving it a random distance in a random direction.

If the move results in a configuration of lower total energy, then the move is completed and another particle is randomly selected to be moved.

If the move results in a configuration of higher energy, a probability for that event is calculated, and a random number between zero and one is generated. If the generated number is smaller than the probability of the event, then the move is accepted; otherwise, another particle is selected and the process is repeated. This process continues until the average total energy of the system no longer decreases; at this point, average configuration data are accumulated to yield the mean density distribution of particles in the system.

Monte Carlo simulations require considerably less computation time than molecular dynamic simulations and can yield the same results; however, neither method provides a really practical way to calculate complete isotherms.

Density Functional Formulation

Density functional theory offers a practical alternative to both molecular dynamic and Monte Carlo simulations. When compared to reference methods based on molecular simulation, this theory provides an accurate method of describing inhomogeneous systems yet requires fewer calculations. Because the density functional theory provides accuracy and a reduced number of calculations, it is the basis embodied in the DFT models.

The system being modeled consists of a single pore represented by two parallel walls separated by a distance H. The pore is open and immersed in a single component fluid (adsorptive) at a fixed temperature and pressure. Under such conditions, the fluid responds to the walls and reaches an equilibrium distribution. In this condition (by the definition of equilibrium), the chemical potential at every point equals the chemical potential of the bulk fluid. The bulk fluid is a homogenous system of constant density; its chemical potential1 is determined by the pressure of the system using well-known equations. The fluid near the walls is not of constant density; its chemical potential is composed of several position-dependent contributions that must total at every point to the same value as the chemical potential of the bulk fluid.

As noted previously, at equilibrium, the whole system has a minimum (Helmholtz) free energy, known thermodynamically as the grand potential energy (GPE). Density functional theory describes the thermodynamic grand potential as a functional of the single-particle density distribution; therefore, calculating the density profile that minimizes the GPE yields the equilibrium density profile. The calculation method requires the solution of a system of complex integral equations that are implicit functions of the density vector. Since analytic solutions are not possible, the problem must be solved using iterative numerical methods. Although calculations using these methods still require supercomputing speed, the calculation of many isotherm pressure points for a wide range of pore sizes is a feasible task. The complete details of the theory and the mathematics can be found in the papers listed under DFT Model References.

The following graphs and accompanying text illustrate the results of using density functional theory to predict the behavior of a model system.

Figure 1 shows the density profile for argon at a carbon surface as calculated by density functional theory for a temperature of 87.3 K and a relative pressure of about 0.5.

Figure 1. Density Profile for Argon on Carbon at 87.3 K and a Relative Pressure of 0.5

This figure represents a cross-section of the region near the surface. Note the layerwise distribution of adsorbate; the first monolayer is sharply defined and a third layer can be distinguished. The area under the profile curve represents the amount adsorbed per unit area at this pressure. The positions of the maxima are separated by a distance determined by the size of the adsorptive atom.

Given the density profile, the amount adsorbed at the stated pressure can be easily calculated as the integral over the profile. Repeating this calculation over a range of pressures yields the adsorption isotherm for the model. If the value of H is very large, the isotherm obtained corresponds to that of an external, or free, surface. If H is smaller, a range of pressures is reached where two minima exist for the grand potential, showing the presence of two metastable phases having different density distributions but the same chemical potential. The phase with the lower GPE is the stable one. As the pressure is increased, a point is reached where the other phase becomes the stable one. This phase transition reflects condensation of adsorbate in the pore; the pressure at which it occurs is called the critical pore-filling pressure. This pressure is analogous to the condensation pressure predicted by the Kelvin equation in the classical model of pore filling.

Figure 2 shows how the profiles change with pressure for a model pore with H = 40 angstroms. The inset shows the density profiles for the corresponding points of the isotherm.

Figure 2. Model Isotherm for Argon at 87.3 K in a 40 Å Slit in a Carbon Substrate

The profiles show the density distribution from one wall to the center of the slit; the other half of the distribution is a mirror image of the profile shown.

As the pressure is first increased from zero, almost all the adsorbed atoms occupy a position close to the surface.

  • Inset a shows the profile corresponding to point a on the isotherm where the surface is about half covered.
  • At point b, the first layer is so full that it is more favorable for atoms to start a new layer.
  • At point c, a third layer is forming. Point c, for this size slit, is the critical pore-filling pressure. In inset c, the profile shows the density decreasing to near zero (actually the bulk gas density) at 4 or 5 molecular diameters from the surface.
  • Inset d shows the profile converging on a density similar to that of bulk liquid argon in the center of the pore, indicating a phase transition.

Note that the adsorption isotherms for pores larger than the one shown in the previous graph is identical up to point c. The lower branch of the isotherm simply continues to a higher pressure for larger pores. This trend is illustrated in the Figure 3, where isotherms for some larger size pores are shown. It is clear that pore size is uniquely characterized by a corresponding critical pore-filling pressure. At large pore sizes, density functional theory produces results for the critical filling pressures that are in good agreement with those produced by the Kelvin equation.

Figure 3. Model Isotherms for Some Larger Pore Widths Argon on Carbon at 87.3 K

Figure 4 shows model isotherms for pores in the micropore size range. Note the logarithmic scale for pressure.

Figure 4. Model Isotherms in the Micropore Size Range of Pore Width Argon on Carbon at 87.3 K

Pores of 4 Å width, barely larger than the argon atom (3.38 Å), fill at pressures below 1 millitorr. Pores below 15 Å fill before a monolayer is completed on the surface of the larger pores. In the micropore size range, the pore volume fills more gradually with pressure and the total shape of the isotherm is important in characterizing the pore size.

Models Included

Non-Local Density Functional Theory with Density-Independent Weights

N2 - DFT Model
AR - DFT Model

Geometry:

Slit

Substrate:

Carbon (graphite)

Category:

Porosity

Method:

Nitrogen at 77 K; Argon at 87 K

Using the methods of non-local density functional theory, two sets of isotherms have been calculated to serve as kernel functions for the characterization of porous solids from adsorption data. The model isotherms are stored in binary format files. These models assume a slit-like pore geometry. The pore size range from 4.0 to 4000 Å is covered in 91 classes in a geometric progression. The class intervals are rounded to the nearest 0.02 molecular diameters. A model for the free or external surface is included to account for unfilled pores. Each of the 92 model isotherms has been calculated at 181 pressure points from near 1×10-6 to near 1.00 relative pressure.

These models are identical to those supplied with the original DOS version of DFT software. Some slight difference from the DOS results may be noted when they are applied to the same data due to improvements in the deconvolution algorithm and better regularization of the current software.

Non-Local Density Functional Theory with Density-Dependent Weights

N2 - Modified Density Functional

Geometry:

Free surface

Substrate:

Surface energy

Method:

Nitrogen at 77K

Using the modified Tarazona prescription described by Olivier (see DFT Model References [items 3 and 4]), model isotherms were calculated for a wide range of adsorptive energies to a relative pressure of 0.6. The model makes no provision for pore filling in the micropore region. If the sample solid contains small mesopores, the isotherm data should be truncated (using the Select Data Points window) to a suitably low relative pressure to avoid trying to fit this region; mesopore filling reports as a large area of low energy in the calculated distribution of adsorptive potential.

The surface energy is reported in terms of the effective Lennard-Jones interaction parameter, i.e., for the adsorptive / adsorbent pair divided by Boltzmann constant. The units are therefore Kelvin.

N2 - Cylindrical Pores - Oxide Surface
AR - Cylindrical Pores - Oxide Surface

Geometry:

Cylinder

Substrate:

Oxide

Category:

Porosity

Method:

Nitrogen at 77K; Argon at 87K

Model isotherms were calculated using a combination of statistical mechanical calculations and experimental observations for macroporous silicas and MCM-41 mesoporous silicas as well as zeolites. The pore-filling pressures were determined as a function of the pore size from adsorption isotherms on MCM-41 materials characterized by X-ray and other techniques. The variation of the pore fluid density with pressure and pore size has been accounted for by density functional theory calculations. The N2 model reports pore sizes ranging from 3.8 to 387 Å and the AR model from 3.8 to over 500 angstroms.

References:

M. Jaroniec, M. Kruk, J.P. Olivier, and S. Koch, “A New Method for the Accurate Pore Size Analysis of MCM-41 and Other Silica-Based Mesoporous Materials,” Proceedings of COPS-V, Heidelberg, Germany (1999).

N2 – Cylindrical Pores – Pillared Clay Surface (Montmorillionite)

Geometry:

Cylinder

Substrate:

Crystalline Silicate

Category:

Porosity

Method:

Nitrogen at 77K

Model isotherms were calculated using a combination of statistical thermodynamic Non-Local Density Functional Theory (NLDFT) calculations and experimental isotherms for reference samples of montmorillionite. The construction method for the hybrid models was analogous to that described in the first reference below (Jaroniec et al,1999). The additional references add additional theoretical details as well as examples of the application of the model to pillared clay catalysts. This model reports pore widths from 3.8 to 387 angstroms.

References:

Mietec Jaroniec, Michal Kruk, James P. Olivier and Stefan Koch, “A New Method for the Characterization of Mesoporous Silicas,” Proceedings of COPS-V, 1999, Studies in Surface Science, Vol 128, Characterization of porous Solids V , Unger, et al, Eds, Elsevier, Amsterdam, 2000.

James P. Olivier and Mario L. Occelli, “Surface Area and Microporosity of a Pillared Interlayered Clay (PILC) from a Hybrid Density Functional Theory (DFT) Method,” The Journal of Physical Chemistry B; 2001, 105(3), 623-629.

M. L. Occelli, J. P. Olivier, J. A. Perdigon-Melon, and A. Auroux, “Surface Area, Pore Volume Distribution, and Acidity in Mesoporous Expanded Clay Catalysts from Hybrid Density Functional Theory (DFT) and Adsorption Microcalorimetry Methods,” Langmuir 2002, 18, 9816-9823.9b.

James P. Olivier, “The Importance of Surface Heterogeneity in Developing Characterization Methods.” 6th International Symposium on the Characterization of Porous Solids, Studies in Surface Science and Catalysis 144, Elsevier, 2002.

James P. Olivier and Mario L. Occelli, “Surface Area and Microporosity of Pillared Rectorite Catalysts from a Hybrid Density Functional Theory Method,” Microporous and Mesoporous Materials 2003, 57, 291-296.

C02 - DFT Model

Geometry:

Slit

Substrate:

Carbon

Category:

Porosity

Method:

Carbon dioxide at 273 K

Model isotherms were calculated using the non-local prescription of Tarazona, employing molecular parameters derived from the known bulk properties of carbon dioxide.

AR - Modified Density Functional Model

Geometry:

Free surface

Substrate:

Any

Category:

Surface energy

Method:

Argon at 87K

This model was produced in the same manner as the N2 Modified Density Functional model listed earlier, except applicable to argon adsorbed at 87.3 K.

N2 - Tarazona NLDFT, Esf = 30.0K

Geometry:

Cylinder

Substrate:

Oxide

Category:

Porosity

Method:

Nitrogen at 77K

Model isotherms were calculated using the prescriptions of Tarazona for density dependent weighting functions and a cylindrical pore geometry. The wall potential used is k = 30 K, typical for a silica or alumina surface.

This model file is particularly useful for sizing zeolites or zeolite containing materials that have substantial micropore volume. The reported pore size range is 3.8 to 387 angstroms.

References:

P. Tarazona, Phys. Rev. A 31: 2672 (1985).
Idem, Phys. Rev. A 32: 3148 (1985).
P. Tarazona, U. M. B. Marconi, and R. Evans, Mol. Phys. 60: 573 (1987).

N2 - Carbon Slit Pores by NLDFT
Ar - Carbon Slit Pores by NLDFT

Geometry:

Slit

Substrate:

Carbon

Category:

Porosity

Method:

Nitrogen at 77K; Argon at 87K

Model isotherms were calculated using the prescriptions of Tarazona for density dependent weighting functions and a slit-like pore geometry. These models are slightly different from N2-DFT and Ar-DFT models that were calculated using NLDFT with density independent weighting functions.

The reported pore size range is from 3.5 to 1000 angstroms.

References:

P. Tarazona, Phys. Rev. A 31: 2672 (1985).
Idem, Phys. Rev. A 32: 3148 (1985).
P. Tarazona, U. M. B. Marconi, and R. Evans, Mol. Phys. 60: 573 (1987).

N2 - Carbon Finite Pores, As=6, 2D-NLDFT
Ar - Carbon Finite Pores, As=6, 2D-NLDFT

Geometry:

Finite Slit

Substrate:

Carbon

Category:

Porosity

Method:

Nitrogen at 77K; Argon at 87K

Model isotherms were calculated using the prescriptions of Tarazona for density dependent weighting functions assuming 2D model of finite slit pores having a diameter-to-width aspect ratio of 6.

This model is particularly useful for microporous carbon materials. The reported pore size range is from 3.5 to 250 angstroms.

References:

Jacek Jagiello and James P. Olivier. “A simple two-dimensional NLDFT model of gas adsorption in finite carbon pores. Application to pore structure analysis.,” The Journal of Physical Chemistry C, 113(45):19382-19385, 2009.

N2 - Carbon Finite Pores, As=12, 2D-NLDFT
Ar - Carbon Finite Pores, As=12, 2D-NLDFT

Geometry:

Finite Slit

Substrate:

Carbon

Category:

Porosity

Method:

Nitrogen at 77K; Argon at 87K

Model isotherms were calculated using the same methods and assumptions that were used in the model above except in this model, the aspect ratio is equal to 12.

These two finite pore models may be used as a research tool in conjunction with independent analytical techniques such as high-resolution transmission electron microscopy (HRTEM) and/or X-ray diffraction (XRD) to obtain comprehensive information about the structure of studied carbon material.

References:

Jacek Jagiello and James P. Olivier. “A simple two-dimensional NLDFT model of gas adsorption in finite carbon pores. Application to pore structure analysis.,” The Journal of Physical Chemistry C, 113(45):19382-19385, 2009.

N2 - Carbon Cylinder, single-wall nanotube by NLDFT
Ar - Argon Cylinder, single-wall nanotube by NLDFT

Geometry:

Cylinder

Substrate:

Carbon

Category:

Porosity

Method:

Nitrogen at 77 K; Argon at 87 K

Model isotherms were calculated using the prescriptions of Tarazona for density dependent weighting functions and cylindrical pore geometry. The pore wall potential is described by the Lennard-Jones potential of interaction between a gas molecule and the graphitic surface of an infinitely long cylinder.

This model is particularly useful for characterizing carbon single-wall nanotubes. The reported pore size range is from 3.5 to 1000 angstroms.

References:

P. Tarazona, Phys. Rev. A 31: 2672 (1985).
Idem, Phys. Rev. A 32: 3148 (1985).
P. Tarazona, U. M. B. Marconi, and R. Evans, Mol. Phys. 60: 573 (1987).

N2 - Carbon Cylinder, multi-wall nanotube by NLDFT
Ar - Argon Cylinder, multi-wall nanotube by NLDFT

Geometry:

Cylinder

Substrate:

Carbon

Category:

Porosity

Method:

Nitrogen at 77 K; Argon at 87 K

Model isotherms were calculated using the prescriptions of Tarazona for density dependent weighting functions and cylindrical pore geometry. The pore wall potential is described by the Lennard-Jones potential of interaction between a gas molecule and multiple concentric graphitic surfaces of infinitely long cylinders.

This model is particularly useful for characterizing carbon multi-wall nanotubes. The reported pore size range is from 3.5 to 1000 angstroms.

References:

P. Tarazona, Phys. Rev. A 31: 2672 (1985).
Idem, Phys. Rev. A 32: 3148 (1985).
P. Tarazona, U. M. B. Marconi, and R. Evans, Mol. Phys. 60: 573 (1987)

Ar - Zeolites H-Form by NLDFT

Geometry:

Cylinder

Substrate:

Zeolite

Category:

Porosity

Method:

Argon at 77 K

Model isotherms were calculated using the prescriptions of Tarazona for density dependent weighting functions and cylindrical pore geometry. The pore wall potential is described by the Lennard-Jones potential of interaction between a gas molecule and the oxide surface of an infinitely long cylinder.

This model is particularly useful for characterizing oxides and H+ and (NH4)+ exchanged zeolites. The reported pore size range is from 3.5 to 300 angstroms.

Ar - Zeolites Me-Form by NLDFT

Geometry:

Cylinder

Substrate:

Zeolite

Category:

Porosity

Method:

Argon at 77 K

Model isotherms were calculated using the prescriptions of Tarazona for density dependent weighting functions and cylindrical pore geometry. The pore wall potential is described by the Lennard-Jones potential of interaction between a gas molecule and the oxide surface of an infinitely long cylinder.

This model is similar to the model above, but it more appropriate is for characterizing alkali metal exchanged zeolites. The reported pore size range is from 3.5 to 300 angstroms.