Phase behavior of patchy colloids confined in patchy porous media

Abstract

A simple model for functionalized disordered porous media is proposed and the effects of confinement on self-association, percolation and phase behavior of a fluid of patchy particles are studied. The media are formed by randomly distributed hard-sphere obstacles fixed in space and decorated by a certain number of off-center square-well sites. The properties of the fluid of patchy particles, represented by the fluid of hard spheres each bearing a set of the off-center square-well sites, are studied using an appropriate combination of the scaled particle theory for the porous media, Wertheim's thermodynamic perturbation theory, and Flory–Stockmayer theory. To assess the accuracy of the theory a set of computer simulations have been performed. In general, predictions of the theory appeared to be in good agreement with the computer simulation results. Confinement and competition between the formation of bonds connecting the fluid particles, and connecting fluid particles and obstacles of the matrix, gave rise to a re-entrant phase behavior with three critical points and two separate regions of the liquid–gas phase coexistence.

---

1. Introduction

Porous materials, both natural (e.g., zeolites, biological tissues, rocks and soil) and artificial (e.g., ceramics, synthetic zeolites, cement, and silica aerogels), are ubiquitous. They are relevant to a broad range of applications, e.g., catalysis and sensing, gas storage and separation, adsorption, drug delivery, biomaterial immobilization, environmental remediation, molecular recognition, energy storage, etc.^1–28 In most of the earlier applications of porous materials, conventional porous materials, such as zeolites, silica aerogels or activated carbon, were used.^29–35 Due to the emergence of new classes of porous materials, which include metal–organic frameworks,^36–40 porous organosilicas,^41,42 porous polymers^43–48 and covalent organic frameworks,^48–50 substantial advances in the field have been achieved over the last few decades. In particular, due to their well-developed porosity and tunable surface chemistry these composite organic–inorganic porous materials can be relatively easily functionalized, i.e. the properties of the inner surface of the pores can be appropriately modified by decorating it with different functional groups, organic and inorganic moieties, etc. Functionalized porous materials offer a number of new opportunities to achieve the desirable properties for a diverse range of industrial applications⁵¹ and thus a detailed and deep understanding of the mechanisms and effects of confinement on the properties of the fluids adsorbed in the porous media is essential for their development.

Over the past several decades substantial progress in the theoretical description of fluids adsorbed in disordered porous media has been achieved.^52–57 Most of the theoretical studies are based on the application of the approach developed by Madden and Glandt^58,59 and Given and Stell.⁶⁰ In these studies a porous medium is modeled as a quenched disordered configuration of the particles and the properties of the fluid adsorbed in the medium are described using the so-called replica Ornstein–Zernike (ROZ) theory.^58–60 Subsequently the theory was extended and applied to the study of the properties of a number of different fluids confined in the porous media, including site–site fluids,^61–63 associating fluids,^64–69 fluids with Coulomb interactions^70–74 and inhomogeneous systems.^75,76 The ROZ theory appears to be very useful and is able to provide relatively accurate predictions for the structure and thermodynamic properties of the fluid adsorbed in disordered porous media. However, the solution of the ROZ equation requires application of numerical methods, i.e. none of the ROZ closures proposed so far are amenable for an analytical solution. In addition, straightforward application of the ROZ approach for phase equilibrium calculations is limited by the absence of a convergent solution to the ROZ equation in a relatively large region of thermodynamic states. Recently, a theoretical approach enabling one to analytically describe the properties of hard-sphere fluid adsorbed in the hard-sphere matrix, formed by the fluid of hard spheres quenched at equilibrium, has been developed.^54,77–81 The theory is based on the appropriate extension of the scaled particle theory.⁸² The hard-sphere fluid is routinely used as a reference system in a number of thermodynamic perturbation theories. In particular, a fluid of hard spheres confined in the hard-sphere matrix was used as a reference system in the framework of the Barker–Henderson perturbation theory,^83,84 high-temperature approximation,⁸⁵ Wertheim's thermodynamic perturbation theory (TPT),^81,86–89 associative mean spherical approximation^90–93 and the collective variable method.^94–96

In the current paper we propose a simple model for functionalized disordered porous media and study the effects of confinement on the clusterization, percolation and phase behavior of a fluid of patchy particles. The media are represented by the disordered matrix of the hard-sphere obstacles bearing a certain number of off-center square-well sites (patches). Due to a strong, short-range attraction between the patches, the obstacles can bond to the fluid particles, and the fluid particles have the ability to form a network. The theoretical model is constructed by combining SPT for the porous media^54,79–81 and Wertheim's TPT for the associating fluids.^97–99 To assess the accuracy of the theory proposed we generated a set of computer simulation data and compared them against our theoretical predictions.

This paper is organized as follows. In section 2 we describe the interaction potential model and present the theory. In section 3 we briefly describe details of the computer simulations. Section 4 contains results and their discussion. Our conclusions are provided in section 5.

2. The model and theory

Patchy colloids are modeled as a one-component hard-sphere fluid with the particles decorated by the set of n₁ square-well bonding sites (patches) located on the surface. The fluid is confined in the porous media represented by the matrix of the patchy hard-sphere obstacles, quenched at hard-sphere fluid equilibrium. Each of the obstacles has on its surface n₀ square-well bonding sites (patches). The pair potential U_ij(12), which describes the interaction between the colloidal particles and between the colloidal particles and obstacles is given bywhere

