Description of a polar molecular liquid in a disordered microporous material with activating chemical groups by a replica RISM theory

We develop a replica generalization of the reference interaction site model (replica RISM) integral equation theory to describe the structure and thermodynamics of a polar molecular liquid sorbed in a quenched disordered porous matrix including polar chemical groups. It provides a successful approach to realistic models of molecular liquids, and properly allows for the effect of a quenched disordered medium on the sorbed liquid. The description can be readily extended to a mobile liquid comprising a mixture of ionic and polar molecular species. The replica RISM integral equations are complemented by the HNC closure and its partial linearization (PLHNC), adequate to ionic and polar molecular liquids. In these approximations, closed expressions for the excess chemical potentials of the quenched-annealed system are derived. We extend the description to the case of soft core interaction potentials between all species of the quenched-annealed system, in which the liquid and matrix equilibrium distributions are characterized in general by two different temperatures. The replica RISM/PLHNC-HNC theory is applied to water sorbed in a quenched matrix roughly modelling porous carboneous material activated with carboxylic (-COOH) groups. The results are in qualitative agreement with experiment for water confined in disordered materials.


Introduction
Microporous materials are of significant engineering as well as scientific interest due to their enhanced adsorption, separation, and catalytic properties [1,2].They can be most generally subdivided into those possessing a crystalline order, such as zeolites, and disordered on a microscopic scale, for instance, activated carbons, silica gels, and porous glasses.The microscopic structure of the latter is formed in the process of gelation or a glass transition, and then assumed to be frozen in space.Theoretical description of liquid adsorbed in the pores of such a quenched random matrix requires specific statistical-mechanical means to allow for a disorder of the amorphous material.Although a reliable modelling can be ultimately done by molecular simulations, integral equation theory offers a very cost-efficient alternative providing fast calculation of the thermodynamic and structural properties of the liquid, which is important in industrial applications.This also allows one to access large-scale systems, as for now still unfeasible to simulations for time and memory limitations.Moreover, the theory is uniquely capable of yielding analytical and asymptotical relations for the fluid properties.
There has been a number of theoretical studies treating the binary system of an annealed, thermally equilibrated fluid confined in a porous matrix of immobile solid obstacles specified as a quenched distribution of equilibrium or nonequilibrium fluid .Madden and Glandt (MG) [3,5] performed a cluster diagram analysis of this quenched-annealed, or partially quenched system and derived the integral equations of Ornstein-Zernike type (OZ) [31].In physical essence, the MG-OZ equations can be viewed as treating the matrix as a single rigid supermolecule immersed in the fluid at infinite dilution [4].Given and Stell [8,9] showed that this set of equations was not entirely adequate to the correct topological specifications of the system and proposed a different formalism and a corrected set, called the replica OZ equations.It is based on the continuum replica method [8] developed to treat systems containing both quenched and annealed degrees of freedom, and starting with the average free energy of the mobile fluid specified as the statistical average over the matrix distribution.This average is difficult because, unlike an entirely annealed system, it requires averaging the logarithm rather than the statistical sum itself.It is alleviated by using the replica identity to replace the logarithm with a limit of powers s.Thereby, the partly quenched system is replaced with an equilibrated, entirely annealed mixture of the matrix species (now mobile) and s replicas of the fluid which interact with the matrix but not with each other.The corresponding OZ equations set for this s + 1-component mixture is then taken in the analytic continuation limit s → 0, leading to the replica OZ equations.As an important advantage, the latter explicitly include the so-called blocking (or disconnected) terms representing the correlation between fluid particles originating from different replicas (noninteracting with each other) due to the presence of immobile matrix particles.In this sense, the MG-OZ set does not distinguish between a quenched and annealed matrix, and its only difference from the equilibrium mixture OZ equations is in the quenched subsystem input, independent on the fluid subsystem [9].One can show [9,10] that the effect of the direct blocking correlations becomes increasingly significant with the density of the matrix obstacles, and is substantial even at low density of adsorbed fluid.It is especially essential for systems with charged species [17].
Thereafter, numerous studies were carried out exploring the replica OZ approach for both fluid and matrix represented with atomic liquid models [10][11][12][13][14][15][16][17][18][19][20][21][22][23] (see also recent reviews in [29,30]).In particular, the questions of the thermodynamics of a quenched-annealed system were thoroughly addressed [11][12][13][14].Different approximations for closures to the replica OZ equations were investigated [7,10], including bridge corrections extracted from a reference system [15] or derived to satisfy thermodynamic and structural consistency [16].A generalization to the inhomogeneous replica OZ theory was elaborated, and spatially inhomogeneous systems of quenched-annealed atomic fluids in spatial confinements of various geometry were studied [21].The replica OZ formalism was extended to describe the structure and the thermodynamic properties of a multicomponent liquid mixture within a porous medium [25,22].However, the case of a molecular liquid sorbed in a quenched matrix is more complicated.
A successful description for realistic molecular liquids of various complexity is provided by the reference interaction site model (RISM) which is an orientational reduction of the molecular OZ integral equation [31].The RISM method was pioneered by Chandler and Andersen [32], and then extended by Hirata and co-workers to polar and quadrupolar liquids [33] and to ions in a molecular polar solvent [34] by adapting the hypernetted chain (HNC) closure.A great advantage of the RISM approach over other integral equation theories for molecular liquids [31] is that it can easily handle the description of solution including polar as well as nonpolar polyatomic molecules and molecular ions, and to properly take into account such chemical specificities as hydrogen bonding [35].
By viewing the porous solid as one huge ubiquitous molecule and using heuristic arguments of the RISM theory, Chandler [4] obtained the RISM equations for polyatomic fluids sorbed in a quenched amorphous material which are similar to the MG-OZ equations.Madden [5] also discussed this approach from the point of view of the MG diagram analysis and suggested the use of the Chandler-Silbey-Ladanyi (CSL) equations [36,31] to provide the interaction site theory to be proper in the diagrammatic sense.Thompson and Glandt [24] combined the polymer RISM theory (PRISM) of Schweizer and Curro [37] with the MG formalism [3,5] to describe a polymer fluid of ideal freely jointed hard sphere chains in a hard sphere matrix.Ford, Thompson, and Glandt [25] further developed the MG-PRISM as well as MG-OZ methods to obtain the thermodynamic properties of fluids confined in disordered porous solids by using the scaled particle theory (SPT) approach.At the level of the MG approximation, however, the direct blocking effects essentially distinguishing a quenched disordered matrix, and important for ionic and polar molecular fluids are again neglected.
In several recent works [26,27], the replica formalism was generalized for Wertheim's multidensity OZ equations [39] describing dimerizing and polymerizing fluids of molecular particles represented by the primitive models of chemical association [40].For instance, a four-site associative model mimicking hydrogen bonding network of water [41] was used to study the phase transition and structure of a networkforming fluid sorbed in a hard sphere matrix [27].The multidensity formalism was also employed to describe the structure of quenched polymerized matrices [28].Al-though providing an advanced description of the most important generic properties of associating fluids [42], this multidensity approach suffers from the quickly arising limitations of the primitive models in representing complex inorganic and polyatomic organic species.The very recent drastic improvement [43] combines Wertheim's multidensity association theory with the molecular OZ equation [31] in the rotational invariant treatment [44].Unfortunately, the latter approach becomes increasingly cumbersome with asphericity of polyatomic molecules.
To our knowledge, so far there has not been elaborated a replica extension of the RISM equations, which would provide an improved description including direct blocking effects of quenched disorder in realistic molecular systems with polar and charged species.In this work, we develop the replica RISM integral equation theory for a polar molecular liquid confined in a quenched disordered matrix of molecular species that can also contain charged groups.At the moment, we postpone the investigation of consistent bridge corrections [16], and complement the replica RISM equations with the HNC closure or its partial linearization (PLHNC) [45][46][47] which are adequate to polar molecular systems [33,34] we intend to study.Within these approximations, we derive a closed analytical form for the excess chemical potential of the liquid, similarly to the entirely annealed molecular mixtures [48].As a simple numerical illustration, the theory is applied to water confined in a quenched matrix modelling a carboneous material with as well as without activating carboxylic groups.