1 and 2 stand for the position and orientation of the particles 1 and 2, the indices i and j take the values (i, j) = (1, 1), (1.0), and (0, 1) and denote either fluid particles (i = 1) or matrix obstacles (i = 0), the indices K and L denote patches and take the values A, B, C, …. U^(hs)_ij(r) is the hard-sphere potential between the particles of the sizes σ_i and σ_j, z₁₂ is the distance between the corresponding sites, and ε_ij and ω_ij are the depth and width of the square-well potential. Note that ε_ij and ω_ij are independent of the type of patches.

The properties of the model are calculated by combining Wertheim's thermodynamic perturbation theory (TPT) for associating fluids^97–99 and the scaled particle theory (SPT), extended to describe a hard-sphere fluid adsorbed in the hard-sphere matrix.^54,79,80 According to Wertheim's TPT Helmholtz free energy of the system A can be written in the following form

A = A_ref + ΔA_as, (3)

where A_ref is the Helmholtz free energy of the reference system, andwhere ρ_i is the number density of the particles of the type i, β = 1/(k_BT), k_B is the Boltzmann constant, T is the temperature, and X_i = X_iK is the fraction of the particles with the attractive site of type K, which belongs to the particles of the type i, non-bonded. Note that all the patches, which belong to the particles of the type i, are equivalent (2). Here X_i is obtained from the solution of the following set of equations^81,97,99and

4πX₀ρ₁σ₀₁³n₁X₁T₀₁^(as)g^(ref)₀₁ + X₀ − 1 = 0, (6)

where σ_ij = (σ_i + σ_j)/2, g^(ref)_ij is the contact value of the radial distribution function of the reference system,and f̃_ij^(as)(r) is the orientationally averaged Mayer function for the square-well site–site potential

The chemical potential, μ, and pressure, P, of the system are calculated using standard thermodynamic relations, i.e.

μ = μ_ref + Δμ_as, P = P_ref + ΔP_as (9)

and are used for the phase equilibrium calculations. Here

Coexisting densities of the low-density (ρ^(l)₁) and high-density (ρ^(h)₁) phases at a certain temperature T follow from the solution of the set of two equations:

μ(T,ρ^(l)₁) = μ(T,ρ^(h)₁), P(T,ρ^(l)₁) = P(T,ρ^(h)₁), (11)

which represents two-phase equilibrium conditions.

The reference system is represented by the hard-sphere fluid confined in the hard-sphere matrix. The properties of this reference system are calculated using corresponding extension of the scaled particle theory.^54,79–81 Closed form analytical expressions for the Helmholtz free energy of the reference system A_ref and the chemical potential μ_ref are presented in the Appendix.

3. Details of computer simulations

Monte Carlo (MC) computer simulations were performed for the model of the patchy particles confined in a patchy hard-sphere matrix, according to the model presented in section 2. We consider the hard-sphere particles of fluid and matrix with diameters of σ₀ and σ₁, respectively. We study the particles with different numbers of patches, which are symmetrically arranged on the surface of particles at a distance of half the diameter from its center. As shown in Fig. 1, the patches in a two-patch particle are located diametrically opposed to each other, in a three-patch particle they form an equilateral triangle, and in four-patch model, the patches are in vertices of a regular tetrahedron. Monte Carlo simulations have been carried out using the conventional Metropolis algorithm in the canonical ensemble (NVT)¹⁰⁰ for a two-component mixture of a patchy matrix and fluid particles taking into account both the translational and rotational moves of the fluid particles, while for the matrix particles, all degrees of freedom have been frozen. To handle the orientations and rotations of the patchy particles the quaternion representation¹⁰⁰ is employed.

Fig. 1 Models of a particle with one, two, three and four patches.

3.1. Calculation of the fraction of m-times bonded fluid

The simulations have been performed in a cubic box of size L with the periodic boundary condition applied. The size of the box has been chosen to be equal to L = 24σ₁, where σ₁ is a diameter of the fluid particles. Each simulation run has been started from a random configuration of fluid particles, both in their positions and orientations. Configurations of the matrix particles were random as well; however, they were the same for each of the systems studied at the given matrix density ρ₀ and for each set of fluid parameters. The maximum trial displacement at each translational move of a fluid particle is limited to the distance of 0.1σ₁. Each rotational move of a particle is performed around a random axis within the maximum rotation angle of 0.1 radians. The trial move ratios were taken to be 50% for translations and 50% for the rotations.

The width ω_ij of the square-well potential is assumed to be 0.119σ₁ for all patch–patch interactions ensuring formation of only one bond per patch. The depth ε_ij is also taken to be the same for all the square-well potentials, assuming that the matrix particles are functionalized with the same groups as the fluid particles, or at least they lead to the association between the fluid and the matrix particles with the energy close to that observed between the fluid particles. In Fig. 2 some snapshots with examples of the simulated systems are shown.

Fig. 2 Models of a patchy fluid in a patchy matrix: one-patch fluid in four-patch matrix L1M4 (left panel) and three-patch fluid in one-patch matrix L3M1 (right panel). The fluid particles are represented by the gray spheres and the matrix obstacles are shown as green spheres.

Each simulation process was conducted in two stages, equilibration and production runs, and both runs took 10⁶ simulation steps, where the simulation step consisted of N trial translational and rotational moves of the fluid particles. It appeared to be sufficient for producing reliable data, which were obtained for two values of the packing fractions of the matrix particles η₀ = πρ₀σ₀³/6, i.e. η₀ = 0.1 and 0.2 and different number densities of the fluid particles ranging from up to , where N₁ is the number of fluid particles and V = L³ is the volume of the simulation box. The packing fraction of the matrix is an inverse characteristic to the porosity φ₀ = 1 − η₀, which is the basic quantity describing a porous medium, ρ₀ = N₀/V, N₀ is the number of matrix particles and the diameter σ₀ = 1.5σ₁ was set larger than that of the fluid particles. In all simulations, the value of the temperature T* = k_Bt/ε₁₁ was T* = 0.12, surpassing the critical temperature of the present model fluids studied in bulk.¹⁰¹ This ensures that the critical temperature of these model fluids, being confined in the matrix, is lower than T* = 0.12.⁸¹ To accelerate the simulations the linked-cells algorithm was used with a cell size of 2.0 σ₁.¹⁰⁰ During the production run a number of bonds between the fluid and matrix particles were collected at each 10 simulation step and averaged along the simulation in order to calculate the fraction of m-times bonded fluid x₁^(m) and matrix x₀^(m) particles and to compare them with our theoretical predictions.

3.2. Estimation of the percolation threshold

The percolation threshold line for the system of patchy particles in a matrix is estimated using MC simulation in an NVT ensemble in the temperature range T* = 0.085–0.12. For these temperatures the associative interaction between patches becomes significant, and the formation of particle clusters of different sizes can appear. Due to a low probability of dissociation events, it is difficult to accurately sample a configurational space for systems using the conventional sets of MC trial moves. Thus, we extended it by an additional type of move consisting of the removal of a particle and its immediate random insertion into the simulation box. Consequently, a particle associated with one cluster can “jump” to the other cluster, if it is energetically favorable. The modified MC algorithm improves the ability to sample the configuration space of the system in the given range of temperatures, which are assumed still higher than the critical temperature. For lower temperatures it is expedient to use schemes employing biased MC moves.^100,102 The rule of the acceptance of the jump trial moves in our scheme was the same as that of the translational or rotational moves. The trial move ratios were determined to be 25% for translation, 50% for rotations, and 25% for jumps. Since the acceptance rate appeared to be small we increased the number of steps in the equilibration and production runs by up to 10⁷. The cluster analysis was performed based on the system configurations stored at each 1000th step of the production run. The depth-first search algorithm was employed to partition the patchy particles into different clusters. Two colloidal particles are assumed to be connected if the distance between their patches is less than ω₁₁. A cluster is defined as an uninterrupted sequence of interconnected particles. We assumed that for a given configuration the state of percolation was achieved if a cluster of the infinite size under the periodic boundary conditions was formed. We calculated the temperature and density plane using a step of Δ_T = 0.005 and Δ_ρ = 0.01, respectively, and the density coordinates of the percolation threshold line at a certain temperature was calculated as . Here and and are the densities with less and more than 50% of the configurations percolating, respectively.

4. Results and discussion

To systematically access accuracy of the theory we studied the properties of several versions of the model. We considered the models with different numbers of the attracting sites and different values of the density of the fluid and matrix particles. These versions are denoted as L_iM_j, where i and j are the number of sites on the fluid and matrix particles. We also considered two values for the matrix packing fraction and one value of the temperature, i.e. η₀ = 0.1 and 0.2 and T* = k_BT/ε₁₁ = 0.12. The other parameters of the model are: σ₀ = 1.5σ₁, ω₁₁ = ω₀₁ = 0.119σ₁ and ε₀₁ = ε₁₁.

In Fig. 3–6 we compared the theoretical and computer simulation results for the fraction of m times bonded fluid x₁^(m) and matrix x₀^(m) particles as a function of the fluid density ρ₁. The latter quantity was calculated using the following equation⁹⁹

Fig. 3 A fraction of the m times bonded fluid particles x₁^(m) as a function of the fluid density of the models L2M0 (top panel), L3M0 (intermediate panel), and L4M0 (bottom panel) at different values of η₀, i.e. η₀ = 0.1 (dashed lines and open circles) and η₀ = 0.2 (solid lines and filled circles). Here m = 1 (red curves and symbols), m = 2 (blue lines and symbols), m = 3 (black lines and symbols), and m = 4 (green line and symbols). The theoretical results are represented by the lines and the computer simulation results are represented by the symbols.

Fig. 4 A fraction of the m times bonded fluid (top panel) and the matrix (bottom panel) particles x_i^(m) as a function of the fluid density of the model L3M1 at different values of η₀, i.e. η₀ = 0.1 (dashed lines and open circles) and η₀ = 0.2 (solid lines and filled circles). Here m = 1 (red and brown curves and symbols), m = 2 (blue lines and symbols) and m = 3 (black lines and symbols). The theoretical results are represented by the lines and the computer simulation results by the symbols.

Fig. 5 Liquid–gas phase diagram (solid lines) and percolation thresholds (dashed lines, squares and triangles) of the model L3M1 with ω₀₁ = 0.1σ₁ and ε₀₁ = 0 (black lines and squares), ε₀₁ = 0.77ε₁₁ (blue lines) and ε₀₁ = 0.825ε₁₁ (red lines and triangles). Here lines and circles (which denote the position of the corresponding critical points) represent the results of the theory and the squares and triangles represent the computer simulation results.

Fig. 6 A fraction of the m times bonded fluid and matrix particles x_i^(m) as a function of the fluid density of the models L1M1 (top panel), L1M2 (intermediate panel), and L1M4 (bottom panel) at different values of η₀, i.e. η₀ = 0.1 (dashed lines and open circles) and η₀ = 0.2 (solid lines and filled circles). Here a fraction of the singly bonded fluid particles (m = 1) are shown by the violet lines and symbols and for the matrix particles we have: m = 1 (red curves and symbols), m = 2 (blue lines and symbols), m = 3 (black lines and symbols), and m = 4 (green line and symbols). The theoretical results are represented by the lines and the computer simulation results are represented by the symbols.