Replica site-site OZ equations
We consider a mixture of two molecular species, one quenched and one annealed, denoted with indices 0 and 1, respectively.Molecular species 0 constituting a porous matrix are equilibrated separately and then frozen in place with the spatial correlations corresponding to a canonical ensemble at a temperature T 0 .A liquid of mobile molecular species 1 in the quenched external field of the matrix is in equilibrium described by a grand canonical distribution at a temperature T .In the statistical-mechanical formalism below, we employ the atomic point of view regarding molecules as clusters of atoms bound by intramolecular potentials of chemical bonds [49], although a molecular approach can equally be used.Following Rosinberg, Tarjus, and Stell [12], we introduce the grand potential of the mobile liquid in this quenched-annealed system, Ω 1 , as a statistical average of the grand potential of the liquid at a particular matrix realization q 0 over the matrix distribution, with the canonical partition function of the matrix and the matrix-dependent grand partition function of the liquid where denotes the positions of M 0 interaction sites in each of N 0 molecules comprising the frozen matrix, and are the site coordinates of N 1 molecules (each with M 1 interaction sites) of the annealed liquid, α=1 µ 1 α is the chemical potential of molecules of the liquid comprising the site contributions µ 1 α , β 0 = 1/(kT 0 ) and β = 1/(kT ) are the inverse temperatures, and H ij stands for the sum of all interactions between species i and j (including intramolecular terms for i = j).Since the system of matrix molecules was equilibrated before freezing, the matrix spatial correlations comprise just the equilibrium configurational contribution, without dynamic terms.Therefore, after scaling the matrix interaction potentials H 00 by the factor λ 0 = T /T 0 , we can treat the matrix correlations to be corresponding to the same temperature as that of the liquid.This is trivial in the case of a matrix of rigid core particles, for which the configurational part is temperature independent.For soft matrix potentials H 00 , the presence of two different temperatures somewhat modifies the thermodynamics of the quenched-annealed system but does not affect its statistical-mechanical description.Then, within the continuum generalization of the replica method [8], we can use the replica identity, to relate the partly quenched average, Ω 1 , to the analytic continuation at s = 0 of the annealed averages of the moments Ξ s 1 , For integer values of s, the annealed average Ξ rep (s) takes the form of the equilibrium partition function of a fully annealed (s + 1)-component molecular liquid mixture comprising the matrix species, now annealed, and s equivalent replicas of the liquid, where is the canonical partition function of the annealed (s + 1)-component molecular system.The matrix component is described in (7) by the canonical ensemble, whereas the s liquid replicas are described by the grand canonical partition with the same chemical potential µ 1 .As a special feature, this mixture is highly nonadditive since there is no interaction between different replicas of the liquid, H ij = 0 for i = j.
On simple rearrangement taking into account that lim s→0 Ξ rep (s)/Z 0 = 1, the averaged grand potential ( 6) is expressed through the thermodynamic function of the annealed replicated system, Ω rep (s) = −kT ln Ξ rep (s), as Then all the thermodynamic properties of the initial partially quenched molecular system are obtained similarly to [12] from those of the replicated, (s + 1)-component annealed molecular liquid by using equation ( 9) in the common assumption [8][9][10][11][12][13][14]24] that the limit s → 0 does not break the permutational symmetry of the liquid replicas.
Owing to the equivalence of ensembles verified for quenched-annealed systems in [12], the description of the liquid sorbed in the matrix can equally be performed in the canonical ensemble [9][10][11][12][13][14].In this case, one considers the average free energy of the molecular liquid, where the free energy of the ideal gas of M 1 N 1 separate noninteracting sites of liquid under averaging, and Λ α = (2πβ/m α ) 1/2 is the de Broglie thermal wavelength of site α with mass m α .Employing the replica identity (5), it is presented in the form through the free energy of the annealed replicated system with the ideal gas part A id rep (s) = A id 0 + sA id 1 comprising the contributions A id 0 from the matrix as well as A id 1 from every liquid replica.The former term disappears in the replica limit (11).
Restricting the present consideration to the case of rigid molecules, we apply Chandler and Andersen's site-site OZ, or in another terminology, RISM integral equations for site-site correlation functions [32] where h ij αγ is the intermolecular part of the total correlation function between interaction sites α and γ of species i and j, c ij αγ is the site-site direct correlation function, the terms h 12  αγ and c 12 αγ mean the site-site correlations between molecules from equivalent but different replicas of the liquid, ω ii αγ is the intramolecular matrix of species i, ρ i is the number density of molecules of species i (ρ 1 = ρ 2 ), and an asterisk * means convolution in direct space and summation over repeating site indices.For rigid molecules with site separations l ii αγ , the intramolecular matrix is specified in reciprocal space as ω ii αγ (k) = sin(kl ii αγ )/(kl ii αγ ), where ω ii αα (k) = 1 for l ii αα (k) = 0. Assuming that there is no replica symmetry breaking in the analytic continuation of equations (13) to s = 0, we obtain the set of the replica RISM integral equations, where we have identified the blocking (or disconnected) site-site correlation functions similarly to the atomic liquid case [8][9][10][11][12][13][14] as the subset of graphs with all paths passing through at least one matrix species site, h The intermolecular site-site total correlation function of the mobile liquid confined in a random matrix, h 11 αγ , is in general defined by the relation through the mean single and pair densities, ρ 1 (r 1 ) = ρ 1 (r 1 ; q 0 ) and ρ 11 αγ (r 1 , r 2 ) = ρ 11 αγ (r 1 , r 2 ; q 0 ), where the overbar denotes a disorder average over all matrix realizations q 0 .In the case of a liquid of rigid molecules, the intramolecular site-site correlation function s 11  αγ has the simple form s 11 αγ (r) = (1 − δ αγ )δ r − l 11 αγ / 4π(l 11 αγ ) 2 .For a statistically homogeneous matrix, the densities are translationally invariant and equation (15) reduces to The site-site density correlation function ( 16) can be broken up, similarly to quenched-annealed atomic liquids [12], into the connected part which is a disorder average of the spatially inhomogeneous site-site density correlation in the matrix of a particular realization q 0 , and the blocking part which is a disorder-averaged correlation between two spatially inhomogeneous single densities at a particular matrix realization q 0 , In the approximation neglecting the blocking part of the site-site direct correlation functions, c (b) αγ = 0, the replica RISM integral equations ( 14) reduce to Chandler's RISM equations for a molecular fluid in a quenched amorphous material [4].This is much like the replica OZ equations for atomic liquids reduce to the Madden-Glandt (MG) OZ theory [3,5] when neglecting the blocking effects [8][9][10][11][12][13][14].
Notice that the replica RISM equations ( 14) are readily generalized to the case of a matrix comprising a mixture of molecular species as well as of a multicomponent mobile liquid.This is done merely by extending the interaction site subscripts to include component indices a, c in the correlation functions in the form h ij aα,cγ .The densities ρ i a also acquire a component index, over which summation is performed.The intramolecular matrix is then modified as ω ii aα,cγ (k) = δ ac sin(kl ii aα,aγ )/(kl ii aα,aγ ), where l ii aα,aγ are the site-site separations in a molecule of component a of the matrix (i = 0) and liquid (i = 1).
Obviously, the replica RISM equations ( 14) inherit the inconsistencies of the RISM approach, such as for instance, inclusion of improper diagrams and imperfect representation of the excluded volume effects [31].It should be emphasized, however, that the classification and the proper treatment of the blocking and connected parts of the correlations, essential for a partially quenched system is retained in equations (14).Despite its deficiencies and owing to effective cancellation of errors, the RISM theory has proved to be successful for a number of chemical and biological systems of practical interest [35].Therefore one can expect that its replica extension is capable of a good qualitative description of molecular liquids in the partially quenched case as well.
For rigid molecular species, the replica RISM equations ( 14) can be seen as a full orientational reduction of the replica molecular OZ integral equations which are obtained in the limit s → 0 from the corresponding molecular OZ equations for the (s + 1)-component annealed molecular mixture.Such a derivation can be readily done similarly to references [45,46,50] by using the implicit assumption of the RISM theory that the molecular direct correlation function is decomposable into partial site contributions, since it has a long-range asymptotics of the molecular interaction potential generally comprising additive site-site terms.It is better, however, to introduce the RISM description before the replica identity limit, at once for the annealed mixture as described by equations (13).This seems to be more consistent in view of the possibility to further extend the replica RISM approach at hand to the case of deformable molecules by treating the intramolecular matrix ω ii αγ (r) as an intramolecular distribution function to be determined within the density functional approach [51,52].Such an approach has been discussed by Chandler [4] at the level of the MG-OZ approximation.Also promising is generalization to a replica polymer RISM (PRISM) theory with the intramolecular distribution functions of flexible polymer chains predetermined, like in the PRISM theory of Schweizer and Curro, in some reasonable form [37] or taken from simulations [38].As compared to the MG-PRISM equations elaborated by Thompson and Glandt [24], such a replica PRISM theory would provide a more consistent description for polymeric liquids sorbed in disordered porous materials, especially for species with polar or charged groups.We postpone these possibilities to further works.

Closures
As usually in integral equation theory of liquids, the replica RISM equations need to be complemented with closures.The exact relations for the site-site correlations can be formally written as where αγ (r) + 1 are as usually the site-site distribution functions, and the intermolecular interactions are specified by additive site-site potentials u ij αγ (r) dependent on separation r between interaction sites α and γ of species i and j.The blocking correlations originate from those between different liquid replicas of the annealed mixture, equation (13e), and thus do not involve an interaction potential in the closure (19b).Different closures use various approximations to the unknown bridge functions b ij αγ (r) and b (b) αγ (r).As appropriately noted by Kierlik, Rosinberg, Tarjus, and Monson [14], the validity of these approximations as well as of the replica symmetry nonbreaking remains partly uncontrollable and can be checked only by comparison of the results for some featured systems against molecular simulations.For quenched-annealed atomic fluids, the standard approximations of liquid state theory [31] have been extensively tested [7,10].As was shown, such linearized closures as the Percus-Yevick (PY) one and the mean spherical approximation (MSA) are quite adequate for quenched-annealed fluids with repulsive core or Lennard-Jones interactions and reproduce their behaviour well.However, with these closures the replica OZ equations reduce to the MG-OZ theory neglecting the direct blocking correlations [9][10][11][12][13][14].It is demonstrated [17] that for electrolyte solutions confined in a hard sphere matrix, the blocking effects become especially essential and can be recapitulated with the hypernetted chain (HNC) closure, most relevant to systems of charged species [31].Similarly to the replica OZ for atomic liquids, the replica RISM equations ( 14) with the site-site PY or MSA closures reduce to Chandler's RISM approximation at the MG level [4].On the other hand, the conventional RISM equations in the HNC approximation neglecting the site-site bridge functions proved to be adequate for modelling ion-molecular nonaqueous as well as aqueous solutions [33][34][35].Therefore, we complement the replica RISM equations with the site-site HNC closures by neglecting in (19) the bridge functions b ij αγ and b Instead of (20), bridge corrections in the closures (19) can be attempted in order to gain thermodynamic and structural consistency, similarly to those elaborated for partially quenched atomic liquids [16].This issue, however, constitutes a challenge for the RISM approach itself, and will be addressed elsewhere.
For molecular systems with strong site-site attraction, in particular, for liquid mixtures in the case of enhanced clustering, the HNC approximation can become divergent.To treat this problem, we have elaborated a partially linearized HNC closure (PLHNC) [45][46][47].For the quenched-annealed system at hand, it is written as The PLHNC approximation (21) combines the HNC exponent for density depletion regions of h < 0 and its linearization for enrichment regions of h > 0. The distribution function and its first derivative are continuous at the joint point, X = 0, by construction.The linearization in the enrichment regions prevents exponential rise of the distribution function in the regions of a strong attractive potential bringing about the divergence.For bulk ambient water [45], simple ions in aqueous solution [46], and molecular ions in a polar molecular solvent [47], the PLHNC approximation as compared to the HNC one somewhat reduces and widens high peaks of the distribution functions but much less affects the coordination numbers of the solvation shells.On the other hand, while using the HNC approximation for the repulsive core range, the PLHNC closure (21), in fact, switches in the enrichment regions of h > 0 to the mean spherical approximation (MSA).For a system of intermediate density close to a critical region, it is an entire tale of the radial distribution function which is long-range and non-oscillating in such a case.This results in better description provided by the PLHNC closure for a liquid near critical and phase transition regions, in contrast to the pathological predictions of the HNC approximation for the stability of fluid phases and the liquid-vapour phase diagram [53].