Fig. 3 shows our results for the models with matrix particles without patches and fluid particles with 2, 3 and 4 patches, i.e. models L2M0, L3P0, and L4P0. The agreement between the theoretical and the computer simulation predictions for the models L2M0 and L3M0 at both values of the matrix packing fraction, η₀ = 0.1 and 0.2, was very good. For the model L4M0 and the lower value of η₀ = 0.1 the results of the theory were less accurate and the agreement with the exact computer simulation results was only semi-quantitative. However with the increase of the matrix packing fraction η₀ = 0.2 the agreement between the theory and the computer simulation became almost as good as in the case of the other models. Our results for the models with patches on both fluid and matrix particles are presented in Fig. 4 and 6. In Fig. 4 we compare our theoretical results against the computer simulation results for the model L3M1. Here one can observe very good agreement between the theory and the simulation for x₁⁽¹⁾x₁⁽²⁾ and x₁⁽³⁾ for both values of the matrix packing fraction η₀. The theoretical predictions for the fraction of singly bonded matrix particles x₀⁽¹⁾ are accurate for the lower values of η₀ = 0.1. At η₀ = 0.2 the theoretical results are less accurate. While computer simulations predicted a slight decrease in x₀⁽¹⁾ upon the increase of η₀, theoretical calculations showed that it increased. As η₀ increased the average distance between the matrix particles decreased and the number of obstacles with the patch blocked by the nearest neighboring obstacles increased. This effect, which is not accounted for in the framework of the present version of the TPT, caused a slight decrease of x₀⁽¹⁾. On the other hand this effect was less pronounced for the models L1M1, L1M2 and L1M4 (Fig. 6) and the corresponding theoretical results were in very good agreement with the computer simulation results. This good agreement was observed in the fraction of the m times bonded fluid and matrix particles for both values of the matrix packing fraction η₀. Thus with the increase of the matrix packing fraction η₀ and the number of the patches on the fluid particles n₁ the first order TPT became less accurate. Further improvement of the theory can be achieved using higher-order versions of the TPT.

In Fig. 5 and 7 we present a liquid–gas phase diagram and percolation threshold lines for the model L3M1 with different values for the ratio between the strength of attractive fluid–fluid and fluid–matrix interactions. In addition in Fig. 5 we show our computer simulation results for the percolation lines. These lines separate the ρ vs. T plane into percolating and non-percolating regions, i.e. in the latter region the particles form finite size clusters and in the former these clusters form an infinite network. Our calculation of the percolation threshold lines is based on the extension of Flory–Stockmayer (FS) theory.^{101,103–105} According to earlier studies^104,106 percolation threshold lines for the model at hand are defined by the equality p₁₁ = 1/2, where p₁₁ is the probability of forming the bonds between the fluid particles, i.e.

Fig. 7 A fraction of the m time bonded fluid and matrix particles x_i^(m) along coexisting lines for the model L3M1 with ω₀₁ = 0.1σ₁ and ε₀₁ = 0 (top panel), ε₀₁ = 0.77ε₁₁ (intermediate panel) and ε₁₁ = 0.825ε₁₁ (bottom panel). x₁⁽⁰⁾ (black lines), x₁⁽³⁾ (blue lines) and x₀⁽⁰⁾ (red lines). Here dashed and solid lines represent gas and liquid branches of the phase diagrams, respectively.