Thermodynamics
The expressions derived within the replica OZ theory for the thermodynamics of quenched-annealed atomic liquids [11][12][13][14]16] can all be extended to the replica RISM approach at hand.In particular, a chemical potential is of primary importance to calculate isotherms of adsorption in a disordered porous matrix.Notice that as has been introduced by Given [13], several natural definitions of the chemical potential of the mobile liquid sorbed in the matrix are possible in terms of a grand canonical Monte Carlo (GCMC) simulation which are in general different quantities.We will utilize the one following from the chemical potential of the replicated system in the limit s → 0 and obeying the standard thermodynamic relations for the sorbed liquid in equilibrium with its bulk phase outside the matrix.The averaged free energy of the mobile liquid defined by equation ( 10) is determined by the canonical ensemble variables (V, T, N 1 ) and the parameters (ρ 0 , T 0 ) specifying the matrix state.By using (11), its full differential, dA 1 (V, T, N 1 ; ρ 0 , T 0 ), is readily expressed through the corresponding change in the free energy (12) of the fully annealed replicated (s + 1)-component mixture with the same number of molecules in all replicas of the liquid, where insertion of a molecule of the liquid in the quenched-annealed system turns upon replication of the liquid into simultaneous insertion of s molecules, one in each replica.The chemical potential in (22) is defined by the usual derivative for the annealed (s + 1)-component mixture at equal numbers of molecules in the replicas but with no constraint for their independent change, The value S rep (s) in (22) has the standard sense of entropy of the common (s + 1)component mixture.The scaling factor λ 0 = T /T 0 in (8) can be regarded as an independent coupling parameter for the matrix-matrix interactions.Then, "switching" the latter by the infinitesimal change dλ 0 contributes additionally to the free energy differential dA rep (s).The corresponding partial λ 0 -derivative obtained by differentiation of the free energy ( 12) is expressed through the thermal average U 00,rep (s) of the matrix-matrix interaction potentials in the annealed replicated mixture with the canonical partition function ( 8), This brings about the last term in (22a).On the other hand, the differential dλ 0 is related to the infinitesimal change of the thermodynamic parameters by the expression dλ 0 = dT /T 0 − T dT 0 /T 2 0 , which casts the free energy change into the form (22b). Applying (11) to (22) yields the infinitesimal change in the averaged free energy of the mobile liquid in the quenched-annealed system as with the following relations for the partial derivatives, where the thermodynamic quantities P 1 , S 1 , µ 1 , ε 1 are defined, following [12], at constant matrix density ρ 0 rather than number of matrix molecules N 0 .This is most appropriate to a quenched-annealed system since the quenched porous material does not stretch with the volume of the sorbed mobile liquid.The second term in the entropy (27) and the presence of the matrix-temperature derivative (30) of the averaged free energy of the liquid are peculiar to a quenched-annealed system with the matrix distributions corresponding to soft core potentials u 00 αγ between matrix particles.In the case of hard core potentials u 00 αγ , these terms vanish due to the zeroth internal energy of the matrix component, U 00,rep (s) = 0, and the entropy reduces to the form derived in [12].This can also be seen as cancellation of the two terms when T 0 is replaced with T for hard core matrices.
Introducing the internal energy of each replica of the liquid in the annealed (s + 1)-component mixture as and applying (11) to the thermodynamic relation A rep (s) = U rep (s) − T S rep (s) for the replicated system with the internal energy U rep (s) = U 00,rep (s) + sU 1,rep (s), we get that it holds for the liquid sorbed in the quenched matrix, where Also, the Gibbs-Duhem equation for the annealed replicated mixture, leads to that for the quenched-annealed system, which together with equation ( 25) allows one, on integrating dA 1 , to recast the averaged free energy (32) into the other standard form For the annealed RISM equations ( 13) complemented with the HNC closure ( 19)- (20), the chemical potentials in excess of the ideal part are readily obtained by extension of Singer and Chandler's closed analytical form [48] to the liquid mixture.With allowance for the symmetry of the s replicas of the liquid, they are written as ∆µ 1,rep (s) = ∆µ 10,rep (s) + ∆µ 11,rep (s) + (s − 1)∆µ 12,rep (s), (37) ∆µ 0,rep (s) = ∆µ 00,rep (s) + s∆µ 01,rep (s) (38) with the partial contribution ∆µ ij,rep (s) of mixture component j written as the familiar HNC expression Inserting ( 37) into (28) and passing to the replica identity limit s → 0, the excess chemical potential of the liquid sorbed in the quenched matrix takes the form where the blocking contribution, ∆µ (b) , is calculated from expression (39) with the blocking site-site correlations, h αγ (r) and c (b) αγ (r).Notice that the excess value (39), and thus (40) result from intermolecular site-site potentials, and are defined over the ideal part of rigid noninteracting molecules, unlike the ideal gas of separate interaction sites used in equation (10).Also, equation ( 29 In the case of the PLHNC approximation ( 21), the excess chemical potential has a closed analytical form as well [45].Instead of (39), the expression to be inserted into equations (40) and ( 41) is then where Θ(x) is the Heaviside step function.We have presented the derivation of the PLHNC form (42) in the appendix of [45].It differs from the HNC one (39) in the term h 2 /2 which is switched off in the regions of the density profile enrichment, h > 0, where the partial linearization is applied.It can be also shown by using the procedure of [45] that in the case of the mixed closure combining the HNC approximation (20b) for the blocking part with the PLHNC one (21a) for all other correlations, the excess chemical potentials are determined by equations ( 40) and ( 41) with the blocking contribution ∆µ (b) in the HNC form (39) and all the other terms given by the PLHNC expression (42).
The pressure (26), entropy (27), and the matrix-state derivatives ( 29) and (30) for the quenched-annealed system cannot be obtained from those of the replicated mixture in a straightforward manner.Rosinberg, Tarjus, and Stell [12] suggested the way to estimate them by using finite differences, for instance, lim s→0 d ds ∆µ 00,rep (s) = ∆µ 00,rep (s = 1) − ∆µ 00,rep with the contributions of order ρ 0 (ρ 1 ) 2 to be neglected.However, there emerges an obstacle in solving the replicated RISM equations ( 13) at s = 1 because of high nonadditivity of the replicated fully annealed mixture.Meroni, Levesque, and Weis [15] have found for atomic liquids that the replicated OZ equations at s = 1 become divergent for higher liquid densities, even though they have the solution in the limit s → 0. As has been noted in [15], such systems resemble by their strong nonadditivity a Widom-Rowlingson mixture known to phase separate [54].On the other hand, the pressure and the matrix-state derivatives for the quenched-annealed system can be expressed [12,13] in terms of the s-derivatives of the correlation functions at s = 0. Instead of (43), Given [13] proposed to obtain the latter from the linear integral equation derived in the limit s → 0 by differentiation of the matrix-matrix OZ equation for the replicated annealed mixture, together with its closure.Along similar lines, for the HNC-type form of the excess chemical potential, equation (39), we write its s-derivative as where the s-derivatives of the matrix correlations are determined from the integral equations derived by differentiating the matrix-matrix RISM equation (13a) with the HNC closure (20a) and then passing to the limit s → 0, dh provided the matrix correlations h 00 αγ and c 00 αγ have been determined before.For the PLHNC approximation (42), the s-derivative of the excess chemical potential is accordingly modified to with the s-derivatives of the correlation functions calculated from the integral equation (45) complemented with the s-derivative of the PLHNC closure, Similarly, the derivative (30) of the free energy of the mobile liquid with respect to the matrix temperature can be obtained from the matrix internal energy (24) as with the s-derivative dh 00 αγ /ds determined from equations (45), and ( 46) or ( 48) in either approximation.
Finally, it is easy to show in the grand canonical ensemble similarly to [12] that the isothermal compressibility of the mobile liquid in the matrix, χ 1 , is related to the fluctuation of the average number of particles as where the brackets and overbar denote respectively the thermal and disorder averages, and N 1 = N 1 .The normalizations in the grand canonical ensemble of the single and pair site densities at a particular matrix realization are written as On averaging them over matrix disorder, the compressibility expresses through the connected part of the site-site total correlation function of the liquid for any pair of sites αγ, By using the replica RISM integral equation (14e) for the connected terms of the site-site correlations in reciprocal space at k → 0, the compressibility (52) can also be recast as

Models and parameters
In order to illustrate the replica RISM theory, we calculated the structure and some thermodynamic properties of water in a quenched disordered matrix roughly modelling microporous carboneous material.We use two matrix models: the first consisting of a quenched liquid of Yukawa spheres associating into clusters to roughly represent carboneous microporous material, and the second including also activating carboxylic groups -COOH.For these systems, we solved the replica RISM equations ( 14) complemented with the PLHNC approximation (21a), more appropriate for calculation of the thermodynamic properties as well as robust for convergence at peculiar thermodynamic conditions.An exception is the closure for the blocking correlations, for which we employ the HNC approximation (20b).Since the blocking part of the total correlation functions, h (b) αγ (r), is positive almost everywhere(see below), the PLHNC closure for the blocking correlations practically reduces to the MSA.The latter, however, neglects the blocking direct correlations functions and thus would reduce the replica equations to the Madden-Glandt approach.The below results of the replica RISM theory are in agreement with the qualitative conclusions found in experiment for the structure of water sorbed in porous solids [57].A comparison with simulations which is, no question, very important to judge about the quantitative precision of the theory proposed will be made in further studies.
The replica RISM/PLHNC-HNC integral equations have been solved on the radial linearly spaced grid of 16384 points with grid resolution of 0.01 Å.The equations have been converged to an accuracy of the root mean square residual of 10 −13 by using the modified direct inversion in the iterative subspace (MDIIS) method elaborated for liquid state theory calculations by Kovalenko, Ten-no, and Hirata [55].Details of the procedure, construction of the residue of the integral equations, and the choice of initial vectors of the correlations have been presented in [46,47,55,56].Notice that we calculate the convolutions in equations ( 14) by employing the linear fast Fourier transform (FFT) technique [58], but not the exponentially spaced nonlinear FFT [59,60] common for the RISM treatment of polar and ionic systems [33][34][35][45][46][47]55,60].Although numerically advantageous and fast in the treatment of long-range electrostatic interactions, the nonlinear FFT suffers from ghost oscillations near the direct and reciprocal grid origin which can become huge because of accumulation of errors in the nonlinear scaling transformations.Such an artifact especially affects the blocking correlations, essential in the repulsive core region and hampers calculation of the thermodynamic properties related to the values of the correlation functions at k = 0.
The potential between matrix particles in the first model is represented by the sum of the repulsive and attractive Yukawa terms, u 00 CC (r) = 4ǫ 00 CC exp 2 Its zero is located at separation d 00 CC , and the minimum of depth ǫ 00 CC is at l 00 as = 1 + δ 00 CC d 00 CC ln 2. For the size, well width and depth parameters chosen as d 00 CC = 2.38 Å, δ 00 CC = 0.119 Å, and ǫ 00 CC = 6.0 kcal/mol, respectively, the potential (54) decays rapidly with particle separation, and thus has a narrow and deep attractive well with the minimum at the distance l 00 as = 2.46 Å corresponding to the separation between adjacent bonds connecting two graphite basal planes.The association energy parameter ǫ 00 CC and the matrix density ρ 0 C are so chosen as to provide formation of interconnected branched chain clusters of matrix "carbon" atoms.The number density of matrix atoms in the first model is taken to be ρ 0 C = 0.010841 Å−3 .This results in matrix porosity p = 0.80 (see below).In the second model, activating carboxylic groups -COOH in molar ratio 1:24 are rigidly grafted to the corresponding fraction of matrix atoms.To keep same porosity, the number density of matrix "carbons" in the second model is somewhat lower, ρ 0 C = 0.009946 Å−3 .Accordingly, their association well is made somewhat deeper: ǫ 00 CC = 6.08 kcal/mol.For both matrix models, the number density of confined water molecules is set as the p-th part of the ambient water density: ρ 1 = 0.03000 Å−3 .
Except for the two-Yukawa potential (54) between matrix particles associating into "carboneous" matrix material, all other site-site interactions in the system are described by the common potential comprising the Coulomb and 12-6 Lennard-Jones (LJ) terms, where the LJ diameter and energy parameters for cross-terms are determined from those of interaction sites by the standard Lorentz-Berthelot mixing rules, The site charges and LJ parameters for both liquid and matrix species are listed in table 1. Water is represented by the extended simple point charge (SPC/E) model of Berendsen, Grigera, and Straatsma [61].The only difference is that a LJ size of 0.4 Å is introduced for the water hydrogen sites, following Pettitt and Rossky [62].This is an adjustable parameter of the RISM theory which does not affect the entire potential of a water molecule since the hydrogens are situated well inside the oxygen core, however allows one to optimize the description for hydrogen bonds [62].Same values are assigned to the LJ size and energy of the hydrogen site entering into the OH chemical bond in -COOH activating groups of the matrix in the second model.The LJ parameters of "carbon" particles constituting the matrix in the two models are the same as for carbon atoms of microporous carbon solids in simulations [63].The geometry of the -COOH group is taken as that [64] of acetic acid (CH 3 COOH), with the methyl group replaced by the matrix "carbon" atom to which the activating group is grafted.Their site charges and LJ parameters are taken from the simulation [63]; a small neutralizing charge is assigned to a grafting "carbon" atom of the matrix.