We considered the following values for the hard-sphere size of the matrix obstacles, i.e. σ₀ = 1.388σ₁ and for the width and depth of the potential well (2): ω₀₁ = 0.1σ₁ and ε₀₁ = 0, 0.77ε₁₁ and 0.825ε₁₁. Our computer simulations for the percolation threshold lines were carried out for the lowest and highest values of the potential well depth, i.e. ε₀₁ = 0 and ε₀₁ = 0.825ε₁₁. This model was studied in the simulation box of size L = 16σ₁ using the scheme described in subsection 2 All the rest of the potential model parameters were the same as those used before. This choice of parameters allowed us to study the effects due to the competition between the formation of the bonds connecting the fluid particles, and the fluid particles and the matrix particles. While the formation of the 3-dimensional network of bonds connecting the fluid particles caused gas–liquid phase separation, the bonding of the fluid and matrix particles suppressed it. For the model with ε₀₁ = 0 the phase diagram was of the usual form, i.e. the difference between the coexisting densities of the gas and liquid phases increased on decreasing temperature (Fig. 5). Here, since the density of the gas phase ρ_{1, gas} lower than the density of the liquid phase ρ_{1, liq} the fraction of the free (or 3 time bonded) particles in the gas phase x⁽⁰⁾_{1, gas} (x⁽³⁾_{1, gas}) was larger (lower) than the corresponding fractions in the liquid phase (Fig. 7, top panel). As one would expect of the fractions of free matrix particles we had x⁽⁰⁾_{0, gas} = x⁽⁰⁾_{0, liq} = 1, where x⁽⁰⁾_{0, gas} and x⁽⁰⁾_{0, liq} denoted fractions of free (nonbonded) particles of the matrix, which confined the fluid particles either in the gas phase or in the liquid phase, respectively. An increase of ε₀₁ caused changes in the phase diagram topology. For the model with ε₀₁ = 0.77ε₁₁ the phase diagram exhibited re-entrant phase behavior, i.e. the decrease of the temperature in the range of approximately 0.065 ≥ T* ≥ 0.046 resulted in the reduction of the difference in the coexisting gas and liquid densities (Fig. 5). This effect was caused by the increased bonding of the fluid particles and matrix obstacles, which destroyed the three-dimensional network formed by the fluid particles and suppressed the phase separation. In this range of temperatures the fraction of free obstacles x⁽⁰⁾₀ in both phases rapidly decreased (Fig. 7, intermediate panel). At the same time, while the fractions of the free and three time bonded fluid particles in the liquid phase did not change much, in the gas phase they approached the liquid phase values (Fig. 7, intermediate panel). This behavior reflects the process of breaking the bonds between the fluid particles and their substitution by the bonds formed between the fluid particles and obstacles of the matrix. Upon reaching a sufficient bonding degree of the matrix obstacles this process slowed down and with further decrease of the temperature (T* < 0.046), the three-dimensional network formed by the fluid particles became stronger. As a result the difference between the coexisting gas and liquid densities increased with the temperature decrease (Fig. 5). Similar behavior can be observed of the model with ε₀₁ = 1.835ε₁₁. However, stronger fluid-matrix attraction completely suppressed the phase separation in the range of the temperature 0.060 ≤ T* ≤ 0.0545. As a result we had the phase diagram with two separate regions of the liquid–gas phase coexistence: the region at higher temperatures had the upper and lower critical points (Fig. 5 and 7) and the region at lower temperatures had the usual liquid–gas type of the critical point. These critical points and the liquid branches of the phase diagrams are located inside the percolating region (Fig. 5). The corresponding percolation line consists of two branches, i.e. high-temperature branch and low-temperature branch. The former branch is located at higher temperatures with the cluster fluid phase above and the percolating fluid phase below the percolation threshold line. The latter branch is located at lower temperatures and connects the high-temperature and low-temperature regions of the phase diagram. Here the percolating fluid phase occurs above the percolation threshold line and the cluster fluid phase occurs below the percolation threshold line. For the model with ε₀₁ = 0 the percolation threshold line was of the usual shape and the increase of ε₀₁ (ε₀₁ = 0077, 0.825) caused only a slight shift in the direction of lower temperature and higher density. The corresponding computer simulation results demonstrated similar qualitative behavior; however, the quantitative agreement between theory and simulation was rather poor: the theoretical predictions for the temperature of percolation transition at lower densities was about 15–20% higher than those of the simulation. According to previous studies the predictions of the FS approach for the percolation threshold line of the patchy colloids in the bulk conditions were relatively accurate.¹⁰⁷ The observed quantitative disagreement between theory and simulation was due to the inability of the present version of the FS approach to account for the effects of heterogeneity caused by the presence of the random porous media. Finally, we note that recently the phase behavior of the patchy colloids adsorbed in the attractive random porous medium has been studied.^88,89 In these studies the medium was represented as a matrix of hard-sphere obstacles with the attractive Yukawa interaction between the centers of the fluid and matrix particles. Similar to the present study, confinement caused the system to undergo re-entrant liquid–gas phase separation. However competition between the Yukawa interaction and bonding did not cause the phase diagram to be separated into two (or more) regions. Taking into account the good agreement between the theoretical and computer simulation predictions for the fractions X₀ and X₁ and the fact that these are the only quantities that enter into expressions for the Helmholtz free energy ΔA_as (eqn (4)) and the probability of forming the bonds p₁₁ (eqn (13)) there are good reasons to believe that our theory is accurate enough to at least qualitatively predict the peculiarities of the phase behavior discussed above. In a subsequent paper we are planning to systematically investigate the phase behavior and percolation properties of the model using an extension of the aggregation-volume-bias MC method^108,109 and the theory developed herein.

5. Conclusions

In this work we proposed a simple model of functionalized disordered porous media represented by a matrix of hard-sphere fluid particles quenched at equilibrium and decorated by a certain number of off-center square-well sites. The model was used to study the effects of confinement on the clusterization, percolation and phase behavior of the fluid of the patchy particles adsorbed in such porous media. The study was carried out by combining Wertheim's TPT for associating fluids, SPT for the porous media and Flory–Stockmayer theory of polymerization. A set of computer simulation data has been generated and used to assess the accuracy of the theory. Very good agreement between theoretical and computer simulation predictions was observed for the fractions of the m-times bonded fluid and matrix particles almost in all the cases studied. The results of the fraction of the m-times bonded matrix particles for the model with fluid particles bearing more than one patch were slightly less accurate. Theoretical predictions for the percolation threshold lines were only qualitatively accurate. Further improvement of the theory can be achieved using higher-order versions of the TPT and correcting the FS approach to account for the heterogeneity of the porous media.

The liquid–gas phase diagram and percolation threshold line for the model with three-patch fluid particles and one-patch matrix particles were calculated and analyzed. It is demonstrated that confinement can substantially change the shape of the phase diagram. Bonding between the fluid particles causes the formation of a network and bonding of the fluid particles to the obstacles suppresses this process. Competition between these two effects define the shape of the phase diagram and gives rise to a re-entrant phase behavior with three critical points and two separate regions of the liquid–gas phase coexistence.

Conflicts of interest

There are no conflicts to declare.

Appendix: Expressions for Helmholtz free energy, chemical potential and pressure of the hard-sphere fluid confined in the hard-sphere matrix