Matrix structure
We complement the matrix-matrix integral equation (14a) with the PLHNC closure (21a) which enables, as has been noted above, the description of matrix "carbon" atoms associating into interconnected branched chain clusters.Figure 1 demonstrates this in terms of the radial distribution function between matrix "carbons", g 00 CC (r), as well as of the running excess coordination number, ∆N 00 CC (r).The latter is defined in general for a molecular mixture as the number of sites γ of component b of species j (j = 0, 1) in excess of their average density, around site α of the labelled molecule of sort ai, As determined by the potential well minimum, atoms strongly associate with the maximal probability g 00 CC (r) at separation l 00 as = 2.46 Å same as that between carbons on the graphite surface (solid lines in figure 1a).In this first peak of the distribution, the excess coordination number ∆N 00 CC (r) gains a value close to 2, which corresponds to association of two neighbours in rolling contact with the labelled atom (solid lines in figure 1b).Then follows a plateau of height about 2 on the distribution g 00 CC (r).Two cusps of ∆N 00 CC (r) at r 1 ≈ 2.7 Å and r 2 ≈ 5.1 Å mark the limits of the second coordination shell.The excess number of atoms in the first and second shells together amounts to ∆N 00 CC (r 2 ) ≈ 4.This can be interpreted as two atoms in the second coordination shell that are associating each with one of the first-shell neighbours so that a bended chain forms.The flat distribution g 00 CC (r) and the according cubic behaviour of ∆N 00 CC (r) in the second coordination shell are determined by the angular correlations of chain bending.The third shell cusp of ∆N 00 CC (r) at r ≈ 7.5 Å is barely visible.The average number of atoms in a cluster, estimated in terms of the static structure factor S 00 CC (k = 0) = 1 + ∆N 00 CC (r → ∞) amounts to about 22 "carbons" (figure 1c).Saturation in the cluster size is gained for the most part at r ≈ 30 Å.The linear slope of ∆N 00 CC (r) within this range is also indicative of the chain-like character of the cluster (for a three-dimensional aggregate, it obviously would increase as r 3 ).The chain bending can be roughly estimated from the slope value of 3.3 particles per association length l 00 as (in both chain directions).Compared to 2 particles per l 00 as in a fully extended conformation, this suggests a helical-like chain conformation with a helicity angle of about 15 degrees.The matrix porosity p defined as the volume fraction not occupied by the matrix particles is estimated for this matrix by assuming the volume of a "carbon" atom to be that of a sphere with diameter equal to its LJ size σ 0 C , and taking into account overlap of chain neighbours associating at the length l 00 as .For the above matrix density and sizes, this yields the value of p = 0.80.
The second matrix model also includes activating -COOH groups in molar concentrational ratio ρ 0 COOH /ρ 0 C = 1 : 24, grafted as described in the previous section to the corresponding part of matrix "carbon" atoms.The matrix density is adjusted accordingly to keep the porosity value of p = 0.80.The mean density of -COOH groups per surface area of matrix chains is about 0.0035 Å−2 , which is of the order of typical surface density of carboxylic groups in activated carbons as studied in experiment [65] and simulations [63].Inclusion of activating groups induces neutralizing charges on matrix grafting "carbons", which brings about additional repulsion between them.This does not change the height of the first association maximum (dashed line in figure 1a), but somewhat rises the interatomic separation in a chain and increases the cluster extent up to 40 Å (dashed lines in figures 1b and 1c).The qualitative picture of matrix association as well as of the "carbon" branched chain conformation remains similar.As expected, the correlations between activating groups are smooth and do not exhibit any clustering behaviour.For example, the pair distribution function of -COOH carbons possess a moderate wide maximum at the separation determined by their LJ radius as well as mutual electrostatic repulsion (dash-dotted line in figure 1a).The running excess coordination number (dash-dotted lines in figures 1b and 1c) saturates at a small value of about 2.9 which is appropriate to the specified concentrational ratio of -COOH groups.

Structure and thermodynamics of the sorbed water
First consider the site-site distribution functions of water, g 11 αγ (r). Figure 2 exhibits the changes in the structure from bulk ambient water (solid lines) upon sorption into the matrix (dashed lines).For all the distributions, evident is considerable enhancement of the short-range structure peaks, in particular, the hydrogenbonding peak of the oxygen-hydrogen distribution g 11 OH (r).The curves, especially the hydrogen-hydrogen distribution g 11 HH (r) are noticeably shifted upwards in the region of the first and second hydration shell.Moreover, the local peak of the oxygen-oxygen distribution g 11 OO (r) at separation r ≈ 4.4 Å corresponding to the "fingerprint" of the tetrahedral bonding structure of bulk ambient water is still present in the structure of sorbed water.This is in good qualitative agreement with neutron scattering experiments [57] revealing that water confined in micro-and mesoporous matrices exhibits enhanced hydrogen bonding, with an increasing effect for lower temperatures and smaller pore sizes.It is also worthwhile noting that the above increase of the site-site distribution peaks of sorbed water resembles the enhancement of the water correlations in water-alcohol binary mixtures at relatively small alcohol molar fractions, observed in molecular simulation [66] and experiment [67].
Since such thermodynamic properties as compressibility of the sorbed liquid are determined by equations ( 52) and (53) in terms of the connected parts of the sitesite correlations, the connected and blocking terms are of as much interest as the full correlations.Radial distributions functions between water interaction sites, following from the replica RISM/PLHNC-HNC theory.The oxygen-oxygen, oxygenhydrogen, and hydrogen-hydrogen distributions (parts (a) to (c), respectively).The full distribution functions, g 11 αγ (r), of bulk ambient water (solid lines) and water sorbed in the matrix of porosity 80% (short-dash lines).The blocking part of the site-site total correlation functions, h (b) αγ (r), following from the replica RISM/PLHNC-HNC theory (long-dash lines), and from that with the Madden-Glandt approximation neglecting the blocking part of the direct correlation functions (dash-dotted lines).HNC theory (long-dash lines).For comparison, we present the blocking part h αγ (r) obtained by the decomposition (14f) from the replica RISM integral equations ( 14) reduced to the Madden-Glandt approximation (dash-dotted lines) by neglecting the blocking part of the site-site direct correlation functions, c (b) αγ = 0 in equation (14g), and complementing the remaining equations with the PLHNC closures (21a).It is apparent that as compared to the MG approximation, the proper treatment provided by the replica RISM approach yields an essential part of the blocking effect, in particular, the high zero-separation values h (b) αγ (r = 0) and the oscillations in the range of the core and the first solvation shell.
The upwards shift of the site-site distributions of sorbed water, g 11 αγ (r), results in the increase of the site-site coordination numbers.Figure 3 presents the running excess coordination numbers for sorbed water, ∆N 11  αγ (r), defined by equation (56).It accounts for the density of each of the two hydrogen sites of a water molecule separately, which is convenient for comparison between the limiting behaviour for different site pairs.According to the decomposition (14f), the excess coordination numbers are broken up into the connected and blocking terms ∆N αγ (r) (dashed lines in upper and lower panels, respectively).The connected part deviates little from ∆N 11 αγ (r) for bulk ambient water (solid lines).In the latter case the blocking part is zero.In sorbed water, the oscillations of ∆N OO (r) (insertion in figure 3b).All they possess a long-range asymptotics saturating at the distance of 30 to 40 Å, which results from the long-range nature of the matrix-matrix correlations manifesting in the matrix excess coordination number, ∆N 00 CC (r) (figure 1).The blocking contribution to the water coordination numbers reaches a value of ∆N (b) αγ (r → ∞) ≈ 7 excess water molecules.Table 2 exhibits the water compressibility calculated from equations (52) or (53).In the former, it is proportional to the mean fluctuation of the connected density,