The reference system is represented by the hard-sphere fluid confined in the hard-sphere matrix. Detailed description of the properties of such system using SPT was carried out earlier.^54,79–81 Based on these studies we have the following expressions for the Helmholtz free energy ΔA_hs and for the contact values of the radial distribution functions g^(hs)₁₀(r) and g^(hs)₁₁(r) for the hard-sphere fluid adsorbed in the hard-sphere matrix, i.e.andwhere N_hs is the number of hard-sphere particles in the fluid,

η ₀ = πρ₀σ₀³/6, φ₀ = 1−η₀, η₁ = πρ₁σ₁³/6 and ϕ = exp(−βμ₁^(ex)).

Here

Z ₀ = (1 + η₀ + η₀²)/(1 − η₀)³, τ_1j = σ₁/σ_j and j = 0, 1.

Acknowledgements

YVK acknowledges the financial support through the MSCA4Ukraine project, which is funded by the European Union, and TP and MH gratefully acknowledge the financial support from the National Research Foundation of Ukraine (project no. 2020.02/0317). YVK expresses his gratitude to the Vanderbilt University, where part of this research was conducted, for the hospitality.

References

1. A. Corma, Chem. Rev., 1997, 97, 2373–2420.
2. A. Vinu, T. Mori and K. Ariga, Sci. Technol. Adv. Mater., 2006, 7, 753–771.
3. I. Slowing, B. Trewyn, S. Giri and V.-Y. Lin, Adv. Funct. Mater., 2007, 17, 1225–1236.
4. Handbook of Heterogeneous Catalysis, ed. G. Ertl, H. Knözinger, F. Schüth and J. Weitkamp, Wiley, 2008.
5. S. P. Adiga, C. Jin, L. A. Curtiss, N. A. Monteiro-Riviere and R. J. Narayan, Wiley Interdiscip. Rev.: Nanomed. Nanobiotechnol., 2009, 1, 568–581.
6. A. Corma, H. García and F. X. L. i Xamena, Chem. Rev., 2010, 110, 4606–4655.
7. P. Stroeve and N. Ileri, Trends Biotechnol., 2011, 29, 259–266.
8. P. Misaelides, Microporous Mesoporous Mater., 2011, 144, 15–18.
9. R. Dawson, A. I. Cooper and D. J. Adams, Prog. Polym. Sci., 2012, 37, 530–563.
10. J. Rouquerol, F. Rouquerol, P. Llewellyn, G. Maurin and K. S. Sing, Adsorption by powders and porous solids: principles, methodology and applications, Academic press, 2013.
11. J. Liu, L. Chen, H. Cui, J. Zhang, L. Zhang and C.-Y. Su, Chem. Soc. Rev., 2014, 43, 6011–6061.
12. T. Zhang and W. Lin, Chem. Soc. Rev., 2014, 43, 5982–5993.
13. L. Canham, Handbook of porous silicon, Springer International Publishing Berlin, Germany, 2014.
14. N. Pal and A. Bhaumik, RSC Adv., 2015, 5, 24363–24391.
15. Y. Li, H. Xu, S. Ouyang and J. Ye, Phys. Chem. Chem. Phys., 2016, 18, 7563–7572.
16. J. Liang, Z. Liang, R. Zou and Y. Zhao, Adv. Mater., 2017, 29, 1701139.
17. T. Kumeria, S. J. P. McInnes, S. Maher and A. Santos, Expert Opin. Drug Delivery, 2017, 14, 1407–1422.
18. G. Danda and M. Drndić, Curr. Opin. Biotechnol., 2019, 55, 124–133.
19. Q. Wang and D. Astruc, Chem. Rev., 2020, 120, 1438–1511.
20. L. Miao, Z. Song, D. Zhu, L. Li, L. Gan and M. Liu, Mater. Adv., 2020, 1, 945–966.
21. Y. Li and J. Yu, Nat. Rev. Mater., 2021, 6, 1156–1174.
22. N. Hosono and T. Uemura, Acc. Chem. Res., 2021, 54, 3593–3603.
23. R. Moretta, L. D. Stefano, M. Terracciano and I. Rea, Sensors, 2021, 21, 1336.
24. G. Singh, R. Bahadur, J. M. Lee, I. Y. Kim, A. M. Ruban, J. M. Davidraj, D. Semit, A. Karakoti, A. H. A. Muhtaseb and A. Vinu, Chem. Eng. J., 2021, 406, 126787.
25. R. Saha, B. Mondal and P. S. Mukherjee, Chem. Rev., 2022, 122, 12244–12307.
26. M. Hadden, D. Martinez-Martin, K.-T. Yong, Y. Ramaswamy and G. Singh, Materials, 2022, 15, 2111.
27. B. Yue, S. Liu, Y. Chai, G. Wu, N. Guan and L. Li, J. Energy Chem., 2022, 71, 288–303.
28. A. Popat, S. B. Hartono, F. Stahr, J. Liu, S. Z. Qiao and G. Q. M. Lu, Nanoscale, 2011, 3, 2801.
29. R. M. Milton, in Zeolite Synthesis, ed. M. L. Occelli and H. E. Robson, American Chemical Society, 1989, ch. 1, pp. 1–10.
30. S. Wang and Y. Peng, Chem. Eng. J., 2010, 156, 11–24.
31. M. Xia, Z. Chen, Y. Li, C. Li, N. M. Ahmad, W. A. Cheema and S. Zhu, RSC Adv., 2019, 9, 20941–20953.
32. W. Gu and G. Yushin, Wiley Interdiscip. Rev.: Energy Environ., 2013, 3, 424–473.
33. S. S. Kistler, Nature, 1931, 127, 741–741.
34. C. T. Kresge, M. E. Leonowicz, W. J. Roth, J. C. Vartuli and J. S. Beck, Nature, 1992, 359, 710–712.
35. D. Zhao, J. Feng, Q. Huo, N. Melosh, G. H. Fredrickson, B. F. Chmelka and G. D. Stucky, Science, 1998, 279, 548–552.
36. H. Li, M. Eddaoudi, M. O'Keeffe and O. M. Yaghi, Nature, 1999, 402, 276–279.
37. M. Wen, G. Li, H. Liu, J. Chen, T. An and H. Yamashita, Environ. Sci.: Nano, 2019, 6, 1006–1025.
38. T. L. Easun, F. Moreau, Y. Yan, S. Yang and M. Schröder, Chem. Soc. Rev., 2017, 46, 239–274.
39. J.-R. Li, R. J. Kuppler and H.-C. Zhou, Chem. Soc. Rev., 2009, 38, 1477.
40. X. Zhang, Z. Chen, X. Liu, S. L. Hanna, X. Wang, R. Taheri-Ledari, A. Maleki, P. Li and O. K. Farha, Chem. Soc. Rev., 2020, 49, 7406–7427.
41. N. Mizoshita, T. Tani and S. Inagaki, Chem. Soc. Rev., 2011, 40, 789–800.
42. A. Modak, J. Mondal, V. K. Aswal and A. Bhaumik, J. Mater. Chem., 2010, 20, 8099.
43. N. B. McKeown, B. Gahnem, K. J. Msayib, P. M. Budd, C. E. Tattershall, K. Mahmood, S. Tan, D. Book, H. W. Langmi and A. Walton, Angew. Chem., Int. Ed., 2006, 45, 1804–1807.
44. D. Wu, F. Xu, B. Sun, R. Fu, H. He and K. Matyjaszewski, Chem. Rev., 2012, 112, 3959–4015.
45. D. Chandra and A. Bhaumik, J. Mater. Chem., 2009, 19, 1901.
46. H. Ou, W. Zhang, X. Yang, Q. Cheng, G. Liao, H. Xia and D. Wang, Environ. Sci.: Nano, 2018, 5, 169–182.
47. Y. Xu, S. Jin, H. Xu, A. Nagai and D. Jiang, Chem. Soc. Rev., 2013, 42, 8012.
48. S. Das, P. Heasman, T. Ben and S. Qiu, Chem. Rev., 2017, 117, 1515–1563.
49. S.-Y. Ding and W. Wang, Chem. Soc. Rev., 2013, 42, 548–568.
50. K. Geng, T. He, R. Liu, S. Dalapati, K. T. Tan, Z. Li, S. Tao, Y. Gong, Q. Jiang and D. Jiang, Chem. Rev., 2020, 120, 8814–8933.
51. R. Yadav, T. Baskaran, A. Kaiprathu, M. Ahmed, S. V. Bhosale, S. Joseph, A. H. Al-Muhtaseb, G. Singh, A. Sakthivel and A. Vinu, Chem. – Asian J., 2020, 15, 2588–2621.
52. M. L. Rosinberg, in New Approaches to Problems in Liquid State Theory. NATO Science Series, ed. C. Caccamo, J. Hansen and G. Stell, Springer, Dordrecht, 1999, vol. 529, pp. 245–278.
53. L. Sarkisov and P. R. V. Tassel, J. Phys.: Condens. Matter, 2008, 20, 333101.
54. M. Holovko, T. Patsahan and W. Dong, Pure Appl. Chem., 2013, 85, 115–133.
55. M. Holovko, V. Shmotolokha and T. Patsahan, in Physics of Liquid Matter: Modern Problems. Springer Proceedings in Physics, ed. L. Bulavin and N. Lebovka, Springer, Cham, 2015, vol. 171, ch. 1, pp. 3–30.
56. W. Dong and X. Chen, Sci. China: Phys., Mech. Astron., 2018, 61, 70501.
57. O. Pizio, in Computational Methods in Surface and Colloid Science, ed. M. Borowko, CRC Press, United States, 2019, ch. 6, pp. 293–345.
58. W. G. Madden and E. D. Glandt, J. Stat. Phys., 1988, 51, 537–558.
59. W. G. Madden, J. Chem. Phys., 1992, 96, 5422–5432.
60. J. A. Given and G. Stell, J. Chem. Phys., 1992, 97, 4573–4574.
61. A. Kovalenko and F. Hirata, J. Chem. Phys., 2001, 115, 8620–8633.
62. D. Chandler, J. Phys.: Condens. Matter, 1991, 3, F1–F8.
63. A. P. Thompson and E. D. Glandt, J. Chem. Phys., 1993, 99, 8325–8329.
64. A. Trokhymchuk, O. Pizio, M. Holovko and S. Sokolowski, J. Chem. Phys., 1997, 106, 200–209.
65. G. A. Orozco, O. Pizio, S. Sokolowski and A. Trokhymchuk, Mol. Phys., 1997, 91, 625–634.
66. O. Pizio, Y. Duda, A. Trokhymchuk and S. Sokolowski, J. Mol. Liq., 1998, 76, 183–194.
67. P. Padilla, O. Pizio, A. Trokhymchuk and C. Vega, J. Phys. Chem. B, 1998, 102, 3012–3017.
68. B. M. Malo, O. Pizio, A. Trokhymchuk and Y. Duda, J. Colloid Interface Sci., 1999, 211, 387–394.
69. T. Urbic, V. Vlachy, O. Pizio and K. Dill, J. Mol. Liq., 2004, 112, 71–80.
70. B. Hribar, O. Pizio, A. Trokhymchuk and V. Vlachy, J. Chem. Phys., 1997, 107, 6335–6341.
71. B. Hribar, O. Pizio, A. Trokhymchuk and V. Vlachy, J. Chem. Phys., 1998, 109, 2480–2489.
72. B. Hribar, V. Vlachy, A. Trokhymchuk and O. Pizio, J. Phys. Chem. B, 1999, 103, 5361–5369.
73. V. Vlachy, H. Dominguez and O. Pizio, J. Phys. Chem. B, 2004, 108, 1046–1055.
74. B. Hribar-Lee, M. Lukšič and V. Vlachy, Annu. Rep. Prog. Chem., Sect. C: Phys. Chem., 2011, 107, 14.
75. O. Pizio and S. Sokolowski, Phys. Rev. E: Stat. Phys., Plasmas, Fluids, Relat. Interdiscip. Top., 1997, 56, R63–R66.
76. A. Kovalenko, S. Sokołowski, D. Henderson and O. Pizio, Phys. Rev. E: Stat. Phys., Plasmas, Fluids, Relat. Interdiscip. Top., 1998, 57, 1824–1831.
77. M. Holovko and W. Dong, J. Phys. Chem. B, 2009, 113, 6360–6365.
78. W. Chen, W. Dong, M. Holovko and X. S. Chen, J. Phys. Chem. B, 2010, 114, 1225–1225.
79. T. Patsahan, M. Holovko and W. Dong, J. Chem. Phys., 2011, 134, 074503.
80. M. Holovko, T. Patsahan and W. Dong, Condens. Matter Phys., 2017, 20, 33602.
81. Y. V. Kalyuzhnyi, M. Holovko, T. Patsahan and P. T. Cummings, J. Phys. Chem. Lett., 2014, 5, 4260–4264.
82. H. Reiss, H. L. Frisch and J. L. Lebowitz, J. Chem. Phys., 1959, 31, 369–380.
83. M. Holovko, T. Patsahan and V. Shmotolokha, Condens. Matter Phys., 2015, 18, 13607.
84. A. Nelson, Y. Kalyuzhnyi, T. Patsahan and C. McCabe, J. Mol. Liq., 2020, 300, 112348.
85. T. V. Hvozd and Y. V. Kalyuzhnyi, Soft Matter, 2017, 13, 1405–1412.
86. T. V. Hvozd, Y. V. Kalyuzhnyi and P. T. Cummings, J. Phys. Chem. B, 2018, 122, 5458–5465.
87. T. Hvozd, Y. V. Kalyuzhnyi and V. Vlachy, Soft Matter, 2020, 16, 8432–8443.
88. T. V. Hvozd, Y. V. Kalyuzhnyi, V. Vlachy and P. T. Cummings, J. Chem. Phys., 2022, 156, 161102.
89. T. Hvozd, Y. V. Kalyuzhnyi and V. Vlachy, Soft Matter, 2022, 18, 9108–9117.
90. M. Holovko and Y. V. Kalyuzhnyi, Mol. Phys., 1991, 73, 1145–1157.
91. H. Krienke, J. Barthel, M. Holovko, I. Protsykevich and Y. Kalyushnyi, J. Mol. Liq., 2000, 87, 191–216.
92. M. Holovko, T. Patsahan and O. Patsahan, J. Mol. Liq., 2017, 228, 215–223.
93. M. Holovko, T. Patsahan and O. Patsahan, J. Mol. Liq., 2017, 235, 53–59.
94. M. F. Holovko, O. Patsahan and T. Patsahan, J. Phys.: Condens. Matter, 2016, 28, 414003.
95. O. V. Patsahan, T. M. Patsahan and M. F. Holovko, Phys. Rev. E, 2018, 97, 022109.
96. O. Patsahan, T. Patsahan and M. Holovko, J. Mol. Liq., 2018, 270, 97–105.
97. M. S. Wertheim, J. Stat. Phys., 1986, 42, 459–476.
98. M. S. Wertheim, J. Stat. Phys., 1986, 42, 477–492.
99. M. S. Wertheim, J. Chem. Phys., 1987, 87, 7323–7331.
100. M. P. Allen and D. J. Tildesley, Computer Simulation of Liquids, Oxford University Press, 2017.
101. E. Bianchi, P. Tartaglia, E. Zaccarelli and F. Sciortino, J. Chem. Phys., 2008, 128, 144504.
102. S. Wierzchowski and D. A. Kofke, J. Chem. Phys., 2001, 114, 8752–8762.
103. E. Bianchi, P. Tartaglia, E. L. Nave and F. Sciortino, J. Phys. Chem. B, 2007, 111, 11765–11769.
104. D. de las Heras, J. M. Tavares and M. M. T. da Gama, Soft Matter, 2011, 7, 5615.
105. J. M. Tavares, P. I. C. Teixeira, M. M. T. da Gama and F. Sciortino, J. Chem. Phys., 2010, 132, 234502.
106. S. Roldán-Vargas, F. Smallenburg, W. Kob and F. Sciortino, J. Chem. Phys., 2013, 139, 244910.
107. Y. Kalyuzhnyi, C. Iacovella, H. Docherty, M. Holovko and P. Cummings, J. Stat. Phys., 2011, 145, 481–506.
108. B. Chen and J. I. Siepmann, J. Phys. Chem. B, 2001, 105, 11275–11282.
109. J. Russo, J. Tavares, P. Teixeira, M. T. Da Gama and F. Sciortino, Phys. Rev. Lett., 2011, 106, 085703.

---