+ ∆N
(c) αγ (r → ∞), that is a matrix-disorder average of spatially inhomogeneous clusters formed in the local cavities of the matrix, according to the definition of the connected density correlation (17).It comprises the labelled water molecule, and ∆N (c) αγ (∞) excess water molecules surrounding it in a connected density fluctuation.In ambient water, the latter quantity is negative and close to unity (figure 3), resulting in a subtle balance for the fluctuation size and thus for the compressibility.The RISM/PLHNC treatment for bulk ambient water yields the compressibility value in pretty good agreement with experiment [64].The 20% decrease of the water density ρ 1 in the matrix leads to the corresponding increase of the compressibility.However, it is almost entirely cancelled with the increased local clustering, h (c) αγ (k = 0), so that the limiting value ∆N (c) αγ (∞) for water sorbed in the matrix is lower than for bulk water by less than 1%.This tiny decrease, however, leads to a significantly larger drop in the magnitude of the mean fluctuations (ρ 1 kT χ 1 ) in the sorbed water.As a result, the compressibility χ 1 for the sorbed water increases against bulk ambient water by 8%.For comparison, we also present the compressibility value following from the Madden-Glandt RISM/PLHNC approximation.The latter predicts a considerably larger increase of the compressibility for the sorbed water.
The matrix effect on the long-range structure of the liquid can be observed alternatively in its partial static structure factors.The intermolecular terms of the partial structure factors for the quenched-annealed molecular system are defined by the usual relation S ij αγ (k) = ρ tot x i x j h ij αγ (k), where ρ tot = M 0 ρ 0 + M 1 ρ 1 is the total number density of atoms in the system, and x i = ρ i /ρ tot is the molar fraction of each site of species i.The partial structure factors of the liquid, S 11  αγ , split up according to (14f) into the connected and blocking terms S αγ .Figure 4 depicts the partial structure factors of water confined in the matrix versus those of bulk water, calculated by means of the replica RISM/PLHNC-HNC approach.Again, the connected parts of the structure factors of the bulk and sorbed water (solid and short-dash lines) are rather close, with the oscillations of the same phase but smaller amplitude in the latter case.Notice that despite the small difference between ∆N (c) αγ (∞) for the sorbed and bulk water (figure 3), the values of S αγ (r).Now proceed to the water-matrix structure.Figure 5 depicts the radial distribution functions between water sites and "carbons" of the non-activated matrix of the first model (dashed lines) as well as matrix atoms at infinite dilution (dashed lines).The water-"carbon" atom distributions at ρ 0 = 0 are, in fact, those between water and a hydrophobic solute, observed in simulations and obtained by using the RISM theory for hydrophobic hydration (see, for instance, literature in [68]).In the matrix of porosity 80%, the relative position and shape of the oxygen-and hydrogen-"carbon" distribution functions remain very similar, except for the oscillations noticeably reduced.On the other hand, the water-matrix distribution functions   acquire a long-range depletion tail.This is also clearly seen from the behaviour of the running excess coordination numbers of water around a matrix "carbon" atom, ∆N 01 CO (r) and ∆N 01 CH (r), shown in figure 6.A single matrix atom induces the depletion extending to about two hydration shells, with the lack of somewhat more than one water molecule on average.In the matrix, in contrast, although their oscillations fade on two hydration shells as well, the excess coordination numbers ∆N 01 Cα (r) linearly fall down to saturate at the level of about 24 water molecules expelled from the depletion region of radius about 30 to 40 Å.Both the linear character and range of decay of ∆N 01 Cα (r) are physically reasonable since they follow from those of the running excess coordination number ∆N 00 CC (r) characterizing the matrix structure.Finally, consider the case of the matrix of the second model including activating carboxylic groups.We do not present the site-site correlations of water sorbed in the activated matrix, as they are very similar to those in the matrix without -COOH groups.The compressibility of sorbed water, however, is much more sensitive to the presence of activating groups.They reduce the compressibility increase relative to the bulk water value by about 20%, as is seen in table 2. Notice also that for the activated matrix, the compressibility value following from the MG-RISM/PLHNC approximation decreases as well, continuing to be proportionally higher than the result of the replica RISM/PLHNC-HNC approach.
The decrease of the water compressibility in the activated matrix (table 2) could be attributed in particular to the specific adsorption of water on activating molecules.Figure 7 exhibits the radial distributions between the interaction sites of water and -COOH molecules.Their short-range structure shows formation of strong hydrogen bonding between water and carboxylic molecules: the -COOH hydrogen bound to the water oxygen (first peak in part d), and the water hydrogens bound to both the oxygens of -COOH (first peak in parts f and g).The strongest hydrogen bonding is to the water oxygen due to its negative charge larger than that of the -COOH oxygens.Although the two -COOH oxygens possess almost the same negative charge, the water hydrogen bonding to the -COOH oxygen of its OH group is weaker because of the screening by the hydrogen positive charge and also because of the steric constraints imposed by the matrix "carbon" chain the -COOH group is grafted to.The other peaks of the distributions g 01 O=,H (r) and g 01 OH (r) do not differ much.The distributions between the water and -COOH oxygens, g 01 O=,O (r) and g 01 OO (r), are very similar too (parts b and c).They resemble the oxygen-oxygen distribution of bulk ambient water, including the peak at 4.6 to 5 Å corresponding to the "fingerprint" of the hydrogen-bonding tetrahedral structure.The distribution between the -COOH and water hydrogens, g 01 HH (r) (part h), is also similar to the hydrogen-hydrogen distribution of water, g 11 HH (r) (figure 2).The first peak of the distribution function between the -COOH carbon and the water oxygen, g 01 CO (r) (figure 7a), is formed to a large extent by all the three arrangements of the water and -COOH molecules making the hydrogen bond.In the carbon-hydrogen distribution, g 01 CH (r) (figure 7e), first is a cusp at r ≈ 2.7 Å on the ascending slope which results from the first peak of the water hydrogen bonding to the -COOH oxygens (largely to O=), smeared due to bonding at different angles.Next is the first maximum corresponding to water hydrogens   Much as in the matrix without activating groups, all the matrix-liquid distribution functions in the activated matrix acquire the same long-range depletion asymptotics, as compared to the limit of vanishing matrix density ρ 0 aα → 0 which reduces to the conventional solute-solvent system at infinite dilution (respectively, dashed and solid lines in figure 7).The distribution functions between "carbon" atoms constituting matrix chains and water oxygen and hydrogen are very similar to those in the non-activated matrix (figure 5).In contrast to the latter, however, the shortrange structure reveals considerable enhancement of the hydrogen bonding between activating carboxylic groups and water molecules as they are sorbed in the matrix.It manifests in the rise of the first hydration shell peaks for almost all the site-site distribution functions between -COOH groups and water, in spite of their long-range depletion.
Table 3 shows the excess chemical potential of sorbed water, µ 1 , calculated by equation (40).A comparison of its constituents reveals that the main change with respect to µ 1 in bulk ambient water comes from the liquid-matrix component, µ 10 .For the non-activated matrix, the liquid-liquid term µ 11 gives an 11% contribution to the change.The blocking term µ (b) amounts to a small value of about 1% for the present system.As activating carboxylic groups are added to the matrix, the term µ 10 somewhat decreases, apparently because of the additional hydrogen bonding between -COOH and water molecules.The blocking term µ (b) remains almost the same.However, the change in the liquid-liquid term µ 11 rises almost twice, which results in noticeable increase of the excess chemical potential with the activation of the matrix.

Conclusion
Based on the replica formalism in integral equation theory for a quenchedannealed system, we have developed a replica RISM integral equation theory for polar molecular liquids sorbed in a quenched disordered microporous material.The description is readily extended to the cases of the mobile liquid comprising a mixture of ion-molecular components as well as of the disordered matrix including ionic and polar species.As compared to the RISM equations in the Madden-Glandt approximation for a molecular quenched-annealed system, the replica RISM approach provides an advanced description of the so-called blocking effects, that is of correlations between mobile liquid particles due to the presence of disordered obstacles frozen in space.The blocking effects essentially distinguish a quenched-annealed system from a fully annealed mixture, and are especially important for systems including ionic and polar molecular species.The replica RISM approach also improves the results for the thermodynamics of the sorbed liquid.
We complement the replica RISM integral equations by the HNC closure as well as its partial linearization (PLHNC), adequate to the case of ionic and polar molecular liquids.As a merit, the PLHNC approximation ensures the existence and smooth convergence of the solutions, and on the other hand, provides better description for liquids near critical and phase transition regions.The HNC closure is employed for the blocking site-site correlations.The use of the HNC and PLHNC approximations allowed us to derive a closed analytical form for the excess chemical potential of a mobile molecular liquid sorbed in a molecular matrix.Further improvement of the closures would imply bridge corrections to enforce thermodynamic and structural consistency, similarly to those known for atomic quenched-annealed systems.
A common description of a quenched-annealed system in the replica integral equation approach is restricted to either of the components comprising hard-core particles, in order to circumvent the presence of two different temperatures describing the equilibrium spatial distributions of the matrix and mobile liquid.To consider more realistic systems, we have extended the thermodynamic description of the mobile liquid to the case of soft core interaction potentials between all species of the quenched-annealed system, in which the liquid and matrix equilibrium distributions are characterized in general by two different temperatures.The pressure and the matrix-state derivatives for the sorbed liquid can be expressed in terms of the s-derivatives of the correlation functions with their analytic continuation at the number of liquid replicas s = 0. We obtain the latter from the linear integral equations derived by differentiation of the replica RISM integral equations together with their closures.
We illustrate the replica RISM/PLHNC-HNC approach by the calculation for water sorbed in a matrix of "carboneous" material with porosity 80%, modeled as chain-like clusters of "carbon" atoms.In the other model, the matrix is activated by carboxylic groups (-COOH) grafted to matrix atoms.The structural and thermodynamic properties of the sorbed water are compared with bulk ambient water, and discussed in view of experimental evidences for water in confined geometry.The matrix confinement considerably enhances the hydrogen bonding of water and the adsorption of water at activating -COOH groups.We postpone a detailed analysis of activation effects as well as the comparison of the results of the replica RISM approach with simulations to further work.
In the description of quenched-annealed systems, modelling the disordered porous material is of paramount importance and, in fact, determines salient physico-chemical properties of the system.A widely used approach consists in representing the disordered porous matrix by means of the random or equilibrium distribution of quenched mesoscopic-size spheres.A more complicated model of porous carbon com-prises randomly distributed graphitic platelets [69].We feel, however, that rather than to introduce a porous matrix model ad hoc, it is more appropriate to involve a microscopic description predicting the structure of pores in disordered material from the real mechanism used in its formation.In the present work, we employ a simple model describing the formation of a porous material consisting of chain-like clusters of "carbon" atoms.Matrix formation by gelation process described within Wertheim's association theory [70] can be beneficially used instead.A new approach to obtain porous materials employs a template removed after the synthesis to form a material with controlled porosity and a pore architecture specific for a given purpose (see literature in [23]).In exploration of such systems of practical interest, Zhang and Van Tassel [23] elaborated very recently the replica OZ theory for adsorption of an atomic liquid in a templated matrix that is formed as a binary atomic mixture with the template particles removed after quenching.An extension of the replica RISM theory to this approach would be of great practical as well as theoretical interest for a description of realistic molecular models of liquids adsorbed in templated porous materials.
) defining the change in the averaged free energy of the sorbed liquid upon adding a molecule to the quenched matrix of constant volume and temperature leads to ∆ν 1 = lim s→0 d ds ∆µ 00,rep (s) + ∆µ 01 .

Figure 1 .
Figure 1.Spatial correlations in the quenched matrix of "carbon" atoms associating into interconnected branched chains.The radial distribution functions g 00 CC (r) (part (a)), and the running excess coordination numbers ∆N 00 CC (r) in the first three coordination shells (part b) and at a distance (partc).Correlations between matrix "carbons" in the pure matrix (solid lines), and in the presence of activating -COOH groups (dashed lines).Intermolecular correlations between carbon atoms of -COOH groups (dash-dotted lines).
Figure 2 also shows the blocking part of the site-site correlation functions of the sorbed water, h (b) αγ (r), following from the replica RISM/PLHNC-

Figure 2 .
Figure 2.Radial distributions functions between water interaction sites, following from the replica RISM/PLHNC-HNC theory.The oxygen-oxygen, oxygenhydrogen, and hydrogen-hydrogen distributions (parts (a) to (c), respectively).The full distribution functions, g 11 αγ (r), of bulk ambient water (solid lines) and water sorbed in the matrix of porosity 80% (short-dash lines).The blocking part of the site-site total correlation functions, h αγ (r) are of a somewhat smaller amplitude, but their phasing remains very similar to∆N 11  αγ (r) in bulk water.The main part of the rise in the coordination numbers of water upon sorption into the matrix appears in their blocking part.The blocking terms ∆N are very close to each other, except for a small difference in the form of tiny oscillations in the first hydration shell, more pronounced for ∆N(b)

Figure 3 .
Figure 3. Water oxygen-oxygen (O-O), oxygen-hydrogen (O-H), and hydrogenhydrogen (H-H) running excess coordination numbers following from the replica RISM/PLHNC-HNC theory.The connected and blocking contributions, ∆N αγ (r) (parts (a) and (b), respectively).Water sorbed in the matrix shortdash lines), and bulk ambient water (solid lines).The coordination numbers are treated separately for each of the two water hydrogen sites.
αγ (k = 0) are noticeably distinct due to the decreased numerical density of water sorbed in the matrix.The blocking parts S (b) αγ (k) of the sorbed water exhibit strong enhancement for k → 0 and decay rapidly beyond the first maximum, in accord with the smooth and long-range increase of the corresponding blocking excess coordination numbers ∆N (b)

Figure 4 .
Figure 4. Water oxygen-oxygen (O-O), oxygen-hydrogen (O-H), and hydrogenhydrogen (H-H) intermolecular partial structure factors following from the replica RISM/PLHNC-HNC theory.The connected and blocking contributions for water sorbed in the matrix, S αγ (k) (short-and long-dash lines, respectively).The intermolecular partial structure factors of bulk ambient water (solid lines).The partial structure factors are treated separately for each of the two water hydrogen sites.

Figure 5 .
Figure 5. Radial distribution functions between the matrix "carbon" atoms, and oxygen and hydrogen of sorbed water (upper and lower panels, respectively), following from the replica RISM/PLHNC-HNC theory.Matrix of porosity 80% and of vanishing density ρ 0 → 0 (dashed and solid lines, respectively).

Figure 6 .
Figure 6.Excess coordination numbers of water oxygen and hydrogen (solid and dashed lines, respectively) around a matrix "carbon" atom, following from the replica RISM/PLHNC-HNC theory.Short-and long-range structure (upper and lower panels, respectively).Matrix of porosity 80% and of vanishing density ρ 0 → 0 (lower and upper curves, respectively).

Figure 7 .
Figure 7. Radial distribution functions between the interaction sites of matrix activating carboxylic groups (-COOH) and sorbed water, following from the replica RISM/PLHNC-HNC theory.Activated matrix of porosity 80% and of vanishing density ρ 0 aα → 0 (dashed and solid lines, respectively).
to this (s + 1)-component molecular mixture.With allowance for the symmetry of the s liquid replicas, they take the form

Table 1 .
Charges and Lennard-Jones parameters of the interaction sites of liquid water and quenched matrix material.

Table 2 .
Compressibility of water sorbed in the activated as well as non-activated matrix against that of bulk ambient water.Results of the replica RISM and MG-RISM approaches.

Table 3 .
(40)ss chemical potential, µ 1 , and its components in equation(40)for water sorbed in the activated as well as non-activated matrix, and for bulk ambient water.Results of the replica RISM/PLHNC-HNC approach.